From 4b79df2b9c8c08bf1e6248683eb09bd4a1fa2044 Mon Sep 17 00:00:00 2001 From: Sifu Luo Date: Fri, 2 Feb 2018 14:51:13 -0600 Subject: [PATCH 01/10] clip at 8 for low pt --- JetAnalyzers/bin/compareJEC.C | 199 +++++++++++++------------ JetAnalyzers/bin/jet_synchtest_x.cc | 22 +-- JetAnalyzers/test/run_JRA_cfg.py | 10 +- JetUtilities/bin/BuildFile.xml | 2 - JetUtilities/interface/ClosureMaker.hh | 2 +- JetUtilities/src/ClosureMaker.cc | 27 ++-- JetUtilities/src/L2Creator.cc | 11 +- 7 files changed, 140 insertions(+), 133 deletions(-) diff --git a/JetAnalyzers/bin/compareJEC.C b/JetAnalyzers/bin/compareJEC.C index c8cbacf4..f6c6e478 100644 --- a/JetAnalyzers/bin/compareJEC.C +++ b/JetAnalyzers/bin/compareJEC.C @@ -67,10 +67,10 @@ double getRho(double mu) { if(TString(_payld).Contains("RunISummer12") || TString(_payld).Contains("Winter14_") ) // Run1 53X MC {p[0]=-0.426 + ue; p[1]=0.504; p[2]=0.0005;} - if(TString(_payld).Contains("Run1Sub53X")) - {p[0]=-0.350 + ue; p[1]=0.526; p[2]=0.001;} + if(TString(_payld).Contains("Run1Sub53X")) + {p[0]=-0.350 + ue; p[1]=0.526; p[2]=0.001;} - if(TString(_payld).Contains("Run1Sub742")) + if(TString(_payld).Contains("Run1Sub742")) {p[0]=-0.855 + ue; p[1]=0.461; p[2]=-0.001;} if(TString(_payld).Contains("50nsRunIISpring15") || TString(_payld).Contains("Summer15") || TString(_payld).Contains("RunIISpring15DR74_bx50") ) @@ -88,8 +88,8 @@ double getRho(double mu) { if( TString(_payld).Contains("Fall15_25ns") ) // {p[0]=0.404 + ue; p[1]=0.374; p[2]=0.0117;} ?? // {p[0]=-1.036 + ue; p[1]=0.705; p[2]=-0.0016;} ?? - - {ue=0.; p[0]=1.084 + ue; p[1]=0.6075; p[2]=0.;} // Alexx; + + {ue=0.; p[0]=1.084 + ue; p[1]=0.6075; p[2]=0.;} // Alexx; // if( TString(_payld).Contains("Summer15_25nsV5ch2") ) // {p[0]=-1.036 + ue; p[1]=0.705; p[2]=-0.0016;} @@ -120,9 +120,9 @@ double getNPV(double mu) { // ------------------------------------------------------------ void setEtaPtE(FactorizedJetCorrector *jec, double eta, double pt, double e, int mu) { - + assert(jec); - + int npv = getNPV(mu); double rho = getRho(mu); @@ -143,11 +143,11 @@ void setEtaPtE(FactorizedJetCorrector *jec, double eta, double pt, double e, if(cone == "7") r=0.7 ; if(cone == "8") r=0.8 ; if(cone == "9") r=0.9 ; - double jeta = TMath::Pi()*r*r; + double jeta = TMath::Pi()*r*r; jec->setJetA(jeta); jec->setRho(rho); - + return; } @@ -163,7 +163,7 @@ Double_t funcCorrPt(Double_t *x, Double_t *p) { double e = pt * cosh(eta); double mu = p[1]; setEtaPtE(_thejec, eta, pt, e, mu); - + return (_thejec->getCorrection() * pt); } @@ -175,15 +175,15 @@ double getEtaPtE(FactorizedJetCorrector *jec, double eta, double pt, double e, setEtaPtE(jec, eta, pt, e, mu); // if using pTgen, need to iterate to solve ptreco - if (_useptgen) { + if (_useptgen) { double ptgen = pt; - _thejec = jec; + _thejec = jec; fCorrPt->SetParameters(eta, mu); // Find ptreco that gives pTreco*JEC = pTgen double ptreco = fCorrPt->GetX(ptgen,5,6500); - + setEtaPtE(jec, eta, ptreco, e, mu); } @@ -203,7 +203,7 @@ double getEtaPtUncert(JetCorrectionUncertainty *unc, unc->setJetEta(eta); unc->setJetPt(pt); - // if not using pTgen, need to solve it + // if not using pTgen, need to solve it if (!_useptgen) { assert(jec); @@ -211,7 +211,7 @@ double getEtaPtUncert(JetCorrectionUncertainty *unc, _thejec = jec; fCorrPt->SetParameters(eta, mu); double ptgen = fCorrPt->Eval(ptreco); - + unc->setJetEta(eta); unc->setJetPt(ptgen); } @@ -222,7 +222,7 @@ double getEtaPtUncert(JetCorrectionUncertainty *unc, // ------------------------------------------------------------ // -void compareJEC(string payld1="Winter14_V8", string payld2="", string payld3="", +void compareJEC(string payld1="Winter14_V8", string payld2="", string payld3="", string algo1="AK5PFchs", string algo2="AK5PFchs", string algo3 ="AK5PFchs", string type1="DATA", string type2="DATA", string type3 ="DATA", bool l1=true, bool l2l3=true, bool res=true) { @@ -254,17 +254,21 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa int maxTries = 7; string strPath; - vector paths = {"CondFormats/JetMETObjects/data/"}; - paths.push_back(string("/fdata/hepx/store/user/aperloff/JEC/80X_Summer16/")+cid1+"/"); - paths.push_back(string("/fdata/hepx/store/user/aperloff/JEC/80X_Summer16/")+cid2+"/"); - paths.push_back(string("/fdata/hepx/store/user/aperloff/JEC/80X_Summer16/")+cid3+"/"); - paths.push_back(string("/home/aperloff/JEC/CMSSW_8_0_20/src/JetMETCorrections/JECDatabase/textFiles/")+cid1+"/"); - paths.push_back(string("/home/aperloff/JEC/CMSSW_8_0_20/src/JetMETCorrections/JECDatabase/textFiles/")+cid2+"/"); - paths.push_back(string("/home/aperloff/JEC/CMSSW_8_0_20/src/JetMETCorrections/JECDatabase/textFiles/")+cid3+"/"); + vector paths = {"","CondFormats/JetMETObjects/data/"}; + paths.push_back(string("/fdata/hepx/store/user/siluo/JEC/94X_Fall17/")+cid1+"/"); + paths.push_back(string("/fdata/hepx/store/user/siluo/JEC/80X_Summer16/")+cid2+"/"); + paths.push_back(string("/fdata/hepx/store/user/siluo/JEC/80X_stitch/")+cid3+"/"); + // paths.push_back(string("/fdata/hepx/store/user/aperloff/JEC/80X_Summer16/")+cid1+"/"); + // paths.push_back(string("/fdata/hepx/store/user/aperloff/JEC/80X_Summer16/")+cid2+"/"); + // paths.push_back(string("/fdata/hepx/store/user/aperloff/JEC/80X_Summer16/")+cid3+"/"); + paths.push_back(string("/home/sifuluo/jec/CMSSW_8_0_20/src/JetMETCorrections/JECDatabase/textFiles/")+cid1+"/"); + paths.push_back(string("/home/sifuluo/jec/CMSSW_8_0_20/src/JetMETCorrections/JECDatabase/textFiles/")+cid2+"/"); + paths.push_back(string("/home/sifuluo/jec/CMSSW_8_0_20/src/JetMETCorrections/JECDatabase/textFiles/")+cid3+"/"); + // JEC1 bool noL2L3_p1 = TString(payld1).Contains("Run1Sub") || TString(payld1).Contains("nsRunIISpring15") || TString(payld1).Contains("RC") ; - // + // // noL2L3_p1 = true ; //ctmp JetCorrectionUncertainty *jecUnc1 = 0; @@ -311,7 +315,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa JetCorrectorParameters *JetCorPar2 =0; JetCorrectionUncertainty *jecUnc2 =0; if(payld2 != ""){ - bool noL2L3_p2 = TString(payld2).Contains("Run1Sub") || TString(payld2).Contains("nsRunIISpring15") || TString(payld2).Contains("RC") ; + bool noL2L3_p2 = TString(payld2).Contains("Run1Sub") || TString(payld2).Contains("nsRunIISpring15") || TString(payld2).Contains("RC") ; strVec.erase(strVec.begin(),strVec.end()); strVec.push_back(Form("%s_L1FastJet_%s.txt",cid2,a2)); @@ -375,11 +379,11 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa } if((type2=="MC" && istr!=3) || type2=="DATA")cout << strPath << endl << flush; if(istr==0) - JetCorPar3L1 = new JetCorrectorParameters(strPath); + JetCorPar3L1 = new JetCorrectorParameters(strPath); else if(istr==1) - JetCorPar3L2 = TString(payld3).Contains("Run1Sub") ? 0 : new JetCorrectorParameters(strPath); + JetCorPar3L2 = TString(payld3).Contains("Run1Sub") ? 0 : new JetCorrectorParameters(strPath); else if(istr==2) - JetCorPar3L3 = TString(payld3).Contains("Run1Sub") ? 0 : new JetCorrectorParameters(strPath); + JetCorPar3L3 = TString(payld3).Contains("Run1Sub") ? 0 : new JetCorrectorParameters(strPath); else if(istr==3) JetCorPar3 = (type3 == "MC" || TString(payld3).Contains("Run1Sub") ? 0 : new JetCorrectorParameters(strPath)); else if(istr==4) @@ -397,16 +401,16 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa // JEC2 vector vParam2; - FactorizedJetCorrector *JEC2 = 0 ; + FactorizedJetCorrector *JEC2 = 0 ; if(payld2 != ""){ if (l1) vParam2.push_back(*JetCorPar2L1); if (l2l3) vParam2.push_back(*JetCorPar2L2); if (l2l3) vParam2.push_back(*JetCorPar2L3); if (res && type2 == "DATA") vParam2.push_back(*JetCorPar2); JEC2 = new FactorizedJetCorrector(vParam2); - } + } // JEC3 - vector vParam3; + vector vParam3; FactorizedJetCorrector *JEC3 = 0 ; if(payld3 != ""){ if (l1) vParam3.push_back(*JetCorPar3L1); @@ -439,7 +443,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa // 507, 592, 686, 790, 905, 1032, 1172, 1327, 1497, 1684, 1890, //1999}; // 2000, 2238, 2500};//, 2787, 3103, 3450}; - {10, 12, 15, 18, 21, 24, + {10, 12, 15, 18, 21, 24, 28, 32, 37, 43, 49, 56, 64, 74, 84, 97, 114, 133, 153, 174, 196, 220, 245, 272, 300, 362, 430, 507, 592, 686, 790, 905, 1032, 1172, 1327, 1497, 1684, 1890, //1999}; @@ -460,7 +464,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa hpt->SetMinimum(0.); hpt->SetMaximum(2.); - if (_paper) { + if (_paper) { if (l1 && !l2l3 && !res) h->SetYTitle("Pileup offset correction"); if (!l1 && l2l3 && !res) h->SetYTitle("Simulated response correction"); if (!l1 && !l2l3 && res) h->SetYTitle("Residual response correction"); @@ -469,7 +473,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa if (!l1 && l2l3 && !res) hpt->SetYTitle("Simulated response correction"); if (!l1 && !l2l3 && res) hpt->SetYTitle("Residual response correction"); // - if (l1 && !l2l3 && !res) h->GetYaxis()->SetRangeUser(0.5,1.4); + if (l1 && !l2l3 && !res) h->GetYaxis()->SetRangeUser(0.5,1.4); if (!l1 && l2l3 && !res) h->GetYaxis()->SetRangeUser(0.85,1.8); // 2.2 if (!l1 && !l2l3 && res) h->GetYaxis()->SetRangeUser(0.85,1.4); @@ -483,7 +487,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa // if (l1 && !l2l3 && !res) hpt->GetYaxis()->SetRangeUser(0.5,1.4); if (!l1 && l2l3 && !res) hpt->GetYaxis()->SetRangeUser(0.85,2.2); // 1.4 - + if (!l1 && !l2l3 && res) hpt->GetYaxis()->SetRangeUser(0.85,1.4); if(_pt<20.){ @@ -524,7 +528,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa TGraph *g2a_mn = new TGraph(0); TGraph *g2d_pl = new TGraph(0); TGraph *g2d_mn = new TGraph(0); - + TGraphErrors *g3a_e = new TGraphErrors(0); TGraphErrors *g3d_e = new TGraphErrors(0); TGraph *g3a_pl = new TGraph(0); @@ -535,7 +539,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa double ptbins[] = {15., 30, 100, 1000, 4000}; // double ptbins[] = {1000, 2000, 3000, 3500, 4000}; - const int npt = sizeof(ptbins)/sizeof(ptbins[0]); + const int npt = sizeof(ptbins)/sizeof(ptbins[0]); TGraphErrors *g21s[npt]; TMultiGraph *mg = new TMultiGraph(); @@ -554,7 +558,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa for (int j = 0; j != npt; ++j) { TGraphErrors *g21 = g21s[j]; double pt = ptbins[j]; - double energy = pt*cosh(eta); + double energy = pt*cosh(eta); if (energy < 6500.) { // if (energy < 3000.) { @@ -565,20 +569,20 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa + getEtaPtE(JEC1, -eta, pt, energy)); _alg = algo2; _payld=payld2; - if(i==1 && j==0) cout << _payld << ": mu and rho " << _mu << " " << getRho(_mu) << endl ; + if(i==1 && j==0) cout << _payld << ": mu and rho " << _mu << " " << getRho(_mu) << endl ; double y2 = 0.5*(getEtaPtE(JEC2, +eta, pt, energy) + getEtaPtE(JEC2, -eta, pt, energy)); g21->SetPoint(g21->GetN(), eta, y2/y1); } // energy < 6500 } // for j - } + } // ***** Fixed pT, versus eta ****** - { + { double pt = _pt ; - double energy = pt*cosh(eta); - + double energy = pt*cosh(eta); + if (energy < 6500.) { // if (energy < 3000.) { @@ -625,7 +629,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa y2 = JEC2 ? getEtaPtE(JEC2, +eta, pt, energy) : 0; _alg = algo3; _payld=payld3; - y3 = JEC3 ? getEtaPtE(JEC3, +eta, pt, energy) : 0; + y3 = JEC3 ? getEtaPtE(JEC3, +eta, pt, energy) : 0; } double e1 = jecUnc1 ? getEtaPtUncert(jecUnc1, JEC1, eta, pt) : 0; @@ -655,26 +659,26 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa // ...... - for (int ieta = -49; ieta<=49; ieta++){ // Loop over eta bins: + for (int ieta = -49; ieta<=49; ieta++){ // Loop over eta bins: // Only one eta value in the main plot // Many eta curves for the _validation_ plot - double eta = ieta*0.1 ; + double eta = ieta*0.1 ; TGraph *gtmp = new TGraph(0); for (int i = 1; i != hpt->GetNbinsX()+1; ++i) { // Loop over pt-bins - + double pt = hpt->GetBinCenter(i); double energy = pt*cosh(eta); - + if (pt>0. && energy < 6500.) { // if (pt>0. && energy < 3000.) { - + // Asymmetric corrections now _alg = algo1; _payld=payld1; double y1 = 0.5*(getEtaPtE(JEC1, +eta, pt, energy) - + getEtaPtE(JEC1, -eta, pt, energy)); + + getEtaPtE(JEC1, -eta, pt, energy)); double y2(0); _alg = algo2; _payld=payld2; @@ -715,7 +719,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa double e1 = jecUnc1 ? getEtaPtUncert(jecUnc1, JEC1, eta, pt) : 0; double e2 = jecUnc2 ? getEtaPtUncert(jecUnc2, JEC2, eta, pt) : 0; // double e3 = jecUnc3 ? getEtaPtUncert(jecUnc3, JEC3, eta, pt) : 0; - if(ieta == _ieta){ + if(ieta == _ieta){ // cout << " pt y1 y2 " << pt << " " << y1 << " " << y2 << endl ; g1d->SetPoint(g1d->GetN(), pt, y1); g2d->SetPoint(g2d->GetN(), pt, y2); @@ -734,19 +738,19 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa gtmp->SetPoint(gtmp->GetN(), pt, y1); if(l1 && !l2l3 && ! res) { if(y1>1.) { - cout << endl ; + cout << endl ; cout << " L1: pt eta cor " << std::setw(7) << pt << " " << std::setw(7) << eta << " " << y1 ; if(y1 > 1.1) - cout << " < ============ " ; + cout << " < ============ " ; cout << endl ; } } } } // END: loop over pt bins gtmp->SetLineWidth(2); -// gtmp->SetLineStyle( ieta>=0 ? 2 : 1); +// gtmp->SetLineStyle( ieta>=0 ? 2 : 1); gtmp->SetLineColor(abs(ieta)/10+1); - if(ieta%10 ==0 && ieta>=0 && ieta<50) + if(ieta%10 ==0 && ieta>=0 && ieta<50) leg10->AddEntry(gtmp,Form("%2.0f<|#eta|<%2.0f",0.1*abs(ieta),0.1*abs(ieta+10)),"l"); mg->Add(gtmp,"l"); } // END: loop over eta bins @@ -790,20 +794,20 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa g2a_pl->Draw("SAMEL"); g2a_mn->SetLineColor(kRed-9); g2a_mn->SetLineStyle(kSolid);//kDotted); - g2a_mn->Draw("SAMEL"); + g2a_mn->Draw("SAMEL"); } if(JEC3){ - g3a_e->SetFillStyle(3003); + g3a_e->SetFillStyle(3003); g3a_e->SetFillColor(kGreen+2); g3a_e->Draw("SAME E3"); g3a_pl->SetLineColor(kGreen-9); g3a_pl->SetLineStyle(kSolid);//kDotted); - g3a_pl->Draw("SAMEL"); + g3a_pl->Draw("SAMEL"); g3a_mn->SetLineColor(kGreen-9); g3a_mn->SetLineStyle(kSolid);//kDotted); g3a_mn->Draw("SAMEL"); - } + } g1a->SetMarkerStyle(kFullSquare); g1a->SetMarkerColor(kBlue); @@ -819,14 +823,14 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa if (JEC3){ g3a->SetMarkerStyle(kOpenSquare); - g3a->SetMarkerColor(kGreen+2); + g3a->SetMarkerColor(kGreen+2); g3a->SetLineColor(kGreen+2); g3a->Draw("SAMEPL"); } TLatex *tex = new TLatex(); tex->SetNDC(); - tex->SetTextSize(0.045); + tex->SetTextSize(0.045); tex->DrawLatex(0.19,0.75,Form("p_{T,%s} = %2.0f GeV",cgen,_pt)); if (l1) tex->DrawLatex(0.19,0.68,Form("#LT#mu#GT = %1.1f",_mu)); @@ -846,8 +850,8 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa if(chs2=="chs") chs2 = "+CHS" ; if(chs2=="Pup") chs2 = "+Puppi" ; if(pf2=="Ca") {pf2=""; chs2 = "Calo" ;} - - + + string cone3 = algo3.substr(2, 1); string pf3 = algo3.substr(3, 2); string chs3 = algo3.substr(5, 3); @@ -879,15 +883,15 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa slg3 = payld3 ; } else { sheader =payld1 ; - sheader_pdf ="_"+payld1; + sheader_pdf ="_"+payld1; } if( (JEC2 && cone2 != cone1) || (JEC3 && cone3 !=cone1) ){ if(slg1 != "") slg1 +=", "; if(slg2 != "") slg2 +=", "; if(slg3 != "") slg3 +=", "; slg1 += "R = 0." + cone1 ; - slg2 += "R = 0." + cone2 ; - slg3 += "R = 0." + cone3 ; + slg2 += "R = 0." + cone2 ; + slg3 += "R = 0." + cone3 ; } else { if(sheader != "") sheader +=", "; sheader += "R = 0." + cone1 ; @@ -897,7 +901,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa if(slg1 != "") slg1 +=", "; if(slg2 != "") slg2 +=", "; if(slg3 != "") slg3 +=", "; - slg1 += pf1 + chs1 ; + slg1 += pf1 + chs1 ; slg2 += pf2 + chs2 ; slg3 += pf3 + chs3 ; } else { @@ -917,7 +921,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa sheader += type1 ; sheader_pdf += "_" + type1 ; } - cout << "string size " << sheader.size() << " " << slg1.size() << " " << slg2.size() << endl ; + cout << "string size " << sheader.size() << " " << slg1.size() << " " << slg2.size() << endl ; if(sheader.size()>25 || slg1.size()>20 || slg2.size()>20 || slg3.size()>20 ) leg1a->SetTextSize(0.031); leg1a->SetTextSize(0.029); @@ -929,7 +933,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa ///// leg1a->AddEntry(g2a,"2016","LPF"); } - + gPad->RedrawAxis(); @@ -944,12 +948,12 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa g1d_pl->SetLineColor(kBlue-9); g1d_pl->SetLineStyle(kSolid); g1d_pl->Draw("SAMEL"); - g1d_mn->SetLineColor(kBlue-9); + g1d_mn->SetLineColor(kBlue-9); g1d_mn->SetLineStyle(kSolid); g1d_mn->Draw("SAMEL"); - + if (JEC2) { - g2d_e->SetFillStyle(3003); + g2d_e->SetFillStyle(3003); g2d_e->SetFillColor(kRed); g2d_e->Draw("SAME E3"); g2d_pl->SetLineColor(kRed-9); @@ -957,15 +961,15 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa g2d_pl->Draw("SAMEL"); g2d_mn->SetLineColor(kRed-9); g2d_mn->SetLineStyle(kSolid); - g2d_mn->Draw("SAMEL"); + g2d_mn->Draw("SAMEL"); } if (JEC3) { g3d_e->SetFillStyle(3003); g3d_e->SetFillColor(kGreen+2); g3d_e->Draw("SAME E3"); - g3d_pl->SetLineColor(kGreen-9); - g3d_pl->SetLineStyle(kSolid); + g3d_pl->SetLineColor(kGreen-9); + g3d_pl->SetLineStyle(kSolid); g3d_pl->Draw("SAMEL"); g3d_mn->SetLineColor(kGreen-9); g3d_mn->SetLineStyle(kSolid); @@ -976,20 +980,20 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa g1d->SetMarkerColor(kBlue); g1d->SetLineColor(kBlue); g1d->Draw("SAMEPL"); - + if (JEC2) { g2d->SetMarkerStyle(kFullCircle); - g2d->SetMarkerColor(kRed); + g2d->SetMarkerColor(kRed); g2d->SetLineColor(kRed); g2d->Draw("SAMEPL"); } if (JEC3) { g3d->SetMarkerStyle(kOpenSquare); - g3d->SetMarkerColor(kGreen+2); + g3d->SetMarkerColor(kGreen+2); g3d->SetLineColor(kGreen+2); g3d->Draw("SAMEPL"); - } + } // tex->DrawLatex(0.19,0.75,"|#eta| = 0"); @@ -1019,11 +1023,11 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa h1ar->SetMaximum(1.55); h1ar->GetYaxis()->SetTitle(Form("%s%s%s%s (%s / %s)",cl1,cl2l3,cpl,cres, slg2.c_str(), slg1.c_str())); -//// h1ar->GetYaxis()->SetTitle(Form("%s%s%s%s (80X / 76X)",cl1,cl2l3,cpl,cres )); // this is tmp +//// h1ar->GetYaxis()->SetTitle(Form("%s%s%s%s (80X / 76X)",cl1,cl2l3,cpl,cres )); // this is tmp h1ar->GetYaxis()->SetTitleSize(0.045); - h1ar->GetYaxis()->SetTitleOffset(1.7); - + h1ar->GetYaxis()->SetTitleOffset(1.7); + TLine *l = new TLine(); l->SetLineStyle(kDashed); @@ -1033,9 +1037,9 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa l->DrawLine(0,1.02,5,1.02); l->DrawLine(0,0.98,5,0.98); - if (l1) tex->DrawLatex(0.19,0.68,Form("#LT#mu#GT = %1.1f",_mu)); + if (l1) tex->DrawLatex(0.19,0.68,Form("#LT#mu#GT = %1.1f",_mu)); + - TLegend *leg = tdrLeg(0.46, 0.57, 0.90, 0.88); if(sheader.size()>25) leg->SetTextSize(0.031); @@ -1045,7 +1049,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa gr->SetLineStyle(styles[i]);//i+1); gr->SetLineWidth(3); gr->Draw("SAME L"); - + leg->SetHeader(sheader.c_str()); leg->AddEntry(gr,Form("p_{T,%s} = %1.0f GeV", cgen,ptbins[i]),"L"); @@ -1054,7 +1058,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa gPad->RedrawAxis(); c1ar->SaveAs(Form("pdf/compareJECversions_%s%s_2over1_%s.pdf",cs,sheader_pdf.c_str(),cgen)); - } + } { c10->cd(); @@ -1078,7 +1082,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa tex->SetTextSize(0.040); tex->DrawLatex(0.19,0.89, - Form("%s %s R = 0.%s, %s%s", payld1.c_str(), type1.c_str(), cone1.c_str(), pf1.c_str(), chs1.c_str())); + Form("%s %s R = 0.%s, %s%s", payld1.c_str(), type1.c_str(), cone1.c_str(), pf1.c_str(), chs1.c_str())); if (l1) tex->DrawLatex(0.19,0.68,Form("#LT#mu#GT = %1.1f",_mu)); @@ -1101,16 +1105,15 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa //______________________________________________________________________________ int main(int argc,char**argv) { - compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK4PFchs", "AK4PFchs", "AK4PFchs", "MC","MC","MC", false, true, false); - compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK4PFchs", "AK4PFchs", "AK4PFchs", "MC","MC","MC", true, false, false); - - compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFchs", "AK8PFchs", "AK8PFchs", "MC","MC","MC", false, true, false); - compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFchs", "AK8PFchs", "AK8PFchs", "MC","MC","MC", true, false, false); + compareJEC("Fall17_25nsV1", "Summer16_25nsV5", "bias2SelectionPow_25nsV1", "AK4PFchs", "AK4PFchs", "AK4PFchs", "MC","MC","MC", false, true, false); + compareJEC("Fall17_25nsV1", "Summer16_25nsV5", "bias2SelectionPow_25nsV1", "AK4PFchs", "AK4PFchs", "AK4PFchs", "MC","MC","MC", true, false, false); - compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK4PFPuppi", "AK4PFPuppi", "AK4PFPuppi", "MC","MC","MC", false, true, false); - compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK4PFPuppi", "AK4PFPuppi", "AK4PFPuppi", "MC","MC","MC", true, false, false); - - compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFPuppi", "AK8PFPuppi", "AK8PFPuppi", "MC","MC","MC", false, true, false); - compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFPuppi", "AK8PFPuppi", "AK8PFPuppi", "MC","MC","MC", true, false, false); + // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFchs", "AK8PFchs", "AK8PFchs", "MC","MC","MC", false, true, false); + // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFchs", "AK8PFchs", "AK8PFchs", "MC","MC","MC", true, false, false); + // + // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK4PFPuppi", "AK4PFPuppi", "AK4PFPuppi", "MC","MC","MC", false, true, false); + // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK4PFPuppi", "AK4PFPuppi", "AK4PFPuppi", "MC","MC","MC", true, false, false); + // + // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFPuppi", "AK8PFPuppi", "AK8PFPuppi", "MC","MC","MC", false, true, false); + // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFPuppi", "AK8PFPuppi", "AK8PFPuppi", "MC","MC","MC", true, false, false); } - diff --git a/JetAnalyzers/bin/jet_synchtest_x.cc b/JetAnalyzers/bin/jet_synchtest_x.cc index 7e16d657..0dc7bf18 100644 --- a/JetAnalyzers/bin/jet_synchtest_x.cc +++ b/JetAnalyzers/bin/jet_synchtest_x.cc @@ -124,7 +124,7 @@ class MatchEventsAndJets { int NBinsNpvRhoNpu; int npvRhoNpuBinWidth; int iIT; - int inpv; + int inpv; int inpv_low; int inpv_high; int irho; @@ -339,9 +339,10 @@ void MatchEventsAndJets::WriteMatchedEventsMaps(map, evtid> outputFilename = "matchedEventsMaps_"+algo1+".root"; outputFilename = outputPath+outputFilename; string name = "mapTree"; - cout << "\tWriting " << name+(noPU ? "NoPU" : "PU") << " to " << outputFilename << " ... " << flush; + cout << "\tWriting " << name+(noPU ? "NoPU" : "PU") << " to " << outputFilename << " ... " << flush; string option = (noPU ? "UPDATE" : "RECREATE"); TFile* mapFile = TFile::Open(outputFilename.c_str(),option.c_str()); + mapFile->cd(); mapFile->WriteObject(&mapTree,(name+(noPU ? "NoPU" : "PU")).c_str()); mapFile->Write(); mapFile->Close(); @@ -440,7 +441,7 @@ void MatchEventsAndJets::DeclareHistograms(bool reduceHistograms) { histograms["g_GenWeight"] = new TH1D("g_GenWeight", "g_GenWeight;log_{10}(GenWeight);Events", 1000,-48,2); histograms["g_pThatWeight"] = new TH1D("g_pThatWeight;log_{10}(pThatWeight);Events","g_pThatWeight", 1000,-48,2); histograms["g_weight"] = new TH1D("g_weight","g_weight;log_{10}(EvtWeight);Events", 1000,-48,2); - histograms["g_pthat"] = new TH1D("g_pthat","g_pthat;#hat{p}_{T}^{PU};Events",(int)vpt[NPtBins]/10.0,vpt[0],vpt[NPtBins]); + histograms["g_pthat"] = new TH1D("g_pthat","g_pthat;#hat{p}_{T}^{PU};Events",(int)vpt[NPtBins]/10.0,vpt[0],vpt[NPtBins]); if(!reduceHistograms) { histograms["g_nj"] = new TH2D("g_nj","g_nj",30,0,30,30,0,30); histograms["g_npv"] = new TH2D("g_npv","g_npv",50,0,50,50,0,50); @@ -795,6 +796,7 @@ void MatchEventsAndJets::LoopOverEvents(bool verbose, bool reduceHistograms, str for (IT::const_iterator it = mapTreePU.begin(); it != mapTreePU.end(); ++it) { if (iftest && nevs >= maxEvts) return; + // if (nevs >= 10000000) return; //if (nevs%10000==0) cout << "\t"<= 0 && j2 >= 0 && j1 < tpu->nref && j2 < tnopu->nref && - tpu->refdrjt->at(j1) < maxDeltaR && tnopu->refdrjt->at(j2) < maxDeltaR && + tpu->refdrjt->at(j1) < maxDeltaR && tnopu->refdrjt->at(j2) < maxDeltaR && fabs(tpu->refpt->at(j1) - tnopu->refpt->at(j2))<0.0001) { jetMap[j1] = j2; } @@ -1186,8 +1188,8 @@ bool MatchEventsAndJets::FillHistograms(bool reduceHistograms) { double avg_offset = 0; double avg_offset_det[NDetectorNames] = {0,0,0,0}; double njet_det[NDetectorNames] = {0,0,0,0}; - - // MATCHING HISTOS. + + // MATCHING HISTOS. // Loop over matched jets int jpu = -1; int jnopu = -1; @@ -1356,7 +1358,7 @@ bool MatchEventsAndJets::FillHistograms(bool reduceHistograms) { dynamic_cast(histograms[hname])->Fill(tpu->refpt->at(jpu),respNopu,weight); hname = Form("p_offAfterOoffBeforeVsrefpt_%s_npv%i_%i",detectorAbbreviation.Data(),inpv_low,inpv_high); dynamic_cast(histograms[hname])->Fill(tpu->refpt->at(jpu),offset/offset_raw,weight); - + //RHO hname = Form("p_resVsrefpt_%s_rho%i_%i",detectorAbbreviation.Data(),irho_low,irho_high); dynamic_cast(histograms[hname])->Fill(tpu->refpt->at(jpu),resp,weight); @@ -1372,7 +1374,7 @@ bool MatchEventsAndJets::FillHistograms(bool reduceHistograms) { dynamic_cast(histograms[hname])->Fill(tpu->refpt->at(jpu),respNopu,weight); hname = Form("p_offAfterOoffBeforeVsrefpt_%s_rho%i_%i",detectorAbbreviation.Data(),irho_low,irho_high); dynamic_cast(histograms[hname])->Fill(tpu->refpt->at(jpu),offset/offset_raw,weight); - + //OTHER hname = Form("p_resVsnpu_%s_pt%.1f_%.1f",detectorAbbreviation.Data(), vpt[JetInfo::getBinIndex(tpu->refpt->at(jpu),vpt,NPtBins)],vpt[JetInfo::getBinIndex(tpu->refpt->at(jpu),vpt,NPtBins)+1]); @@ -1389,7 +1391,7 @@ bool MatchEventsAndJets::FillHistograms(bool reduceHistograms) { avg_offset_det[idet]+=offset; njet_det[idet]+=1.; - } // for matched jets + } // for matched jets avg_offset /= jetMap.size(); for (int det=0;det/.root', + # 'root://cmsxrootd.fnal.gov//store/mc//.root' + 'root://cmsxrootd.fnal.gov//store/mc/RunIIFall17DRPremix/QCD_Pt-15to7000_TuneCP5_Flat_13TeV_pythia8/AODSIM/94X_mc2017_realistic_v10-v1/50000/00304636-1BDB-E711-B6F3-FA163ECE02A9.root', + 'root://cmsxrootd.fnal.gov//store/mc/RunIIFall17DRPremix/QCD_Pt-15to7000_TuneCP5_Flat_13TeV_pythia8/AODSIM/94X_mc2017_realistic_v10-v1/50000/0036C92E-DFDB-E711-952A-008CFAFC05DE.root', ) process.source = cms.Source("PoolSource", fileNames = inputFiles ) diff --git a/JetUtilities/bin/BuildFile.xml b/JetUtilities/bin/BuildFile.xml index bea2eaff..585839b7 100644 --- a/JetUtilities/bin/BuildFile.xml +++ b/JetUtilities/bin/BuildFile.xml @@ -12,5 +12,3 @@ - - diff --git a/JetUtilities/interface/ClosureMaker.hh b/JetUtilities/interface/ClosureMaker.hh index 96dbcaa0..b202e41f 100644 --- a/JetUtilities/interface/ClosureMaker.hh +++ b/JetUtilities/interface/ClosureMaker.hh @@ -83,7 +83,7 @@ public: void checkResponse(); pair determineCanvasRange(double xmin, double xmax); void makeCanvases(); - void makeMergedCanvas(); + void makeMergedCanvas(bool finemerge); void writeToFile(); void makeClosure(const VARIABLES::Variable ivar = VARIABLES::refpt); diff --git a/JetUtilities/src/ClosureMaker.cc b/JetUtilities/src/ClosureMaker.cc index 77dc7878..c6e1cb24 100644 --- a/JetUtilities/src/ClosureMaker.cc +++ b/JetUtilities/src/ClosureMaker.cc @@ -88,7 +88,7 @@ void ClosureMaker::openInputFile() { ifile = TFile::Open(path+filename+"_"+flavor+".root","READ"); else ifile = TFile::Open(path+filename+".root","READ"); - + if(ifile == 0) { cout << "ERROR::ClosureMaker::openInputFiles() Could not open the file " << path << endl; std::terminate(); @@ -102,7 +102,7 @@ void ClosureMaker::getHistograms(TDirectoryFile* idir) { } if(var == VARIABLES::refpt || var == VARIABLES::ptclpt) { - hl.load_objects(idir,"RelRspVsRefPt:JetEta"); + hl.load_objects(idir,"RelRspVsRefPt:JetEta"); objects_loaded = true; if(hl.nobjects()!=NDetectorNames) { cout << "One or more of the histogram pointers from file " << path << " is NULL." << endl; @@ -125,7 +125,7 @@ void ClosureMaker::getHistograms(TDirectoryFile* idir) { void ClosureMaker::openOutputFile() { // // Open/create the output directory and file - // + // TString ofname = Form("%s/ClosureVs%s.root",outputDir.Data(),getVariableTitleString(var).c_str()); if(!flavor.IsNull()) ofname = Form("%s/ClosureVs%s_%s.root",outputDir.Data(), getVariableTitleString(var).c_str(),flavor.Data()); @@ -209,7 +209,7 @@ void ClosureMaker::loopOverAlgorithms() { if (strcmp(dirKey->GetClassName(),"TDirectoryFile")!=0) continue; TDirectoryFile* idir = (TDirectoryFile*)dirKey->ReadObj(); alg = idir->GetName(); if (!JetInfo::contains(algs,alg)) continue; - + algIndex++; cout << alg << " ... " << endl; @@ -258,7 +258,7 @@ void ClosureMaker::loopOverAlgorithms() { // makeCanvases(); if(var == VARIABLES::refpt || var == VARIABLES::jtpt || var == VARIABLES::ptclpt) { - makeMergedCanvas(); + makeMergedCanvas(false); } writeToFile(); @@ -421,7 +421,7 @@ void ClosureMaker::checkResponse() { cout << "\tWARNING Closure for " << hClosure[ih]->GetName() << " at " << hClosure[ih]->GetBinLowEdge(ibin) << " < " << getVariableTitleString(var) << " < " << hClosure[ih]->GetBinLowEdge(ibin+1) - <<" has relresp as low as " << binCont << " +/- " << binErr << endl; + <<" has relresp as low as " << binCont << " +/- " << binErr << endl; } } } @@ -491,7 +491,7 @@ void ClosureMaker::makeCanvases() { tdrLeg(0.58,0.16,0.9,0.4))); if((var == VARIABLES::refpt || var == VARIABLES::jtpt || var == VARIABLES::ptclpt) && ih<3) canvases_legends.back().first->GetPad(0)->SetLogx(); - + // // Format and draw the pave // @@ -533,10 +533,10 @@ void ClosureMaker::makeCanvases() { } //______________________________________________________________________________ -void ClosureMaker::makeMergedCanvas() { +void ClosureMaker::makeMergedCanvas(bool finemerge = false) { TString name = Form("ClosureVs%s_Overview_%s",getVariableTitleString(var).c_str(),alg.Data()); if(!flavor.IsNull()) name+="_"+flavor; - + if (finemerge) name+="_Fine"; // // Setup the frame, canvas, legend, and pave // @@ -548,10 +548,11 @@ void ClosureMaker::makeMergedCanvas() { frame->GetXaxis()->SetLimits(XminCalo[0],Xmax[0]); frame->GetXaxis()->SetMoreLogLabels(); frame->GetXaxis()->SetNoExponent(); - pair range = determineCanvasRange(frame->GetXaxis()->GetXmin(),frame->GetXaxis()->GetXmax()); - //frame->GetYaxis()->SetRangeUser(0.95,1.05); + pair range = determineCanvasRange(frame->GetXaxis()->GetXmin(),frame->GetXaxis()->GetXmax()); + if (range.second - range.first > 0.2 && !finemerge) makeMergedCanvas(true); + if (finemerge) frame->GetYaxis()->SetRangeUser(0.95,1.05); //frame->GetYaxis()->SetRangeUser(0.35,1.35); - frame->GetYaxis()->SetRangeUser(range.first,range.second); + else frame->GetYaxis()->SetRangeUser(range.first,range.second); frame->GetXaxis()->SetTitle(getVariableAxisTitleString(var).c_str()); frame->GetYaxis()->SetTitle("Response"); canvases_legends.push_back(make_pair(tdrCanvas(name,frame,14,11,true), @@ -635,5 +636,3 @@ void ClosureMaker::makeClosure(const VARIABLES::Variable ivar) { loopOverAlgorithms(); closeFiles(); } - - diff --git a/JetUtilities/src/L2Creator.cc b/JetUtilities/src/L2Creator.cc index de67892a..fd6d234a 100644 --- a/JetUtilities/src/L2Creator.cc +++ b/JetUtilities/src/L2Creator.cc @@ -232,7 +232,7 @@ void L2Creator::loopOverEtaBins() { // only add points to the graphs if the current histo is not empty // the current setting might be a little high // - if (hrsp->GetEntries() > 4) {//hrsp->Integral()!=0) { + if (hrsp->GetEntries() > 400) {//hrsp->Integral()!=0) { //TF1* frsp = (TF1*)hrsp->GetListOfFunctions()->Last(); //std::cout << "hrspName = " << hrsp->GetName() << ": frsp = " << frsp << std::endl; @@ -1313,6 +1313,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //For eta-dependent spline clipping int pt_limit = 70; + int pt_clip = 8; unsigned int vector_size = 0; vector_size = vabscor_eta.size(); @@ -1345,6 +1346,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { bool abovePtLimit = false; bool lastLine = false; + bool firstline = true; for(int isection=0; isectiongetNSections(); isection++) { if(lastLine) continue; @@ -1364,7 +1366,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //if(isection==spline->getNSections()-1) { // bounds.second = 6500; //} - + if(bounds.second < pt_clip) continue; if(isection==spline->getNSections()-1) lastLine = true; if(bounds.second >= pt_limit) { abovePtLimit = true; @@ -1373,16 +1375,17 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //For expediency of Summer16_25nsV5_MC do eta-dependent clipping fout<getNpar()+2) //Number of parameters + 2 - <setParameters(isection); for(int p=0; pgetNpar(); p++) { fout<GetParameter(p); } fout< Date: Mon, 5 Feb 2018 17:15:59 -0600 Subject: [PATCH 02/10] Adding functions of clipping pt in lower pt and statistical threshold for L2 --- JetAnalyzers/test/run_JRA_cfg.py | 6 +- JetUtilities/interface/ClosureMaker.hh | 7 +- JetUtilities/interface/L2Creator.hh | 3 +- JetUtilities/src/ClosureMaker.cc | 3 +- JetUtilities/src/L2Creator.cc | 7 +- JetUtilities/src/SynchFittingProcedure.hh | 94 +++++++++++------------ 6 files changed, 64 insertions(+), 56 deletions(-) diff --git a/JetAnalyzers/test/run_JRA_cfg.py b/JetAnalyzers/test/run_JRA_cfg.py index 222679e5..340a2db4 100644 --- a/JetAnalyzers/test/run_JRA_cfg.py +++ b/JetAnalyzers/test/run_JRA_cfg.py @@ -77,8 +77,10 @@ print "Couldn't open the external list of files from DAS. If you just checkout out the JetResponseAnalyzer package you will need to make this file yourself. Currently Falling back to opening the list hard-coded in run_JRA_cfg.py. This is not a bad action as long as it is what you intended to have happen." inputFiles = cms.untracked.vstring( # 'root://cmsxrootd.fnal.gov//store/mc//.root' - 'root://cmsxrootd.fnal.gov//store/mc/RunIIFall17DRPremix/QCD_Pt-15to7000_TuneCP5_Flat_13TeV_pythia8/AODSIM/94X_mc2017_realistic_v10-v1/50000/00304636-1BDB-E711-B6F3-FA163ECE02A9.root', - 'root://cmsxrootd.fnal.gov//store/mc/RunIIFall17DRPremix/QCD_Pt-15to7000_TuneCP5_Flat_13TeV_pythia8/AODSIM/94X_mc2017_realistic_v10-v1/50000/0036C92E-DFDB-E711-952A-008CFAFC05DE.root', + # 'root://cmsxrootd.fnal.gov//store/mc/RunIIFall17DRPremix/QCD_Pt-15to7000_TuneCP5_Flat_13TeV_pythia8/AODSIM/94X_mc2017_realistic_v10-v1/50000/00304636-1BDB-E711-B6F3-FA163ECE02A9.root', + # 'root://cmsxrootd.fnal.gov//store/mc/RunIIFall17DRPremix/QCD_Pt-15to7000_TuneCP5_Flat_13TeV_pythia8/AODSIM/94X_mc2017_realistic_v10-v1/50000/0036C92E-DFDB-E711-952A-008CFAFC05DE.root', + 'root://cms-xrd-global.cern.ch//store/user/kirschen/QCD_Pt-15to7000_TuneCP5_Flat2017_13TeV_pythia8/crab_pickEvents/180201_223255/0000/pickevents_1.root' + # /afs/cern.ch/user/k/kirschen/public/forJERC/forMCTruthDebugging/pickevents_NoPU.root ) process.source = cms.Source("PoolSource", fileNames = inputFiles ) diff --git a/JetUtilities/interface/ClosureMaker.hh b/JetUtilities/interface/ClosureMaker.hh index b202e41f..3b809829 100644 --- a/JetUtilities/interface/ClosureMaker.hh +++ b/JetUtilities/interface/ClosureMaker.hh @@ -89,14 +89,14 @@ public: private: bool objects_loaded, draw_guidelines; - double CMEnergy, nsigma; + double CMEnergy, nsigma; TString path, filename, outputDir, outputFilename, flavor, alg, histMet; vector algs, outputFormat; JetInfo *ji; TFile *ifile, *ofile; ObjectLoader hl; - vector h; - vector func; + vector h; + vector func; vector hClosure; TF1 *line, *linePlus, *lineMinus; vector > canvases_legends; @@ -104,6 +104,7 @@ private: VARIABLES::Variable var; TDirectoryFile *odir; HistUtil::HistogramMetric histogramMetric; + int statTh; }; #endif diff --git a/JetUtilities/interface/L2Creator.hh b/JetUtilities/interface/L2Creator.hh index 35785efc..345f5171 100644 --- a/JetUtilities/interface/L2Creator.hh +++ b/JetUtilities/interface/L2Creator.hh @@ -99,7 +99,7 @@ private: TString output, outputDir, l2calofit, l2pffit; vector formats, algs; bool l2l3, delphes; - int maxFitIter; + int maxFitIter, statTh; HistUtil::HistogramMetric histogramMetric; TFile* ofile; TFile* ifile; @@ -116,6 +116,7 @@ private: vector vabscor_eta; vector vrelcor_eta; vector vabscor_eta_spline; + float ptclip; }; #endif diff --git a/JetUtilities/src/ClosureMaker.cc b/JetUtilities/src/ClosureMaker.cc index c6e1cb24..f26db55c 100644 --- a/JetUtilities/src/ClosureMaker.cc +++ b/JetUtilities/src/ClosureMaker.cc @@ -57,6 +57,7 @@ ClosureMaker::ClosureMaker(CommandLine& cl) { histMet = cl.getValue ("histMet", "mu_h"); histogramMetric = HistUtil::getHistogramMetricType(string(histMet)); bool help = cl.getValue ("help", false); + statTh = cl.getValue ("statTh", 4); if (help) {cl.print(); return;} if (!cl.partialCheck()) return; @@ -343,7 +344,7 @@ void ClosureMaker::loopOverBins(TH2F* hvar, unsigned int iVarBin) { varBins[ibin].c_str(),varBins[ibin+1].c_str()); h.push_back(hvar->ProjectionY(name,ibin+1,ibin+1,"e")); - if (h.back()->GetEntries()>4) { + if (h.back()->GetEntries()> statTh) { if(histogramMetric==HistUtil::mu_f || histogramMetric==HistUtil::mpv) { int nbins = 50;//100; TSpectrum *spec = new TSpectrum(10); diff --git a/JetUtilities/src/L2Creator.cc b/JetUtilities/src/L2Creator.cc index fd6d234a..7becd59d 100644 --- a/JetUtilities/src/L2Creator.cc +++ b/JetUtilities/src/L2Creator.cc @@ -46,6 +46,9 @@ L2Creator::L2Creator(CommandLine& cl) { histMet = cl.getValue ("histMet", "mu_h"); histogramMetric = HistUtil::getHistogramMetricType(histMet); + ptclip = cl.getValue ("ptclip", 0.); + statTh = cl.getValue ("statTh", 4); + if (!cl.partialCheck()) return; cl.print(); } @@ -232,7 +235,7 @@ void L2Creator::loopOverEtaBins() { // only add points to the graphs if the current histo is not empty // the current setting might be a little high // - if (hrsp->GetEntries() > 400) {//hrsp->Integral()!=0) { + if (hrsp->GetEntries() > statTh) {//hrsp->Integral()!=0) { //TF1* frsp = (TF1*)hrsp->GetListOfFunctions()->Last(); //std::cout << "hrspName = " << hrsp->GetName() << ": frsp = " << frsp << std::endl; @@ -1313,7 +1316,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //For eta-dependent spline clipping int pt_limit = 70; - int pt_clip = 8; + float pt_clip = ptclip; unsigned int vector_size = 0; vector_size = vabscor_eta.size(); diff --git a/JetUtilities/src/SynchFittingProcedure.hh b/JetUtilities/src/SynchFittingProcedure.hh index a4f5cbc0..6340a317 100644 --- a/JetUtilities/src/SynchFittingProcedure.hh +++ b/JetUtilities/src/SynchFittingProcedure.hh @@ -102,7 +102,7 @@ TCanvas * getCanvasFromFittingProcedure(TString cname , TProfile2D * prof, TStri int S_comp_ = cname.CompareTo("ParametersVsGenSumPtOA"); //cout <Fit("pol2","0QS","",50,1800); } @@ -114,7 +114,7 @@ TCanvas * getCanvasFromFittingProcedure(TString cname , TProfile2D * prof, TStri { fr = hoff->Fit("pol2","0QS","",5,35); } - + // Skip if fit failed if (fr->Status()){ @@ -126,7 +126,7 @@ TCanvas * getCanvasFromFittingProcedure(TString cname , TProfile2D * prof, TStri cout<<" ERROR getCanvasFromFittingProcedure() NPAR != than fr->NPar()"<NPar();p++){ aux[p]->SetBinContent(eb,fr->Parameter(p)); @@ -147,12 +147,12 @@ TCanvas * getCanvasFromFittingProcedure(TString cname , TProfile2D * prof, TStri //clean up delete hoff; - }//for + }//for // Close the file fout->Write(); fout->Close(); - + // Put all histos into the canvas aux[0]->GetYaxis()->SetRangeUser(-8,8); aux[1]->GetYaxis()->SetRangeUser(-1.5,5); @@ -164,7 +164,7 @@ TCanvas * getCanvasFromFittingProcedure(TString cname , TProfile2D * prof, TStri } return c; - + }//getCanvasFromFittingProcedure @@ -189,12 +189,12 @@ TH1 * getResolutionHistoFromHisto(TString cname, TString title, TH2 * histo_in){ //double ptreferr = 0.5*(histo->GetXaxis()->GetBinLowEdge(nb)+histo->GetXaxis()->GetBinUpEdge(nb)); double val = 0; double valerr = 0; - + TH1 * aux = histo_in->ProjectionY("_py",nb,nb); if (aux->GetEntries() > 0) { TFitResultPtr fr = aux->Fit("gaus","0qS"); - + // Skip if fit failed if (fr.Get() && !fr->Status()){ double mean = fr->Parameter(1); @@ -239,12 +239,12 @@ TH1 * getResolutionHistoFromHisto_v3(TString cname, TString title, TH2 * histo_i //double ptreferr = 0.5*(histo->GetXaxis()->GetBinLowEdge(nb)+histo->GetXaxis()->GetBinUpEdge(nb)); double val = 0; double valerr = 0; - + TH1 * aux = histo_in->ProjectionY("_py",nb,nb); if (aux->GetEntries() > 0) { TFitResultPtr fr = aux->Fit("gaus","0qS"); - + // Skip if fit failed if (!fr->Status()){ //double mean = fr->Parameter(1); @@ -289,7 +289,7 @@ TH1 * getResolutionHistoFromHisto_v2(TString cname, TString title, TH2 * histo_i //double ptreferr = 0.5*(histo->GetXaxis()->GetBinLowEdge(nb)+histo->GetXaxis()->GetBinUpEdge(nb)); double val = 0; double valerr = 0; - + TH1 * aux = histo_in->ProjectionY("_py",nb,nb); TH1 * aux2= off_in->ProjectionY("_py",nb,nb); if (aux->GetEntries() > 0 && aux2->GetEntries()>0) { @@ -345,7 +345,7 @@ TH1 * getMeanHistoFromHisto(TString cname, TString title, TH2 *off_in, double & histo->Reset(); //histo->Clear(); histo->GetYaxis()->SetTitle(title); - + // Now loop over the entries of prof and set the histo for (int nb = 1 ; nb <= histo->GetXaxis()->GetNbins() ; nb++){ @@ -353,7 +353,7 @@ TH1 * getMeanHistoFromHisto(TString cname, TString title, TH2 *off_in, double & //double ptreferr = 0.5*(histo->GetXaxis()->GetBinLowEdge(nb)+histo->GetXaxis()->GetBinUpEdge(nb)); double val = 0; double valerr = 0; - + TH1 * aux= off_in->ProjectionY("_py",nb,nb); if (aux->GetEntries() > 0) { @@ -406,7 +406,7 @@ TH1 * getMeanOverBinCenterHistoFromHisto(TString cname, TString title, TH2 *off_ histo->Reset(); //histo->Clear(); histo->GetYaxis()->SetTitle(title); - + // Now loop over the entries of prof and set the histo for (int nb = 1 ; nb <= histo->GetXaxis()->GetNbins() ; nb++){ @@ -414,7 +414,7 @@ TH1 * getMeanOverBinCenterHistoFromHisto(TString cname, TString title, TH2 *off_ double ptreferr = 0.5*(histo->GetXaxis()->GetBinLowEdge(nb)+histo->GetXaxis()->GetBinUpEdge(nb)); double val = 0; double valerr = 0; - + TH1 * aux= off_in->ProjectionY("_py",nb,nb); if (aux->GetEntries() > 0) { @@ -503,7 +503,7 @@ TCanvas * getCanvasResolution(TString cname, TString algo, TString title, vector cout<<"\t Doing fits for Resolution "<Clear(); histod->GetYaxis()->SetTitle(""); histod->SetTitle(ctitle+" vs. p_{T}^{GEN}"); - + //---to check, produce n/d - + TH1 * histocheck = off->ProjectionX(cname); histocheck->Reset(); histocheck->Clear(); histocheck->GetYaxis()->SetTitle("#sigma(p_{T}^{PU}-p_{T}^{noPU})/"); histocheck->SetTitle(ctitle+" #sigma(p_{T}^{PU}-p_{T}^{noPU})/ vs. p_{T}^{GEN}"); - + double maxy1=0; double maxy2=0; double maxycheck=0; @@ -717,7 +717,7 @@ TCanvas * getResolutionNumDenom(TString cname, TString ctitle, TString algo, TH2 rmserr2 = fr2->ParError(2); val = rms2/mean; valerr = val * sqrt( pow(rmserr2/rms2,2) + pow(meanerr/mean,2)); - + if (maxy1cd(); ccheck->SetLogx(); histocheck->Draw("E"); - + histon->GetYaxis()->SetRangeUser(0,30); histod->GetYaxis()->SetRangeUser(0,maxy1*1.2); @@ -749,7 +749,7 @@ TCanvas * getResolutionNumDenom(TString cname, TString ctitle, TString algo, TH2 c->cd(2); histod->Draw("E"); c->cd(); - + return c; }//getResolutionNumDenom @@ -764,7 +764,7 @@ TCanvas * getGausMeanOffset(TString cname, TString ctitle, TString algo, vector< cout<<"\t Doing fits for Mean "<Draw("sameE"); tdrDraw(hh,"E",kOpenCircle,1,kSolid,1); - + TLegend* leg = (TLegend*)baseCanvas->GetPrimitive(cname+"_leg"); int NPV_Rho; if (cname.Contains("npv",TString::kIgnoreCase)) @@ -971,8 +971,8 @@ TCanvas * getGausMeanOffsetOverPtref(TString cname, TString ctitle, TString algo cout<<"\t Doing fits for Mean "<SetLogx(); @@ -1060,8 +1060,8 @@ TCanvas * getGausMeanOffsetScale(TString cname, TString ctitle, TString algo, ve cout<<"\t Doing fits for Mean "<SetLogx(); TLegend * leg = new TLegend(0.15,0.72,0.4,0.99); @@ -1085,7 +1085,7 @@ TCanvas * getGausMeanOffsetScale(TString cname, TString ctitle, TString algo, ve histo->Reset(); //histo->Clear(); histo->Sumw2(); - + double maxy = 0; double miny = 0; @@ -1094,14 +1094,14 @@ TCanvas * getGausMeanOffsetScale(TString cname, TString ctitle, TString algo, ve TString hname = cname; hname += Form("_%i",j); hh[j] = getMeanHistoFromHisto(hname, ctitle, off[j], miny, maxy); - //hh[j]->Sumw2(); + //hh[j]->Sumw2(); for (int nb = 1 ; nb <= histo->GetXaxis()->GetNbins() ; nb++) { histo->SetBinContent(nb,hh[j]->GetBinContent(scaleNo)); histo->SetBinError(nb,hh[j]->GetBinError(scaleNo)); - } + } hh[j]->Divide(histo); } @@ -1162,7 +1162,7 @@ TCanvas * getOffsetStack(TString cname, TString ctitle, TString algo, vectorGetXaxis()->SetRangeUser(10.0,300.0); hh->GetYaxis()->SetRangeUser(0,0.8); hh->Draw("sameEX0"); - + TLegend* leg = (TLegend*)baseCanvas->GetPrimitive(cname+"_leg"); int NPV_Rho; if (cname.Contains("npv",TString::kIgnoreCase)) @@ -1450,8 +1450,8 @@ TCanvas * getCanvasResolution_v2(TString cname, TString algo, TString title, vec cout<<"\t Doing fits for Resolution "<SetLogx(); TLegend * leg = new TLegend(0.2,0.56,0.45,0.85); @@ -1606,9 +1606,9 @@ TH1 * getIntegralHistoFromHisto(TString cname, TString title,TProfile *off_in){ for (int nb = 1 ; nb <= histo->GetXaxis()->GetNbins() ; nb++){ double valerr; double val = off_in->IntegralAndError(nb,binmax,valerr); - + //cout <GetXaxis()->GetBinCenter(nb)<<" "<GetBinContent(nb)<<" "<SetBinContent(nb,val); histo->SetBinError(nb,valerr); if (val>maxy) maxy = val; @@ -1635,8 +1635,8 @@ TCanvas * getCanvasAverage(TString cname, TString algo, TString title, vectorSetFillColor(0); leg->SetBorderSize(0); - //vector hh(prof.size(),(TH1*)0); - vector hh; + //vector hh(prof.size(),(TH1*)0); + vector hh; for (unsigned int j=0;jIntegralAndError(nb,binmax,inteErr); double jetEneAvg = 0; double jetEneAvgErr = 0; - + for (int nb2 = nb; nb2<=binmax;nb2++) { jetEneAvg += off_in->GetXaxis()->GetBinCenter(nb)*off_in->GetBinContent(nb); jetEneAvgErr += off_in->GetBinError(nb)*off_in->GetXaxis()->GetBinCenter(nb) + off_in->GetBinContent(nb)* 0.5*(off_in->GetXaxis()->GetBinWidth(nb)); } - + if (inte>0) { jetEneAvgErr = sqrt(inte*inte*jetEneAvgErr*jetEneAvgErr+jetEneAvg*jetEneAvg*inteErr*inteErr)/inteErr/inteErr; - + jetEneAvg /= inte; - + histo->SetBinContent(nb,jetEneAvg); histo->SetBinError(nb,jetEneAvgErr); if (jetEneAvg>maxy) maxy = jetEneAvg; @@ -1777,7 +1777,7 @@ double findNonOverlappingYmax(TCanvas* c, vector hists, TLegend* leg, bool vcmin[0] = leg->GetY1NDC(); //legend always at the bottom vcmin[1] = tbound; //just compare to top of plot (margin + ticks) on the other side - + //loop over histos double bh[2]; //height of highest bin + error (legend) bh[0] = bh[1] = 0; @@ -1792,11 +1792,11 @@ double findNonOverlappingYmax(TCanvas* c, vector hists, TLegend* leg, bool for(int ibin = 1; ibin <= hists[s]->GetNbinsX(); ibin++) { if(hists[s]->GetBinLowEdge(ibin)GetXmin() || hists[s]->GetBinLowEdge(ibin)>=x1->GetXmax()) unseenBins++; } - + for(int i = 0; i < 2; i++){ //check each side of plot //new bin #s for limited range int xbmin, xbmax; - + xbmin = logxy.first ? x1->FindBin(x1->GetXmin()*pow(x1->GetXmax()/x1->GetXmin(), (ucmin[i] - pad->GetLeftMargin())/gx)) : x1->FindBin((ucmin[i] - pad->GetLeftMargin())*(x1->GetXmax() - x1->GetXmin())/gx + x1->GetXmin()); @@ -1805,7 +1805,7 @@ double findNonOverlappingYmax(TCanvas* c, vector hists, TLegend* leg, bool ? x1->FindBin(x1->GetXmin()*pow(x1->GetXmax()/x1->GetXmin(), (ucmax[i] - pad->GetLeftMargin())/gx)) : x1->FindBin((ucmax[i] - pad->GetLeftMargin())*(x1->GetXmax() - x1->GetXmin())/gx + x1->GetXmin()); if(xbmax < hists[s]->GetNbinsX()) xbmax += 1; //include partial overlap - + if(debug) { cout << "hist name = " << hists[s]->GetName() << endl << "\tucmin[" << i << "] = " << ucmin[i] << endl From 691f9be9a38577a684acd92d82354d04a2e8eb3e Mon Sep 17 00:00:00 2001 From: Sifu Luo Date: Sat, 10 Feb 2018 19:06:19 -0600 Subject: [PATCH 03/10] Fit included in HistogramUtility --- JetAnalyzers/bin/jet_synchplot_x.cc | 208 +++++++++---------- JetUtilities/interface/HistogramUtilities.hh | 11 +- JetUtilities/src/HistogramUtilities.cc | 112 +++++++++- JetUtilities/src/SynchFittingProcedure.hh | 11 +- 4 files changed, 225 insertions(+), 117 deletions(-) diff --git a/JetAnalyzers/bin/jet_synchplot_x.cc b/JetAnalyzers/bin/jet_synchplot_x.cc index b93f72a9..6e92e107 100644 --- a/JetAnalyzers/bin/jet_synchplot_x.cc +++ b/JetAnalyzers/bin/jet_synchplot_x.cc @@ -71,7 +71,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 TString algo1(calgo1); TString algo2(calgo2); TString algo12 = algo1+"_"+algo2; - if (algo1.EqualTo(algo2)) + if (algo1.EqualTo(algo2)) algo12 = algo1; TString algo(algo12); TString filename=Form("%s/output_%s.root",inputDir.Data(),algo.Data()); @@ -114,20 +114,20 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 return; */ - // Event-Matching performance + // Event-Matching performance if(histograms.find("m_refpt_diff")!=histograms.end()) { c = new TCanvas("RefPtDiff","RefPtDiff"); c->SetLogy(); histograms["m_refpt_diff"]->Draw(); } - - // Event-Matching performance + + // Event-Matching performance if(histograms.find("m_refpdgid_diff")!=histograms.end()) { c = new TCanvas("RefPdgidDiff","RefPdgidDiff"); c->SetLogy(); histograms["m_refpdgid_diff"]->Draw(); } - + // Sanity check: g_pthat if(histograms.find("m_deltaPthat")!=histograms.end()) { c = new TCanvas("PthatDiff","PthatDiff"); @@ -146,7 +146,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 histograms["m_njet_pt_pu"]->GetYaxis()->SetRangeUser(0,4.5e6); histograms["m_njet_pt_pu"]->Draw("E"); histograms["m_njet_pt_nopu"]->Draw("sameE"); - + leg = new TLegend(0.7,0.4,0.9,0.6); leg->SetFillColor(0); leg->SetBorderSize(0); @@ -164,7 +164,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 ratio->GetYaxis()->SetRangeUser(0,10); ratio->Divide(histograms["m_njet_pt_nopu"]); ratio->Draw("E"); - + leg = new TLegend(0.7,0.4,0.9,0.6); leg->SetFillColor(0); leg->SetBorderSize(0); @@ -181,14 +181,14 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 ratio_GenPtCut->GetYaxis()->SetRangeUser(0,10); ratio_GenPtCut->Divide(histograms["m_njet_pthigh_nopu"]); ratio_GenPtCut->Draw("E"); - + leg = new TLegend(0.7,0.4,0.9,0.6); leg->SetFillColor(0); leg->SetBorderSize(0); leg->AddEntry(ratio_GenPtCut, "#frac{PU sample}{NoPU sample}","lep"); leg->Draw(); } - + // njet vs npv if(histograms.find("m_all_nj_npv")!=histograms.end() && histograms.find("m_matched_nj_npv")!=histograms.end() && @@ -201,7 +201,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 histograms["m_all_nj_npv"]->Draw("E"); histograms["m_matched_nj_npv"]->Draw("sameE"); histograms["m_unmatched_nj_npv"]->Draw("sameE"); - + leg = new TLegend(0.7,0.75,0.9,0.95); leg->SetFillColor(0); leg->SetBorderSize(0); @@ -223,7 +223,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 histograms["m_all_jtpt_npv"]->Draw("E"); histograms["m_matched_jtpt_npv"]->Draw("sameE"); histograms["m_unmatched_jtpt_npv"]->Draw("sameE"); - + leg = new TLegend(0.7,0.75,0.9,0.95); leg->SetFillColor(0); leg->SetBorderSize(0); @@ -232,7 +232,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 leg->AddEntry(histograms["m_unmatched_jtpt_npv"], "UnMatched jets","lep"); leg->Draw(); } - + // Fraction of Matched Jets c = new TCanvas("FractionMatchedJetsNoPU","FractionMatchedJets NoPU Sample"); @@ -279,7 +279,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 } leg->Draw(); - + // Fraction of RG-Matched Jets c = new TCanvas("FractionRGMatchedJetsNoPU","FractionRGMatchedJets NoPU Sample"); c->SetLogx(); @@ -300,8 +300,8 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 histograms[hname]->GetYaxis()->SetRangeUser(0.5,1.1); leg->AddEntry(histograms[hname],detector_regions[det],"lep"); } - leg->Draw(); - + leg->Draw(); + // Fraction of RG-Matched Jets c = new TCanvas("FractionRGMatchedJetsPU","FractionRGMatchedJets PU Sample"); c->SetLogx(); @@ -326,9 +326,9 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 leg->Draw(); // Fraction of Matched Jets-PU-NPV - if(histograms.find("m_frac_nj_pt_b_match_pu_npv10")!=histograms.end() && - histograms.find("m_frac_nj_pt_b_match_pu_npv20")!=histograms.end() && - histograms.find("m_frac_nj_pt_b_match_pu_npv30")!=histograms.end() && + if(histograms.find("m_frac_nj_pt_b_match_pu_npv10")!=histograms.end() && + histograms.find("m_frac_nj_pt_b_match_pu_npv20")!=histograms.end() && + histograms.find("m_frac_nj_pt_b_match_pu_npv30")!=histograms.end() && histograms.find("m_frac_nj_pt_b_match_pu_npvO")!=histograms.end()) { c = new TCanvas("FractionMatchedJetsPU_NPV","FractionMatchedJets vs. npv PU Sample"); c->SetLogx(); @@ -342,7 +342,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 histograms["m_frac_nj_pt_b_match_pu_npv20"]->Draw("same"); histograms["m_frac_nj_pt_b_match_pu_npv30"]->Draw("same"); histograms["m_frac_nj_pt_b_match_pu_npvO"]->Draw("same"); - + leg = new TLegend(0.7,0.4,0.9,0.6); leg->SetFillColor(0); leg->SetBorderSize(0); @@ -369,7 +369,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 histograms["m_frac_nj_pt_b_match_RG_pu_npv20"]->Draw("same"); histograms["m_frac_nj_pt_b_match_RG_pu_npv30"]->Draw("same"); histograms["m_frac_nj_pt_b_match_RG_pu_npvO"]->Draw("same"); - + leg = new TLegend(0.7,0.4,0.9,0.6); leg->SetFillColor(0); leg->SetBorderSize(0); @@ -396,7 +396,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 histograms["m_frac_nj_pt_b_match_nopu_npv2"]->Draw("same"); histograms["m_frac_nj_pt_b_match_nopu_npv3"]->Draw("same"); histograms["m_frac_nj_pt_b_match_nopu_npvO"]->Draw("same"); - + leg = new TLegend(0.7,0.4,0.9,0.6); leg->SetFillColor(0); leg->SetBorderSize(0); @@ -434,7 +434,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 // INTEGRAL of PU distribution of jets per event per NPV // Number of unmatched jets in pu sample minus #of unmatch jets in nopu sample - + vector pnJetPt(npvRhoNpuBins.size(),(TProfile*)0); for(unsigned int ibin=0; ibinSetLogx(); c->Draw(); - + c = getCanvasAverage("AvgPUJetEneDistribution",algo,"Avg. Jet Ene (GeV)",pnJetPt,npvRhoNpuBins); c->SetLogx(); c->Draw(); - + // - // 2D histogram of jtarea diff. vs. refpt + // 2D histogram of jtarea diff. vs. refpt // if(histograms.find("p_areaVsrefpt")!=histograms.end()) { c = new TCanvas("areaVsrefpt","areaVsrefpt"); histograms["p_areaVsrefpt"]->Draw("colZ"); - c->SetLogx(); - + c->SetLogx(); + TProfile *p_areaVsrefpt_prof = dynamic_cast(histograms["p_areaVsrefpt"])->ProfileX(); c = new TCanvas("areaVsrefptProf","areaVsrefptProf"); p_areaVsrefpt_prof->GetYaxis()->SetTitle(""); @@ -466,7 +466,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 } // - // 2D histogram of jtarea diff. vs. refpt + // 2D histogram of jtarea diff. vs. refpt // if(histograms.find("p_areaVsoffset_1000")!=histograms.end()) { c = new TCanvas("areaVsOffset_1000","areaVsOffset_1000"); @@ -479,10 +479,10 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 p_areaVsoffset_1000_prof->GetYaxis()->SetRangeUser(-0.3,0.3); p_areaVsoffset_1000_prof->Draw(""); } - - + + // - // 2D histogram of jtarea diff. vs. refpt + // 2D histogram of jtarea diff. vs. refpt // if(histograms.find("p_areaVsoffset_30_50")!=histograms.end()) { c = new TCanvas("areaVsOffset_30_50","areaVsOffset_30_50"); @@ -496,22 +496,22 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 p_areaVsoffset_30_50_prof->Draw(""); } - - // + + // // 2D histogram of NPV vs Rho with 15Draw("COLZ"); + c = new TCanvas("npvVsRhoOffset1515h","npvVsRhoOffset1515h"); + histograms["p_npvVsRho_offset_15_15h"]->Draw("COLZ"); } - - - // + + + // // 2D histogram of NPV vs Rho with 15Draw("COLZ"); + c = new TCanvas("npvVsRhoOffset1515hPeak","npvVsRhoOffset1515hPeak"); + histograms["p_npvVsRho_offset_15_15h"]->Draw("COLZ"); //p_npvVsRho_offset_15_15h->ShowPeaks(2,"nodraw",0.2); //TList *functions = p_npvVsRho_offset_15_15h->GetListOfFunctions(); //TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker"); @@ -529,26 +529,26 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 { p_npvVsRho_offset_15_15h_y->SetBinContent(it,histograms["p_npvVsRho_offset_15_15h"]->GetBinContent((int)peakX,it)); } - c = new TCanvas("npvVsRhoOffset1515hX","npvVsRhoOffset1515hX"); + c = new TCanvas("npvVsRhoOffset1515hX","npvVsRhoOffset1515hX"); p_npvVsRho_offset_15_15h_x->Draw(); p_npvVsRho_offset_15_15h_x->Fit("gaus"); - c = new TCanvas("npvVsRhoOffset1515hY","npvVsRhoOffset1515hY"); + c = new TCanvas("npvVsRhoOffset1515hY","npvVsRhoOffset1515hY"); p_npvVsRho_offset_15_15h_y->Draw(); p_npvVsRho_offset_15_15h_y->Fit("gaus"); } - + // - // Profile of dr vs. refpt for the matched jets + // Profile of dr vs. refpt for the matched jets // if(histograms.find("p_drVsrefpt")!=histograms.end()) { - c = new TCanvas("drVsrefptMatchedJets","drVsrefptMatchedJets"); + c = new TCanvas("drVsrefptMatchedJets","drVsrefptMatchedJets"); dynamic_cast(histograms["p_drVsrefpt"])->SetErrorOption("s"); c->SetLogx(); histograms["p_drVsrefpt"]->Draw(); histograms["p_drVsrefpt"]->GetYaxis()->SetRangeUser(0,0.3); histograms["p_drVsrefpt"]->GetYaxis()->SetTitle("<#DeltaR> #pm #sigma(#DeltaR)"); } - + // // Profile of npv and rho vs offset PU // @@ -556,15 +556,15 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 c = new TCanvas("npvrhoVsOffset","npvrhoVsOffset"); TProfile *p_npvVsOff_prof = dynamic_cast(histograms["p_npvVsOff"])->ProfileX(); TProfile *p_rhoVsOff_prof = dynamic_cast(histograms["p_rhoVsOff"])->ProfileX(); - p_npvVsOff_prof->SetErrorOption("s"); + p_npvVsOff_prof->SetErrorOption("s"); setHistoColor(p_npvVsOff_prof,colPU); p_npvVsOff_prof->Draw("E1"); p_npvVsOff_prof->GetYaxis()->SetRangeUser(0,45); p_npvVsOff_prof->GetYaxis()->SetTitle(" #pm #sigma(X)"); setHistoColor(p_rhoVsOff_prof,colNoPU); - p_rhoVsOff_prof->SetErrorOption("s"); + p_rhoVsOff_prof->SetErrorOption("s"); p_rhoVsOff_prof->Draw("sameE1"); - + leg = new TLegend(0.2,0.72,0.4,0.92); leg->SetFillColor(0); leg->SetBorderSize(0); @@ -572,7 +572,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 leg->AddEntry(p_rhoVsOff_prof," Rho ","lep"); leg->Draw(); } - + // Profiles of npv, rho, tnpu vs themselves if(histograms.find("p_npvVsNpv")!=histograms.end() && histograms.find("p_rhoVsRho")!=histograms.end() && @@ -630,7 +630,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 // 2D plot of npv vs offset PU // if(histograms.find("p_npvVsOff")!=histograms.end()) { - c = new TCanvas("npvVsOffset2D","npvVsOffset2D"); + c = new TCanvas("npvVsOffset2D","npvVsOffset2D"); histograms["p_npvVsOff"]->SetTitle("2D Histogram of N_{PV} and _{jets}, LogZ"); histograms["p_npvVsOff"]->GetYaxis()->SetRangeUser(0,45); //p_npvVsOff->GetXaxis()->SetRangeUser(0,45); @@ -639,12 +639,12 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 histograms["p_npvVsOff"]->Draw("CONTZ"); c->SetLogz(); } - + // // 2D plot of rho vs offset PU // if(histograms.find("p_rhoVsOff")!=histograms.end()) { - c = new TCanvas("rhoVsOffset2D","rhoVsOffset2D"); + c = new TCanvas("rhoVsOffset2D","rhoVsOffset2D"); histograms["p_rhoVsOff"]->SetTitle("2D Histogram of N_{PV} and _{jets}, LogZ"); histograms["p_rhoVsOff"]->GetYaxis()->SetRangeUser(0,45); //p_rhoVsOff->GetXaxis()->SetRangeUser(0,45); @@ -653,11 +653,11 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 histograms["p_rhoVsOff"]->Draw("CONTZ"); c->SetLogz(); } - + // // Profile of npv vs offset PU breakdown into detector parts // - c = new TCanvas("npvVsOffset","npvVsOffset"); + c = new TCanvas("npvVsOffset","npvVsOffset"); TProfile * hnpvOff_prof[NDetectorNames]; for (int det=0;det(histograms[hname])->ProfileX(); - hnpvOff_prof[det]->SetErrorOption("s"); + hnpvOff_prof[det]->SetErrorOption("s"); setHistoColor(hnpvOff_prof[det],colDet[det]); hnpvOff_prof[det]->GetYaxis()->SetRangeUser(0,45); } @@ -688,7 +688,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 // // Profile of rho vs offset PU breakdown into detector parts // - c = new TCanvas("rhoVsOffset","rhoVsOffset"); + c = new TCanvas("rhoVsOffset","rhoVsOffset"); TProfile * hrhoOff_prof[4]; for (int det=0;det(histograms[hname])->ProfileX(); - hrhoOff_prof[det]->SetErrorOption("s"); + hrhoOff_prof[det]->SetErrorOption("s"); setHistoColor(hrhoOff_prof[det],colDet[det]); hrhoOff_prof[det]->GetYaxis()->SetRangeUser(0,45); } @@ -714,7 +714,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 leg->AddEntry(hrhoOff_prof[det],detector_names[det],"lep"); } leg->Draw(); - + // // Jet Energy Resolution (sigma(pt/ptref)/mean(pt/ptref) vs. ptref) PU breakdown into detector parts // @@ -727,32 +727,32 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 } c = getCanvasResponseResolution("PUresponseResolutionVsptref",algo, "#sigma(p_{T}^{PU}/p_{T}^{GEN})/",hresResPt); c->Draw(); - + for(int det=0; det(histograms[hname]); - } + } c = getCanvasResponseResolution("NoPUresponseResolutionVsptref",algo, "#sigma(p_{T}^{noPU}/p_{T}^{GEN})/",hresResPt); c->Draw(); - - + + // // profile # of matchedjet vs offset // - if(histograms.find("p_matchedjet_off")!=histograms.end()) { - c = new TCanvas("MatchedJetOffset","MatchedJetOffset"); - dynamic_cast(histograms["p_matchedjet_off"])->SetErrorOption("s"); + if(histograms.find("p_matchedjet_off")!=histograms.end()) { + c = new TCanvas("MatchedJetOffset","MatchedJetOffset"); + dynamic_cast(histograms["p_matchedjet_off"])->SetErrorOption("s"); histograms["p_matchedjet_off"]->Draw(); } // Offset PT energy distribution, constructed from pt(pu)-pt(nopu) VS NPV - if(histograms.find("p_off_etaVsNpv")!=histograms.end()) { + if(histograms.find("p_off_etaVsNpv")!=histograms.end()) { c = new TCanvas("OffsetDistributionVsNPV","OffsetDistributionVsNPV"); histograms["p_off_etaVsNpv"]->GetYaxis()->SetNdivisions(6); histograms["p_off_etaVsNpv"]->Draw("lego2"); @@ -760,22 +760,22 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 } // Offset PT energy distribution, constructed from pt(pu)-pt(nopu) VS Rho - if(histograms.find("p_off_etaVsRho")!=histograms.end()) { + if(histograms.find("p_off_etaVsRho")!=histograms.end()) { c = new TCanvas("OffsetDistributionVsRho","OffsetDistributionVsRho"); histograms["p_off_etaVsRho"]->GetYaxis()->SetNdivisions(6); histograms["p_off_etaVsRho"]->Draw("lego2"); histograms["p_off_etaVsRho"]->GetZaxis()->SetRangeUser(-10,50); } - + // Offset PT energy distribution, constructed from pt(pu)-pt(nopu) VS PUEff - if(histograms.find("p_off_etaVsPUEff")!=histograms.end()) { + if(histograms.find("p_off_etaVsPUEff")!=histograms.end()) { c = new TCanvas("OffsetDistributionVsPUEff","OffsetDistributionVsPUEff"); histograms["p_off_etaVsPUEff"]->Draw("lego2"); histograms["p_off_etaVsPUEff"]->GetZaxis()->SetRangeUser(-10,50); } // Offset PT energy distribution, constructed from pt(pu)-pt(nopu) VS GenSumPtOA - if(histograms.find("p_off_etaVsGenSumPtOA")!=histograms.end()) { + if(histograms.find("p_off_etaVsGenSumPtOA")!=histograms.end()) { c = new TCanvas("OffsetDistributionVsGenSumPtOA","OffsetDistributionVsGenSumPtOA"); histograms["p_off_etaVsGenSumPtOA"]->GetYaxis()->SetTitle("SumPt/jetArea"); histograms["p_off_etaVsGenSumPtOA"]->GetYaxis()->SetNdivisions(6); @@ -785,7 +785,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 } // Offset PT energy distribution, constructed from pt(pu)-pt(nopu) VS JetPt - if(histograms.find("p_off_etaVsJetPt")!=histograms.end()) { + if(histograms.find("p_off_etaVsJetPt")!=histograms.end()) { c = new TCanvas("OffsetDistributionVsJetPt","OffsetDistributionVsJetPt"); histograms["p_off_etaVsJetPt"]->Draw("lego2"); c->SetLogy(); @@ -793,14 +793,14 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 } // OffsetOverArea PT energy distribution, constructed from (pt(pu)-pt(nopu))/area(pu) VS JetPt - if(histograms.find("p_offOverA_etaVsJetPt")!=histograms.end()) { + if(histograms.find("p_offOverA_etaVsJetPt")!=histograms.end()) { c = new TCanvas("OffsetOverAreaDistributionVsJetPt","OffsetOverAreaDistributionVsJetPt"); histograms["p_offOverA_etaVsJetPt"]->Draw("lego2"); c->SetLogy(); histograms["p_offOverA_etaVsJetPt"]->GetZaxis()->SetRangeUser(-40,100); } - // do the fitting in each eta range and return the parameters. + // do the fitting in each eta range and return the parameters. // the last parameter is the name of the file name with which all functions are saved to. c = getCanvasFromFittingProcedure("ParametersVsNpv",dynamic_cast(histograms["p_off_etaVsNpv"]),"fittingFunctionsNpv_"+algo+".root"); if(c) @@ -812,23 +812,23 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 if(c) c->Draw(); fin->cd(); - - // do the fitting in each eta range and return the parameters. + + // do the fitting in each eta range and return the parameters. // the last parameter is the name of the file name with which all functions are saved to. c = getCanvasFromFittingProcedure("ParametersVsPUEff",dynamic_cast(histograms["p_off_etaVsPUEff"]),"fittingFunctionsPUEff_"+algo+".root"); if(c) c->Draw(); fin->cd(); - - - // do the fitting in each eta range and return the parameters. + + + // do the fitting in each eta range and return the parameters. // the last parameter is the name of the file name with which all functions are saved to. c = getCanvasFromFittingProcedure("ParametersVsGenSumPtOA",dynamic_cast(histograms["p_off_etaVsGenSumPtOA"]),"fittingFunctionsGenSumPtOA_"+algo+".root"); if(c) c->Draw(); fin->cd(); - // do the fitting in each eta range and return the parameters. + // do the fitting in each eta range and return the parameters. // the last parameter is the name of the file name with which all functions are saved to. c = getCanvasFromFittingProcedure("ParametersOffOverAVsJetPt",dynamic_cast(histograms["p_offOverA_etaVsJetPt"]),"fittingFunctionsOffOverAJetPt_"+algo+".root"); // c->SetLogy(); @@ -1036,7 +1036,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 pOffPF[ipf] = dynamic_cast(histograms[hname]); } c = getCanvasResolution("ResolutionOffRefPF_BB",algo,"#sigma(p_{T}^{PU}-p_{T}^{noPU})",hResRho,1,npvRhoNpuBins); - c->Draw(); + c->Draw(); c = getGausMeanOffset("MeanOffRefPF_BB","",algo,hResRho,fixedRange,npvRhoNpuBins); c->Draw(); c = getGausMeanOffsetWithSum("MeanOffRefPFWithSum_BB","",algo,hResRho,dynamic_cast(histograms["p_offResVsrefpt_bb_all"]),fixedRange,npvRhoNpuBins,make_pair(minNpvRhoNpu,maxNpvRhoNpu)); @@ -1196,7 +1196,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 c = getGausMeanOffsetOverPtref("OffMeanOverPttnpuRef_EI","/p_{T}^{GEN}",algo,hOffRho,fixedRange,npvRhoNpuBins); c->Draw(); clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,npvRhoNpuBins.size()); - + for(unsigned int ibin=0; ibin(histograms[hname]); @@ -1209,8 +1209,8 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 c->Draw(); c = getGausMeanOffsetWithSum("OffMeanrhoRefWithSum_EO"," (GeV)",algo,hOffRho,dynamic_cast(histograms[Form("p_offresVsrefpt_eo_rho%i_%i",minNpvRhoNpu,maxNpvRhoNpu)]),fixedRange,npvRhoNpuBins,make_pair(minNpvRhoNpu,maxNpvRhoNpu)); c->Draw(); - clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,npvRhoNpuBins.size()); - + clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,npvRhoNpuBins.size()); + for(unsigned int ibin=0; ibin(histograms[hname]); @@ -1230,7 +1230,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 hOffRho[ibin] = dynamic_cast(histograms[hname]); } c = getCanvasResolution_v2("OffResolutionRhoRef_FF",algo,"#sigma(p_{T}^{PU}-p_{T}^{noPU})/",hResRho,hOffRho,npvRhoNpuBins); - c->Draw(); + c->Draw(); c = getGausMeanOffset("OffMeanrhoRef_FF"," (GeV)",algo,hOffRho,fixedRange,npvRhoNpuBins); c->Draw(); c = getGausMeanOffsetWithSum("OffMeanrhoRefWithSum_FF"," (GeV)",algo,hOffRho,dynamic_cast(histograms[Form("p_offresVsrefpt_ff_rho%i_%i",minNpvRhoNpu,maxNpvRhoNpu)]),fixedRange,npvRhoNpuBins,make_pair(minNpvRhoNpu,maxNpvRhoNpu)); @@ -1248,7 +1248,7 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 c = getGausMeanOffsetOverPtref("OffMeanOverPttnpuRef_FF","/p_{T}^{GEN}",algo,hOffRho,fixedRange,npvRhoNpuBins); c->Draw(); clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,npvRhoNpuBins.size()); - + //Resolution of response for PU // get the canvas from the resolution for bb @@ -1306,30 +1306,30 @@ void SynchPlots(TString inputDir="./",TString calgo1="ak5pf",TString calgo2="ak5 c = getGausMeanOffsetScale("OffMeannpvRef_BB_3035","/",algo,hOffRho,binNum3035,fixedRange,npvRhoNpuBins); c->Draw(); - + c = getGausMeanOffsetScale("OffMeannpvRef_BB_2023","/",algo,hOffRho,binNum2023,fixedRange,npvRhoNpuBins); c->Draw(); clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,npvRhoNpuBins.size()); - + for(unsigned int ibin=0; ibin(histograms[hname]); hname = Form("p_offresVsrefpt_ei_npv%i_%i",npvRhoNpuBins[ibin].first,npvRhoNpuBins[ibin].second); hOffRho[ibin] = dynamic_cast(histograms[hname]); } - + c = getCanvasResolution_v2("OffResolutionnpvRef_EI",algo,"#sigma(p_{T}^{PU}-p_{T}^{noPU})/",hResRho,hOffRho,npvRhoNpuBins); c->Draw(); - + c = getGausMeanOffset("OffMeannpvRef_EI"," (GeV)",algo,hOffRho,fixedRange,npvRhoNpuBins); c->Draw(); - + c = getGausMeanOffsetScale("OffMeannpvRef_EI_3035","/",algo,hOffRho,binNum3035,fixedRange,npvRhoNpuBins); c->Draw(); - + c = getGausMeanOffsetScale("OffMeannpvRef_EI_2023","/",algo,hOffRho,binNum2023,fixedRange,npvRhoNpuBins); c->Draw(); - clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,npvRhoNpuBins.size()); + clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,npvRhoNpuBins.size()); for(unsigned int ibin=0; ibin(histograms[hname]); } - + c = getCanvasResolution_v2("OffResolutionnpvRef_EO",algo,"#sigma(p_{T}^{PU}-p_{T}^{noPU})/",hResRho,hOffRho,npvRhoNpuBins); c->Draw(); c = getGausMeanOffset("OffMeannpvRef_EO"," (GeV)",algo,hOffRho,fixedRange,npvRhoNpuBins); - c->Draw(); - + c->Draw(); + c = getGausMeanOffsetScale("OffMeannpvRef_EO_3035","/",algo,hOffRho,binNum3035,fixedRange,npvRhoNpuBins); c->Draw(); - + c = getGausMeanOffsetScale("OffMeannpvRef_EO_2023","/",algo,hOffRho,binNum2023,fixedRange,npvRhoNpuBins); c->Draw(); - clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,npvRhoNpuBins.size()); - + clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,npvRhoNpuBins.size()); + for(unsigned int ibin=0; ibin",hResRho,hOffRho,npvRhoNpuBins); - c->Draw(); + c->Draw(); c = getGausMeanOffset("OffMeannpvRef_FF"," (GeV)",algo,hOffRho,fixedRange,npvRhoNpuBins); c->Draw(); - + c = getGausMeanOffsetScale("OffMeannpvRef_FF_3035","/",algo,hOffRho,binNum3035,fixedRange,npvRhoNpuBins); c->Draw(); - + c = getGausMeanOffsetScale("OffMeannpvRef_FF_2023","/",algo,hOffRho,binNum2023,fixedRange,npvRhoNpuBins); c->Draw(); - clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,NPDGIDcat); + clearHistograms(hResRho,hOffRho,hOffPdgid,pOffPF,NPDGIDcat); for(int ibin=0; ibin @@ -48,12 +49,12 @@ namespace HistUtil{ // (2) The mean (mu_f or mpv) and RMS (RMS_f or sigma_f) of a fit to the histogram (i.e. a Gaussian fit) // (3) The median of the histogram enum HistogramMetric {none, mu_h, RMS_h, mu_f, mpv, RMS_f, sigma_f, median}; - static const unsigned int nHistogramMetric = 8; + static const unsigned int nHistogramMetric = 8; enum AxisDirection {X, Y, Z}; static const unsigned int nAxisDirections = 3; - // A routine that returns the string given the HistogramMetric + // A routine that returns the string given the HistogramMetric string getHistogramMetricString(HistogramMetric); // A routine that returns the HistogramMetric type given the string @@ -64,7 +65,7 @@ namespace HistUtil{ // A routine that returns a given HistogramMetric and its error (if any) pair getHistogramMetric1D(HistogramMetric, TH1*, double fallback_threshold = 0.05, bool verbose = false); - pair getHistogramMetric1D(string, TH1*, double fallback_threshold = 0.05, bool verbose = false); + pair getHistogramMetric1D(string, TH1*, double fallback_threshold = 0.05, bool verbose = false); // A routine that returns the median of a given histogram and an error equal to the width of the bin that contains the histogram pair getHistogramMedian1D(TH1*, bool debug = false); @@ -91,6 +92,10 @@ namespace HistUtil{ // (i.e. if the AxisDirection=Y, then the metric will be plotted along X) vector > getHistogramMetric2D(HistogramMetric, TH2*, AxisDirection, const vector&); vector > getHistogramMetric2D(string, TH2*, AxisDirection, const vector&); + + void adjust_fitrange(TH1* h,double& min,double& max); + int number_filled_bins(TH1* h, double min, double max); + TF1* fit_gaussian(TH1*& hrsp, const double nsigma, const double jtptmin, const int niter, const int verbose); } #endif diff --git a/JetUtilities/src/HistogramUtilities.cc b/JetUtilities/src/HistogramUtilities.cc index bd4c23a4..224ac72a 100644 --- a/JetUtilities/src/HistogramUtilities.cc +++ b/JetUtilities/src/HistogramUtilities.cc @@ -14,7 +14,7 @@ namespace HistUtil { //______________________________________________________________________________ - // A routine that returns the string given the HistogramMetric + // A routine that returns the string given the HistogramMetric string getHistogramMetricString(HistogramMetric type){ if (type == mu_h) return "mu_h"; else if (type == RMS_h) return "RMS_h"; @@ -134,7 +134,7 @@ namespace HistUtil { << "\tabs(high_edge-median): " << std::abs(high_edge-median) << endl << "\tmedian-low_edge: " << median-low_edge << endl << "\thigh_edge-median: " << high_edge-median << endl; - } + } delete [] x; delete [] y; delete [] b; @@ -228,14 +228,14 @@ namespace HistUtil { // c: Median bar width if(debug) { cout << "L_m = " << boundaries[ind[jl]] << endl - << "N = " << 2.*sumTot2 << endl + << "N = " << 2.*sumTot2 << endl << "F_(m-1) = " << sumToJl << endl << "f_m = " << w[ind[jl]] << endl << "c = " << (2.*(a[ind[jl]]-boundaries[ind[jl]])) << endl; } if (n%2 == 1) median = boundaries[ind[jl]]+(((((2.*sumTot2)+1.)/2.)-sumToJl)/(w[ind[jl]]))*(2.*(a[ind[jl]]-boundaries[ind[jl]])); - else + else median = boundaries[ind[jl]]+(((sumTot2)-sumToJl)/(w[ind[jl]]))*(2.*(a[ind[jl]]-boundaries[ind[jl]])); } else { // Traditional ROOT method @@ -314,7 +314,7 @@ namespace HistUtil { if(!use_names) delete h1; - } + } return ret; }//getHistogramMetric2D @@ -323,4 +323,106 @@ namespace HistUtil { return getHistogramMetric2D(getHistogramMetricType(type),hist,axis,names); }//getHistogramMetric2D + void adjust_fitrange(TH1* h,double& min,double& max) + { + int imin=1; while (h->GetBinLowEdge(imin)GetBinLowEdge(imax)1) {imin--; min = h->GetBinCenter(imin); } + if (imaxGetNbinsX()-1) { imax++; max=h->GetBinCenter(imax); } + } + }//adjust fitting range + + int number_filled_bins(TH1* h, double min, double max) + { + int first_bin = h->FindBin(min); + int last_bin = h->FindBin(max); + int counter(0); + + for(int ibin=first_bin; ibin<=last_bin; ibin++) { + if(h->GetBinContent(ibin)>0) counter++; + } + + return counter; + } + + TF1* fit_gaussian(TH1*& hrsp, const double nsigma, const double jtptmin, const int niter, const int verbose) + { + if (0==hrsp) { + cout<<"ERROR: Empty pointer to fit_gaussian()"<GetName(); + double mean = hrsp->GetMean(); + double rms = hrsp->GetRMS(); + + double norm = hrsp->GetMaximumStored(); + double peak = mean; + int nbins = 50;//100; + TSpectrum *spec = new TSpectrum(8); + if(nbins < 100) spec->Search(hrsp,6,"nobackground nodraw goff"); //turn off background removal when nbins too small + else spec->Search(hrsp,6,"nodraw goff"); + Double_t* xpos = spec->GetPositionX(); + //Double_t* ypos = spec->GetPositionY(); + Double_t p = xpos[0]; + //Double_t ph = ypos[0]; + //std::cout << "peak: " << p << std::endl; + //std::cout << "peak height: " << ph << std::endl; + peak = p; + + double sigma = rms; + + double xmin = hrsp->GetXaxis()->GetXmin(); + double xmax = hrsp->GetXaxis()->GetXmax(); + TF1* fitfnc(0); int fitstatus(-1); + for (int iiter=0;iiter vv; + vv.push_back(xmin); + vv.push_back(peak-nsigma*sigma); + double fitrange_min = *std::max_element(vv.begin(),vv.end()); + double fitrange_max = std::min(xmax,peak+nsigma*sigma); + adjust_fitrange(hrsp,fitrange_min,fitrange_max); + fitfnc = new TF1("fgaus","gaus",fitrange_min,fitrange_max); + fitfnc->SetParNames("N","#mu","#sigma"); + fitfnc->SetParameter(0,norm); + fitfnc->SetParameter(1,peak); + fitfnc->SetParameter(2,sigma); + + //Total entries in hrsp could be >0, but the number of entries in the fit range <=1. + //In that case ROOT would print "Warning in : Fit data is empty" + //So instead check that the number of entries (unweighted integral) in the fit range is >1 + double entries_in_range = hrsp->Integral(hrsp->FindBin(fitrange_min),hrsp->FindBin(fitrange_max)); + int filled_bins_in_range = number_filled_bins(hrsp,fitrange_min,fitrange_max); + if(entries_in_range<=1 || filled_bins_in_range<=1) { + if(verbose>0) + cout << "Warning in : Fit data is empty" << endl + << "hrsp->GetName(): " << histname << endl; + return 0; + } + else { + fitstatus = hrsp->Fit(fitfnc,"RQ0"); + } + delete fitfnc; + fitfnc = hrsp->GetFunction("fgaus"); + //fitfnc->ResetBit(TF1::kNotDraw); + if (fitfnc) { + norm = fitfnc->GetParameter(0); + peak = fitfnc->GetParameter(1); + sigma = fitfnc->GetParameter(2); + } + } + if(hrsp->GetFunction("fgaus")==0) + { + cout << "No function recorded in histogram " << hrsp->GetName() << endl; + } + if (0!=fitstatus){ + cout<<"fit_gaussian() to "<GetName() + <<" failed. Fitstatus: "<GetListOfFunctions()->Delete(); + return 0; + } + TF1* returnf = dynamic_cast(hrsp->GetListOfFunctions()->Last()); + return returnf; + } }// end of namespace HistUtil diff --git a/JetUtilities/src/SynchFittingProcedure.hh b/JetUtilities/src/SynchFittingProcedure.hh index 6340a317..50d45d74 100644 --- a/JetUtilities/src/SynchFittingProcedure.hh +++ b/JetUtilities/src/SynchFittingProcedure.hh @@ -12,6 +12,7 @@ #include #include "JetMETAnalysis/JetUtilities/interface/Style.h" +#include "JetMETAnalysis/JetUtilities/interface/HistogramUtilities.hh" using namespace std; @@ -358,13 +359,14 @@ TH1 * getMeanHistoFromHisto(TString cname, TString title, TH2 *off_in, double & TH1 * aux= off_in->ProjectionY("_py",nb,nb); if (aux->GetEntries() > 0) { - TFitResultPtr fr = aux->Fit("gaus","0qS"); + // TFitResultPtr fr = aux->Fit("gaus","0qS"); //cout << cname << "sfsg1\tnb=" << nb << endl; + TF1* fr = HistUtil::fit_gaussian(aux, 1.5, 1.0, 10, false); // Skip if fit failed - if (fr.Get() && !fr->Status()){ - double mean = fr->Parameter(1); - double meanerr = fr->ParError(1); + if (fr){ + double mean = fr->GetParameter(1); + double meanerr = fr->GetParError(1); //double rms = fr->Parameter(2); //double rmserr = fr->ParError(2); @@ -372,7 +374,6 @@ TH1 * getMeanHistoFromHisto(TString cname, TString title, TH2 *off_in, double & if (val>maxy && meanerr<0.4) maxy=val; if (val Date: Mon, 12 Feb 2018 23:06:05 -0600 Subject: [PATCH 04/10] Update run_JRA_cfg.py --- JetAnalyzers/test/run_JRA_cfg.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/JetAnalyzers/test/run_JRA_cfg.py b/JetAnalyzers/test/run_JRA_cfg.py index 340a2db4..ac82a415 100644 --- a/JetAnalyzers/test/run_JRA_cfg.py +++ b/JetAnalyzers/test/run_JRA_cfg.py @@ -77,10 +77,6 @@ print "Couldn't open the external list of files from DAS. If you just checkout out the JetResponseAnalyzer package you will need to make this file yourself. Currently Falling back to opening the list hard-coded in run_JRA_cfg.py. This is not a bad action as long as it is what you intended to have happen." inputFiles = cms.untracked.vstring( # 'root://cmsxrootd.fnal.gov//store/mc//.root' - # 'root://cmsxrootd.fnal.gov//store/mc/RunIIFall17DRPremix/QCD_Pt-15to7000_TuneCP5_Flat_13TeV_pythia8/AODSIM/94X_mc2017_realistic_v10-v1/50000/00304636-1BDB-E711-B6F3-FA163ECE02A9.root', - # 'root://cmsxrootd.fnal.gov//store/mc/RunIIFall17DRPremix/QCD_Pt-15to7000_TuneCP5_Flat_13TeV_pythia8/AODSIM/94X_mc2017_realistic_v10-v1/50000/0036C92E-DFDB-E711-952A-008CFAFC05DE.root', - 'root://cms-xrd-global.cern.ch//store/user/kirschen/QCD_Pt-15to7000_TuneCP5_Flat2017_13TeV_pythia8/crab_pickEvents/180201_223255/0000/pickevents_1.root' - # /afs/cern.ch/user/k/kirschen/public/forJERC/forMCTruthDebugging/pickevents_NoPU.root ) process.source = cms.Source("PoolSource", fileNames = inputFiles ) From 4e7ddc37b2571d8bd4f85e1989c47c8ecb965510 Mon Sep 17 00:00:00 2001 From: Sifu Luo Date: Mon, 12 Feb 2018 23:11:27 -0600 Subject: [PATCH 05/10] reversed run_JRA_cfg.py to default --- JetAnalyzers/test/run_JRA_cfg.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/JetAnalyzers/test/run_JRA_cfg.py b/JetAnalyzers/test/run_JRA_cfg.py index 340a2db4..9a2a4e1f 100644 --- a/JetAnalyzers/test/run_JRA_cfg.py +++ b/JetAnalyzers/test/run_JRA_cfg.py @@ -76,11 +76,7 @@ except ImportError: print "Couldn't open the external list of files from DAS. If you just checkout out the JetResponseAnalyzer package you will need to make this file yourself. Currently Falling back to opening the list hard-coded in run_JRA_cfg.py. This is not a bad action as long as it is what you intended to have happen." inputFiles = cms.untracked.vstring( - # 'root://cmsxrootd.fnal.gov//store/mc//.root' - # 'root://cmsxrootd.fnal.gov//store/mc/RunIIFall17DRPremix/QCD_Pt-15to7000_TuneCP5_Flat_13TeV_pythia8/AODSIM/94X_mc2017_realistic_v10-v1/50000/00304636-1BDB-E711-B6F3-FA163ECE02A9.root', - # 'root://cmsxrootd.fnal.gov//store/mc/RunIIFall17DRPremix/QCD_Pt-15to7000_TuneCP5_Flat_13TeV_pythia8/AODSIM/94X_mc2017_realistic_v10-v1/50000/0036C92E-DFDB-E711-952A-008CFAFC05DE.root', - 'root://cms-xrd-global.cern.ch//store/user/kirschen/QCD_Pt-15to7000_TuneCP5_Flat2017_13TeV_pythia8/crab_pickEvents/180201_223255/0000/pickevents_1.root' - # /afs/cern.ch/user/k/kirschen/public/forJERC/forMCTruthDebugging/pickevents_NoPU.root + 'root://cmsxrootd.fnal.gov//store/mc//.root', ) process.source = cms.Source("PoolSource", fileNames = inputFiles ) From ca0d9b0d02b520532b4f12170a1337f127e0d972 Mon Sep 17 00:00:00 2001 From: Sifu Luo Date: Sun, 4 Mar 2018 18:17:17 -0600 Subject: [PATCH 06/10] instance with nentries <= 4 will not be shown in closure --- JetUtilities/src/ClosureMaker.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/JetUtilities/src/ClosureMaker.cc b/JetUtilities/src/ClosureMaker.cc index f26db55c..3e09676c 100644 --- a/JetUtilities/src/ClosureMaker.cc +++ b/JetUtilities/src/ClosureMaker.cc @@ -386,10 +386,10 @@ void ClosureMaker::loopOverBins(TH2F* hvar, unsigned int iVarBin) { continue; } } - else if(h.back()->GetEntries()<=4 && h.back()->GetEntries()>1) { - hClosure.back()->SetBinContent(ibin+1,h.back()->GetMean()); - hClosure.back()->SetBinError(ibin+1,h.back()->GetMeanError()); - } + // else if(h.back()->GetEntries()<=4 && h.back()->GetEntries()>1) { + // hClosure.back()->SetBinContent(ibin+1,h.back()->GetMean()); + // hClosure.back()->SetBinError(ibin+1,h.back()->GetMeanError()); + // } else { continue; } From b7b393d3b56225a791d7ac1964414ac729bcdc74 Mon Sep 17 00:00:00 2001 From: Sifu Luo Date: Mon, 2 Jul 2018 15:58:40 -0500 Subject: [PATCH 07/10] updated --- JetAnalyzers/bin/compareJEC.C | 12 +--- JetUtilities/interface/ClosureMaker.hh | 77 +++++++++++++------------- JetUtilities/interface/L2Creator.hh | 4 +- JetUtilities/src/L2Creator.cc | 54 +++++++++++------- 4 files changed, 79 insertions(+), 68 deletions(-) diff --git a/JetAnalyzers/bin/compareJEC.C b/JetAnalyzers/bin/compareJEC.C index f6c6e478..37d0575e 100644 --- a/JetAnalyzers/bin/compareJEC.C +++ b/JetAnalyzers/bin/compareJEC.C @@ -255,15 +255,9 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa int maxTries = 7; string strPath; vector paths = {"","CondFormats/JetMETObjects/data/"}; - paths.push_back(string("/fdata/hepx/store/user/siluo/JEC/94X_Fall17/")+cid1+"/"); - paths.push_back(string("/fdata/hepx/store/user/siluo/JEC/80X_Summer16/")+cid2+"/"); - paths.push_back(string("/fdata/hepx/store/user/siluo/JEC/80X_stitch/")+cid3+"/"); - // paths.push_back(string("/fdata/hepx/store/user/aperloff/JEC/80X_Summer16/")+cid1+"/"); - // paths.push_back(string("/fdata/hepx/store/user/aperloff/JEC/80X_Summer16/")+cid2+"/"); - // paths.push_back(string("/fdata/hepx/store/user/aperloff/JEC/80X_Summer16/")+cid3+"/"); - paths.push_back(string("/home/sifuluo/jec/CMSSW_8_0_20/src/JetMETCorrections/JECDatabase/textFiles/")+cid1+"/"); - paths.push_back(string("/home/sifuluo/jec/CMSSW_8_0_20/src/JetMETCorrections/JECDatabase/textFiles/")+cid2+"/"); - paths.push_back(string("/home/sifuluo/jec/CMSSW_8_0_20/src/JetMETCorrections/JECDatabase/textFiles/")+cid3+"/"); + paths.push_back(string("")+cid1+"/"); + paths.push_back(string("")+cid2+"/"); + paths.push_back(string("")+cid3+"/"); // JEC1 diff --git a/JetUtilities/interface/ClosureMaker.hh b/JetUtilities/interface/ClosureMaker.hh index 3b809829..29a6501d 100644 --- a/JetUtilities/interface/ClosureMaker.hh +++ b/JetUtilities/interface/ClosureMaker.hh @@ -62,49 +62,48 @@ using namespace std; //////////////////////////////////////////////////////////////////////////////// class ClosureMaker { public: - ClosureMaker(); - ClosureMaker(CommandLine& cl); - - // - // member functions - // - void openInputFile(); - void getHistograms(TDirectoryFile* idir); - void openOutputFile(); - void closeFiles(); - void makeLines(); - void loopOverDirectories(); - void loopOverAlgorithms(); - void resetForNextAlgorithm(); - void loopOverHistograms(); - void etaClosureMergedPtBins(TDirectoryFile* idir); - void loopOverBins(TH2F* hvar, unsigned int iVarBin); - void adjust_fitrange(TH1* h,double& min,double& max); - void checkResponse(); - pair determineCanvasRange(double xmin, double xmax); - void makeCanvases(); - void makeMergedCanvas(bool finemerge); - void writeToFile(); - void makeClosure(const VARIABLES::Variable ivar = VARIABLES::refpt); + ClosureMaker(); + ClosureMaker(CommandLine& cl); + // + // member functions + // + void openInputFile(); + void getHistograms(TDirectoryFile* idir); + void openOutputFile(); + void closeFiles(); + void makeLines(); + void loopOverDirectories(); + void loopOverAlgorithms(); + void resetForNextAlgorithm(); + void loopOverHistograms(); + void etaClosureMergedPtBins(TDirectoryFile* idir); + void loopOverBins(TH2F* hvar, unsigned int iVarBin); + void adjust_fitrange(TH1* h,double& min,double& max); + void checkResponse(); + pair determineCanvasRange(double xmin, double xmax); + void makeCanvases(); + void makeMergedCanvas(bool finemerge); + void writeToFile(); + void makeClosure(const VARIABLES::Variable ivar = VARIABLES::refpt); private: - bool objects_loaded, draw_guidelines; - double CMEnergy, nsigma; - TString path, filename, outputDir, outputFilename, flavor, alg, histMet; - vector algs, outputFormat; - JetInfo *ji; - TFile *ifile, *ofile; - ObjectLoader hl; + bool objects_loaded, draw_guidelines; + double CMEnergy, nsigma; + TString path, filename, outputDir, outputFilename, flavor, alg, histMet; + vector algs, outputFormat; + JetInfo *ji; + TFile *ifile, *ofile; + ObjectLoader hl; vector h; vector func; - vector hClosure; - TF1 *line, *linePlus, *lineMinus; - vector > canvases_legends; - vector pave; - VARIABLES::Variable var; - TDirectoryFile *odir; - HistUtil::HistogramMetric histogramMetric; - int statTh; + vector hClosure; + TF1 *line, *linePlus, *lineMinus; + vector > canvases_legends; + vector pave; + VARIABLES::Variable var; + TDirectoryFile *odir; + HistUtil::HistogramMetric histogramMetric; + int statTh; }; #endif diff --git a/JetUtilities/interface/L2Creator.hh b/JetUtilities/interface/L2Creator.hh index 345f5171..b52bd92e 100644 --- a/JetUtilities/interface/L2Creator.hh +++ b/JetUtilities/interface/L2Creator.hh @@ -99,7 +99,7 @@ private: TString output, outputDir, l2calofit, l2pffit; vector formats, algs; bool l2l3, delphes; - int maxFitIter, statTh; + int maxFitIter, statThreshold; HistUtil::HistogramMetric histogramMetric; TFile* ofile; TFile* ifile; @@ -116,6 +116,8 @@ private: vector vabscor_eta; vector vrelcor_eta; vector vabscor_eta_spline; + vector ptclipcones; + vector ptclips; float ptclip; }; diff --git a/JetUtilities/src/L2Creator.cc b/JetUtilities/src/L2Creator.cc index 7becd59d..7e7e0b7d 100644 --- a/JetUtilities/src/L2Creator.cc +++ b/JetUtilities/src/L2Creator.cc @@ -31,23 +31,25 @@ L2Creator::L2Creator(CommandLine& cl) { // // evaluate command-line / configuration file options // - input = cl.getValue ("input"); - era = cl.getValue ("era"); - l3input = cl.getValue ("l3input", "l3.root"); - output = cl.getValue ("output", "l2.root"); - outputDir = cl.getValue ("outputDir", "./"); - formats = cl.getVector ("formats", ""); - algs = cl.getVector ("algs", ""); - l2l3 = cl.getValue ("l2l3", true); - l2calofit = cl.getValue ("l2calofit", "standard"); - l2pffit = cl.getValue ("l2pffit", "standard"); - delphes = cl.getValue ("delphes", false); - maxFitIter = cl.getValue ("maxFitIter", 30); - histMet = cl.getValue ("histMet", "mu_h"); + input = cl.getValue ("input"); + era = cl.getValue ("era"); + l3input = cl.getValue ("l3input", "l3.root"); + output = cl.getValue ("output", "l2.root"); + outputDir = cl.getValue ("outputDir", "./"); + formats = cl.getVector ("formats", ""); + algs = cl.getVector ("algs", ""); + l2l3 = cl.getValue ("l2l3", true); + l2calofit = cl.getValue ("l2calofit", "standard"); + l2pffit = cl.getValue ("l2pffit", "standard"); + delphes = cl.getValue ("delphes", false); + maxFitIter = cl.getValue ("maxFitIter", 30); + histMet = cl.getValue ("histMet", "mu_h"); histogramMetric = HistUtil::getHistogramMetricType(histMet); - ptclip = cl.getValue ("ptclip", 0.); - statTh = cl.getValue ("statTh", 4); + ptclipcones = cl.getVector ("ptclipcones", ""); + ptclips = cl.getVector ("ptclips", ""); + ptclip = cl.getValue ("ptclip", 0.); + statThreshold = cl.getValue ("statThreshold", 4); if (!cl.partialCheck()) return; cl.print(); @@ -142,6 +144,20 @@ void L2Creator::loopOverAlgorithms(string makeCanvasVariable) { // ji = new JetInfo(alg); + // Set ptclip + if (ptclipcones.size() != 0){ + if (ptclips.size() != ptclipcones.size()){ + ptclips.resize(ptclipcones.size(),ptclip); + cout << "ptclips and ptclipcones have different size, ptclips is filled with default ptclip = " << ptclip <getConeSize().Atoi() == ptclipcones[conesit]) { + ptclip = ptclips[conesit]; + cout << "for " << alg << " ptclip is set to " << ptclip <GetEntries() > statTh) {//hrsp->Integral()!=0) { + if (hrsp->GetEntries() > statThreshold) {//hrsp->Integral()!=0) { //TF1* frsp = (TF1*)hrsp->GetListOfFunctions()->Last(); //std::cout << "hrspName = " << hrsp->GetName() << ": frsp = " << frsp << std::endl; @@ -626,12 +642,12 @@ bool L2Creator::checkFormulaEvaluator() { // // Check that ipt is not outside [ptmin,ptmax] // - if (iptptmax) continue; + if (ipt < ptmin || ipt > ptmax) continue; // // Check that the ipt is not outside the pt clipping area // - if(ipt > pt_limit) continue; + if (ipt < ptclip || ipt > pt_limit) continue; // // Set the inputs for the FactorizedJetCorrector @@ -1381,7 +1397,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { <getNpar()+2) //Number of parameters + 2 - < bounds.first) ? pt_clip : bounds.first) <setParameters(isection); for(int p=0; pgetNpar(); p++) { From d7e884766905bfb1b1e98ba7b62ef966c46055ab Mon Sep 17 00:00:00 2001 From: Sifu Luo Date: Thu, 5 Jul 2018 16:57:04 -0500 Subject: [PATCH 08/10] minor update ready to merge. --- JetAnalyzers/bin/compareJEC.C | 10 +--------- JetAnalyzers/bin/jet_synchtest_x.cc | 2 -- 2 files changed, 1 insertion(+), 11 deletions(-) diff --git a/JetAnalyzers/bin/compareJEC.C b/JetAnalyzers/bin/compareJEC.C index 37d0575e..bb2863d6 100644 --- a/JetAnalyzers/bin/compareJEC.C +++ b/JetAnalyzers/bin/compareJEC.C @@ -225,7 +225,7 @@ double getEtaPtUncert(JetCorrectionUncertainty *unc, void compareJEC(string payld1="Winter14_V8", string payld2="", string payld3="", string algo1="AK5PFchs", string algo2="AK5PFchs", string algo3 ="AK5PFchs", string type1="DATA", string type2="DATA", string type3 ="DATA", - bool l1=true, bool l2l3=true, bool res=true) { + bool l1=true, bool l2l3=true, bool res=true) { //gROOT->ProcessLine(".L tdrstyle_mod14_ia.C"); setTDRStyle(); @@ -1102,12 +1102,4 @@ int main(int argc,char**argv) compareJEC("Fall17_25nsV1", "Summer16_25nsV5", "bias2SelectionPow_25nsV1", "AK4PFchs", "AK4PFchs", "AK4PFchs", "MC","MC","MC", false, true, false); compareJEC("Fall17_25nsV1", "Summer16_25nsV5", "bias2SelectionPow_25nsV1", "AK4PFchs", "AK4PFchs", "AK4PFchs", "MC","MC","MC", true, false, false); - // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFchs", "AK8PFchs", "AK8PFchs", "MC","MC","MC", false, true, false); - // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFchs", "AK8PFchs", "AK8PFchs", "MC","MC","MC", true, false, false); - // - // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK4PFPuppi", "AK4PFPuppi", "AK4PFPuppi", "MC","MC","MC", false, true, false); - // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK4PFPuppi", "AK4PFPuppi", "AK4PFPuppi", "MC","MC","MC", true, false, false); - // - // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFPuppi", "AK8PFPuppi", "AK8PFPuppi", "MC","MC","MC", false, true, false); - // compareJEC("Summer16_25nsV5", "Summer16_25nsV4", "Spring16_25nsV6", "AK8PFPuppi", "AK8PFPuppi", "AK8PFPuppi", "MC","MC","MC", true, false, false); } diff --git a/JetAnalyzers/bin/jet_synchtest_x.cc b/JetAnalyzers/bin/jet_synchtest_x.cc index 0dc7bf18..11bacc82 100644 --- a/JetAnalyzers/bin/jet_synchtest_x.cc +++ b/JetAnalyzers/bin/jet_synchtest_x.cc @@ -796,9 +796,7 @@ void MatchEventsAndJets::LoopOverEvents(bool verbose, bool reduceHistograms, str for (IT::const_iterator it = mapTreePU.begin(); it != mapTreePU.end(); ++it) { if (iftest && nevs >= maxEvts) return; - // if (nevs >= 10000000) return; - //if (nevs%10000==0) cout << "\t"< Date: Wed, 8 Aug 2018 03:47:28 -0500 Subject: [PATCH 09/10] compareJEC now works with commandline configurations and can generate dummyl3 files upon request. --- JetAnalyzers/bin/compareJEC.C | 122 +++++++++++++++++++--------- JetUtilities/interface/L2Creator.hh | 5 +- JetUtilities/src/L2Creator.cc | 48 ++++++----- 3 files changed, 108 insertions(+), 67 deletions(-) diff --git a/JetAnalyzers/bin/compareJEC.C b/JetAnalyzers/bin/compareJEC.C index bb2863d6..1bc8344f 100644 --- a/JetAnalyzers/bin/compareJEC.C +++ b/JetAnalyzers/bin/compareJEC.C @@ -22,7 +22,10 @@ #include #include #include +#include +#include +#include "JetMETAnalysis/JetUtilities/interface/CommandLine.h" #include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h" #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" @@ -30,7 +33,7 @@ //#include "CondFormats/JetMETObjects/interface/SimpleJetCorrector.h" //#include "CondFormats/JetMETObjects/interface/SimpleJetCorrectionUncertainty.h" -#include "FWCore/ParameterSet/interface/FileInPath.h" +// #include "FWCore/ParameterSet/interface/FileInPath.h" #include "JetMETAnalysis/JetUtilities/interface/Style.h" @@ -219,10 +222,16 @@ double getEtaPtUncert(JetCorrectionUncertainty *unc, return (unc->getUncertainty(true)); } // getEtaPtUncert +// Test if file of given path exists +bool fileexist (const std::string& name) { + struct stat buffer; + return (stat (name.c_str(), &buffer) == 0); +} // ------------------------------------------------------------ // -void compareJEC(string payld1="Winter14_V8", string payld2="", string payld3="", +void compareJEC(string path1 ="", string path2 ="", string path3 ="", + string payld1="Winter14_V8", string payld2="", string payld3="", string algo1="AK5PFchs", string algo2="AK5PFchs", string algo3 ="AK5PFchs", string type1="DATA", string type2="DATA", string type3 ="DATA", bool l1=true, bool l2l3=true, bool res=true) { @@ -252,12 +261,12 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa const char *cid3 = sid3.c_str(); const char *a3 = algo3.c_str(); - int maxTries = 7; + // int maxTries = 7; string strPath; - vector paths = {"","CondFormats/JetMETObjects/data/"}; - paths.push_back(string("")+cid1+"/"); - paths.push_back(string("")+cid2+"/"); - paths.push_back(string("")+cid3+"/"); + vector paths; + paths.push_back(path1); + paths.push_back(path2); + paths.push_back(path3); // JEC1 @@ -279,16 +288,17 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa //strVec.push_back(Form("%s_Uncertainty_%s.txt",cid1,a1)); for(unsigned int istr=0;istr ("basepath", "/"); + string path1 = cl.getValue ("path1", ""); + string path2 = cl.getValue ("path2", ""); + string path3 = cl.getValue ("path3", ""); + vector eras = cl.getVector ("eras", "Fall17_25nsV1"); + vector algos = cl.getVector ("algos", "AK4PF"); + vector types = cl.getVector ("types", ""); + bool dummyl3 = cl.getValue ("dummyl3", true); + + if(!cl.check()) return 0; + cl.print(); + if (path1.substr(0,1) != "/") path1 = basepath + path1; + if (path2.substr(0,1) != "/") path2 = basepath + path2; + if (path3.substr(0,1) != "/") path3 = basepath + path3; + cout << "paths are: "<< endl << path1 < paths = {path1,path2,path3}; + if (algos.size() != 3) algos.resize(3,algos[0]); + if (eras.size() != 3) eras.resize(3,eras[0]); + if (types.size() != 3) types.resize(3,"MC"); + + if (dummyl3){ + for (unsigned ip = 0; ip < paths.size(); ++ip){ + ofstream l3f(Form("%s%s_%s_L3Absolute_%s.txt", paths[ip].c_str(), eras[ip].c_str(),types[ip].c_str(), algos[ip].c_str())); + l3f << "{1 JetEta 1 JetPt 1 Correction L3Absolute}"<& collection,const string& element); @@ -92,7 +92,7 @@ public: void setAndFitFLogAndFGaus(TGraphErrors* gabscor, TF1* flog, TF1* fgaus, double xmin); void setOfflinePFParameters(TGraphErrors* gabscor, TF1* fabscor, double xmin, double xmax); void writeTextFileForCurrentAlgorithm(); - void writeTextFileForCurrentAlgorithm_spline(); + void writeTextFileForCurrentAlgorithm_spline(float ptclip); private: string input, era, l3input, histMet; @@ -118,7 +118,6 @@ private: vector vabscor_eta_spline; vector ptclipcones; vector ptclips; - float ptclip; }; #endif diff --git a/JetUtilities/src/L2Creator.cc b/JetUtilities/src/L2Creator.cc index 7e7e0b7d..62c762db 100644 --- a/JetUtilities/src/L2Creator.cc +++ b/JetUtilities/src/L2Creator.cc @@ -47,8 +47,7 @@ L2Creator::L2Creator(CommandLine& cl) { histogramMetric = HistUtil::getHistogramMetricType(histMet); ptclipcones = cl.getVector ("ptclipcones", ""); - ptclips = cl.getVector ("ptclips", ""); - ptclip = cl.getValue ("ptclip", 0.); + ptclips = cl.getVector ("ptclips", "8."); statThreshold = cl.getValue ("statThreshold", 4); if (!cl.partialCheck()) return; @@ -143,21 +142,6 @@ void L2Creator::loopOverAlgorithms(string makeCanvasVariable) { // Make the JetInfo object // ji = new JetInfo(alg); - - // Set ptclip - if (ptclipcones.size() != 0){ - if (ptclips.size() != ptclipcones.size()){ - ptclips.resize(ptclipcones.size(),ptclip); - cout << "ptclips and ptclipcones have different size, ptclips is filled with default ptclip = " << ptclip <getConeSize().Atoi() == ptclipcones[conesit]) { - ptclip = ptclips[conesit]; - cout << "for " << alg << " ptclip is set to " << ptclip <getConeSize().Atoi() == ptclipcones[conesit]) { + ptclip = ptclips[conesit]; + break; + } + } + } + cout << "for " << alg << " ptclip is set to " << ptclip < pt_to_check = {10.0,30.0,100.0,500.0,1000.0,2000.0,3000.0,4000.0}; unsigned int vector_size = 0; @@ -1324,7 +1323,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm() { } //______________________________________________________________________________ -void L2Creator::writeTextFileForCurrentAlgorithm_spline() { +void L2Creator::writeTextFileForCurrentAlgorithm_spline(float ptclip) { TString txtfilename = outputDir + era + "_L2Relative_" + ji->getAlias() + ".txt"; ofstream fout(txtfilename); fout.setf(ios::right); @@ -1332,7 +1331,6 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //For eta-dependent spline clipping int pt_limit = 70; - float pt_clip = ptclip; unsigned int vector_size = 0; vector_size = vabscor_eta.size(); @@ -1385,7 +1383,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //if(isection==spline->getNSections()-1) { // bounds.second = 6500; //} - if(bounds.second < pt_clip) continue; + if(bounds.second < ptclip) continue; if(isection==spline->getNSections()-1) lastLine = true; if(bounds.second >= pt_limit) { abovePtLimit = true; @@ -1397,7 +1395,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { <getNpar()+2) //Number of parameters + 2 - < bounds.first) ? pt_clip : bounds.first) + < bounds.first) ? ptclip : bounds.first) <setParameters(isection); for(int p=0; pgetNpar(); p++) { From b408bdfb347c92320d6fa8c0a49aa767c098bfc9 Mon Sep 17 00:00:00 2001 From: Sifu Luo Date: Wed, 8 Aug 2018 03:47:28 -0500 Subject: [PATCH 10/10] compareJEC now works with commandline configurations and can generate dummyl3 files upon request. --- JetAnalyzers/bin/compareJEC.C | 122 +++++++++++++++++++--------- JetUtilities/interface/L2Creator.hh | 6 +- JetUtilities/src/L2Creator.cc | 51 ++++++------ 3 files changed, 111 insertions(+), 68 deletions(-) diff --git a/JetAnalyzers/bin/compareJEC.C b/JetAnalyzers/bin/compareJEC.C index bb2863d6..1bc8344f 100644 --- a/JetAnalyzers/bin/compareJEC.C +++ b/JetAnalyzers/bin/compareJEC.C @@ -22,7 +22,10 @@ #include #include #include +#include +#include +#include "JetMETAnalysis/JetUtilities/interface/CommandLine.h" #include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h" #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" @@ -30,7 +33,7 @@ //#include "CondFormats/JetMETObjects/interface/SimpleJetCorrector.h" //#include "CondFormats/JetMETObjects/interface/SimpleJetCorrectionUncertainty.h" -#include "FWCore/ParameterSet/interface/FileInPath.h" +// #include "FWCore/ParameterSet/interface/FileInPath.h" #include "JetMETAnalysis/JetUtilities/interface/Style.h" @@ -219,10 +222,16 @@ double getEtaPtUncert(JetCorrectionUncertainty *unc, return (unc->getUncertainty(true)); } // getEtaPtUncert +// Test if file of given path exists +bool fileexist (const std::string& name) { + struct stat buffer; + return (stat (name.c_str(), &buffer) == 0); +} // ------------------------------------------------------------ // -void compareJEC(string payld1="Winter14_V8", string payld2="", string payld3="", +void compareJEC(string path1 ="", string path2 ="", string path3 ="", + string payld1="Winter14_V8", string payld2="", string payld3="", string algo1="AK5PFchs", string algo2="AK5PFchs", string algo3 ="AK5PFchs", string type1="DATA", string type2="DATA", string type3 ="DATA", bool l1=true, bool l2l3=true, bool res=true) { @@ -252,12 +261,12 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa const char *cid3 = sid3.c_str(); const char *a3 = algo3.c_str(); - int maxTries = 7; + // int maxTries = 7; string strPath; - vector paths = {"","CondFormats/JetMETObjects/data/"}; - paths.push_back(string("")+cid1+"/"); - paths.push_back(string("")+cid2+"/"); - paths.push_back(string("")+cid3+"/"); + vector paths; + paths.push_back(path1); + paths.push_back(path2); + paths.push_back(path3); // JEC1 @@ -279,16 +288,17 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa //strVec.push_back(Form("%s_Uncertainty_%s.txt",cid1,a1)); for(unsigned int istr=0;istr ("basepath", "/"); + string path1 = cl.getValue ("path1", ""); + string path2 = cl.getValue ("path2", ""); + string path3 = cl.getValue ("path3", ""); + vector eras = cl.getVector ("eras", "Fall17_25nsV1"); + vector algos = cl.getVector ("algos", "AK4PF"); + vector types = cl.getVector ("types", ""); + bool dummyl3 = cl.getValue ("dummyl3", true); + + if(!cl.check()) return 0; + cl.print(); + if (path1.substr(0,1) != "/") path1 = basepath + path1; + if (path2.substr(0,1) != "/") path2 = basepath + path2; + if (path3.substr(0,1) != "/") path3 = basepath + path3; + cout << "paths are: "<< endl << path1 < paths = {path1,path2,path3}; + if (algos.size() != 3) algos.resize(3,algos[0]); + if (eras.size() != 3) eras.resize(3,eras[0]); + if (types.size() != 3) types.resize(3,"MC"); + + if (dummyl3){ + for (unsigned ip = 0; ip < paths.size(); ++ip){ + ofstream l3f(Form("%s%s_%s_L3Absolute_%s.txt", paths[ip].c_str(), eras[ip].c_str(),types[ip].c_str(), algos[ip].c_str())); + l3f << "{1 JetEta 1 JetPt 1 Correction L3Absolute}"<& collection,const string& element); @@ -92,7 +92,7 @@ public: void setAndFitFLogAndFGaus(TGraphErrors* gabscor, TF1* flog, TF1* fgaus, double xmin); void setOfflinePFParameters(TGraphErrors* gabscor, TF1* fabscor, double xmin, double xmax); void writeTextFileForCurrentAlgorithm(); - void writeTextFileForCurrentAlgorithm_spline(); + void writeTextFileForCurrentAlgorithm_spline(float ptclip); private: string input, era, l3input, histMet; @@ -116,9 +116,9 @@ private: vector vabscor_eta; vector vrelcor_eta; vector vabscor_eta_spline; + float ptclip; vector ptclipcones; vector ptclips; - float ptclip; }; #endif diff --git a/JetUtilities/src/L2Creator.cc b/JetUtilities/src/L2Creator.cc index 7e7e0b7d..edf61e89 100644 --- a/JetUtilities/src/L2Creator.cc +++ b/JetUtilities/src/L2Creator.cc @@ -46,9 +46,9 @@ L2Creator::L2Creator(CommandLine& cl) { histMet = cl.getValue ("histMet", "mu_h"); histogramMetric = HistUtil::getHistogramMetricType(histMet); - ptclipcones = cl.getVector ("ptclipcones", ""); - ptclips = cl.getVector ("ptclips", ""); ptclip = cl.getValue ("ptclip", 0.); + ptclipcones = cl.getVector ("ptclipcones" ""); + ptclips = cl.getVector ("ptclips" ""); statThreshold = cl.getValue ("statThreshold", 4); if (!cl.partialCheck()) return; @@ -143,21 +143,6 @@ void L2Creator::loopOverAlgorithms(string makeCanvasVariable) { // Make the JetInfo object // ji = new JetInfo(alg); - - // Set ptclip - if (ptclipcones.size() != 0){ - if (ptclips.size() != ptclipcones.size()){ - ptclips.resize(ptclipcones.size(),ptclip); - cout << "ptclips and ptclipcones have different size, ptclips is filled with default ptclip = " << ptclip <getConeSize().Atoi() == ptclipcones[conesit]) { - ptclip = ptclips[conesit]; - cout << "for " << alg << " ptclip is set to " << ptclip <getConeSize().Atoi() == ptclipcones[conesit]) { + ptcliping = ptclips[conesit]; + break; + } + } + } + cout << "for " << alg << " ptcliping is set to " << ptcliping < pt_to_check = {10.0,30.0,100.0,500.0,1000.0,2000.0,3000.0,4000.0}; unsigned int vector_size = 0; @@ -647,7 +647,7 @@ bool L2Creator::checkFormulaEvaluator() { // // Check that the ipt is not outside the pt clipping area // - if (ipt < ptclip || ipt > pt_limit) continue; + if (ipt < ptcliping || ipt > pt_limit) continue; // // Set the inputs for the FactorizedJetCorrector @@ -1324,7 +1324,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm() { } //______________________________________________________________________________ -void L2Creator::writeTextFileForCurrentAlgorithm_spline() { +void L2Creator::writeTextFileForCurrentAlgorithm_spline(float ptcliping) { TString txtfilename = outputDir + era + "_L2Relative_" + ji->getAlias() + ".txt"; ofstream fout(txtfilename); fout.setf(ios::right); @@ -1332,7 +1332,6 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //For eta-dependent spline clipping int pt_limit = 70; - float pt_clip = ptclip; unsigned int vector_size = 0; vector_size = vabscor_eta.size(); @@ -1385,7 +1384,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //if(isection==spline->getNSections()-1) { // bounds.second = 6500; //} - if(bounds.second < pt_clip) continue; + if(bounds.second < ptcliping) continue; if(isection==spline->getNSections()-1) lastLine = true; if(bounds.second >= pt_limit) { abovePtLimit = true; @@ -1397,7 +1396,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { <getNpar()+2) //Number of parameters + 2 - < bounds.first) ? pt_clip : bounds.first) + < bounds.first) ? ptcliping : bounds.first) <setParameters(isection); for(int p=0; pgetNpar(); p++) {