diff --git a/JetAnalyzers/bin/compareJEC.C b/JetAnalyzers/bin/compareJEC.C index c8cbacf4..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" @@ -67,10 +70,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 +91,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 +123,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 +146,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 +166,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 +178,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 +206,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 +214,7 @@ double getEtaPtUncert(JetCorrectionUncertainty *unc, _thejec = jec; fCorrPt->SetParameters(eta, mu); double ptgen = fCorrPt->Eval(ptreco); - + unc->setJetEta(eta); unc->setJetPt(ptgen); } @@ -219,13 +222,19 @@ 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) { + bool l1=true, bool l2l3=true, bool res=true) { //gROOT->ProcessLine(".L tdrstyle_mod14_ia.C"); setTDRStyle(); @@ -252,19 +261,17 @@ 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("/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; + paths.push_back(path1); + paths.push_back(path2); + paths.push_back(path3); + // JEC1 bool noL2L3_p1 = TString(payld1).Contains("Run1Sub") || TString(payld1).Contains("nsRunIISpring15") || TString(payld1).Contains("RC") ; - // + // // noL2L3_p1 = true ; //ctmp JetCorrectionUncertainty *jecUnc1 = 0; @@ -281,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 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 +449,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 +470,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 +479,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 +493,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 +534,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 +545,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 +564,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 +575,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 +635,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 +665,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 +725,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 +744,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 +800,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 +829,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 +856,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 +889,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 +907,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 +927,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 +939,7 @@ void compareJEC(string payld1="Winter14_V8", string payld2="", string pa ///// leg1a->AddEntry(g2a,"2016","LPF"); } - + gPad->RedrawAxis(); @@ -944,12 +954,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 +967,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 +986,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 +1029,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 +1043,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 +1055,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 +1064,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 +1088,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 +1111,39 @@ 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); + CommandLine cl; + if (!cl.parse(argc,argv)) return 0; + string basepath = cl.getValue ("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}"<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, 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); @@ -796,7 +797,6 @@ void MatchEventsAndJets::LoopOverEvents(bool verbose, bool reduceHistograms, str if (iftest && nevs >= maxEvts) 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 +1186,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 +1356,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 +1372,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 +1389,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 - - diff --git a/JetUtilities/interface/ClosureMaker.hh b/JetUtilities/interface/ClosureMaker.hh index 96dbcaa0..29a6501d 100644 --- a/JetUtilities/interface/ClosureMaker.hh +++ b/JetUtilities/interface/ClosureMaker.hh @@ -62,48 +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(); - 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; - vector h; - vector func; - vector hClosure; - TF1 *line, *linePlus, *lineMinus; - vector > canvases_legends; - vector pave; - VARIABLES::Variable var; - TDirectoryFile *odir; - HistUtil::HistogramMetric histogramMetric; + 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; }; #endif diff --git a/JetUtilities/interface/HistogramUtilities.hh b/JetUtilities/interface/HistogramUtilities.hh index 150cb666..b8802eb9 100644 --- a/JetUtilities/interface/HistogramUtilities.hh +++ b/JetUtilities/interface/HistogramUtilities.hh @@ -15,6 +15,7 @@ #include "TMath.h" #include "TH2.h" #include "TF1.h" +#include "TSpectrum.h" //C++ libraries #include @@ -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/interface/L2Creator.hh b/JetUtilities/interface/L2Creator.hh index 35785efc..4985f93b 100644 --- a/JetUtilities/interface/L2Creator.hh +++ b/JetUtilities/interface/L2Creator.hh @@ -68,7 +68,7 @@ public: L2Creator(); L2Creator(CommandLine& cl); - bool checkFormulaEvaluator(); + bool checkFormulaEvaluator(float ptclip); void closeFiles(); /// check if a vector of strings contains a certain element bool contains(const vector& collection,const string& element); @@ -92,14 +92,14 @@ 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; TString output, outputDir, l2calofit, l2pffit; vector formats, algs; bool l2l3, delphes; - int maxFitIter; + int maxFitIter, statThreshold; HistUtil::HistogramMetric histogramMetric; TFile* ofile; TFile* ifile; @@ -116,6 +116,9 @@ private: vector vabscor_eta; vector vrelcor_eta; vector vabscor_eta_spline; + float ptclip; + vector ptclipcones; + vector ptclips; }; #endif diff --git a/JetUtilities/src/ClosureMaker.cc b/JetUtilities/src/ClosureMaker.cc index 77dc7878..3e09676c 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; @@ -88,7 +89,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 +103,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 +126,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 +210,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 +259,7 @@ void ClosureMaker::loopOverAlgorithms() { // makeCanvases(); if(var == VARIABLES::refpt || var == VARIABLES::jtpt || var == VARIABLES::ptclpt) { - makeMergedCanvas(); + makeMergedCanvas(false); } writeToFile(); @@ -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); @@ -385,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; } @@ -421,7 +422,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 +492,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 +534,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 +549,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 +637,3 @@ void ClosureMaker::makeClosure(const VARIABLES::Variable ivar) { loopOverAlgorithms(); closeFiles(); } - - 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/L2Creator.cc b/JetUtilities/src/L2Creator.cc index de67892a..edf61e89 100644 --- a/JetUtilities/src/L2Creator.cc +++ b/JetUtilities/src/L2Creator.cc @@ -31,21 +31,26 @@ 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.); + ptclipcones = cl.getVector ("ptclipcones" ""); + ptclips = cl.getVector ("ptclips" ""); + statThreshold = cl.getValue ("statThreshold", 4); + if (!cl.partialCheck()) return; cl.print(); } @@ -138,7 +143,6 @@ void L2Creator::loopOverAlgorithms(string makeCanvasVariable) { // Make the JetInfo object // ji = new JetInfo(alg); - // // Get the response from the l3 file only if l2l3 is set to false; // @@ -177,13 +181,28 @@ void L2Creator::loopOverAlgorithms(string makeCanvasVariable) { // doRelCorFits(); } - + // Set ptcliping + float ptcliping = ptclip; + if (ptclipcones.size() != 0 ){ + if (ptclips.size() != ptclipcones.size() ){ + if (ptclips.size() < ptclipcones.size()) cout << "ptclips have more less entries than ptclipcones, missing conesize now using default ptcliping = " << ptclip <getConeSize().Atoi() == ptclipcones[conesit]) { + ptcliping = ptclips[conesit]; + break; + } + } + } + cout << "for " << alg << " ptcliping is set to " << ptcliping <GetEntries() > 4) {//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; @@ -556,7 +575,7 @@ void L2Creator::loopOverEtaBins() { } //______________________________________________________________________________ -bool L2Creator::checkFormulaEvaluator() { +bool L2Creator::checkFormulaEvaluator(float ptcliping) { vector pt_to_check = {10.0,30.0,100.0,500.0,1000.0,2000.0,3000.0,4000.0}; unsigned int vector_size = 0; @@ -623,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 < ptcliping || ipt > pt_limit) continue; // // Set the inputs for the FactorizedJetCorrector @@ -1305,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); @@ -1345,6 +1364,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 +1384,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //if(isection==spline->getNSections()-1) { // bounds.second = 6500; //} - + if(bounds.second < ptcliping) continue; if(isection==spline->getNSections()-1) lastLine = true; if(bounds.second >= pt_limit) { abovePtLimit = true; @@ -1373,16 +1393,17 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //For expediency of Summer16_25nsV5_MC do eta-dependent clipping fout<getNpar()+2) //Number of parameters + 2 - < bounds.first) ? ptcliping : bounds.first) <setParameters(isection); for(int p=0; pgetNpar(); p++) { fout<GetParameter(p); } fout< #include "JetMETAnalysis/JetUtilities/interface/Style.h" +#include "JetMETAnalysis/JetUtilities/interface/HistogramUtilities.hh" using namespace std; @@ -102,7 +103,7 @@ TCanvas * getCanvasFromFittingProcedure(TString cname , TProfile2D * prof, TStri int S_comp_ = cname.CompareTo("ParametersVsGenSumPtOA"); //cout <Fit("pol2","0QS","",50,1800); } @@ -114,7 +115,7 @@ TCanvas * getCanvasFromFittingProcedure(TString cname , TProfile2D * prof, TStri { fr = hoff->Fit("pol2","0QS","",5,35); } - + // Skip if fit failed if (fr->Status()){ @@ -126,7 +127,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 +148,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 +165,7 @@ TCanvas * getCanvasFromFittingProcedure(TString cname , TProfile2D * prof, TStri } return c; - + }//getCanvasFromFittingProcedure @@ -189,12 +190,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 +240,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 +290,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 +346,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,18 +354,19 @@ 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) { - 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 (valReset(); //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 +415,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 +504,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 +718,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 +750,7 @@ TCanvas * getResolutionNumDenom(TString cname, TString ctitle, TString algo, TH2 c->cd(2); histod->Draw("E"); c->cd(); - + return c; }//getResolutionNumDenom @@ -764,7 +765,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 +972,8 @@ TCanvas * getGausMeanOffsetOverPtref(TString cname, TString ctitle, TString algo cout<<"\t Doing fits for Mean "<SetLogx(); @@ -1060,8 +1061,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 +1086,7 @@ TCanvas * getGausMeanOffsetScale(TString cname, TString ctitle, TString algo, ve histo->Reset(); //histo->Clear(); histo->Sumw2(); - + double maxy = 0; double miny = 0; @@ -1094,14 +1095,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 +1163,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 +1451,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 +1607,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 +1636,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 +1778,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 +1793,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 +1806,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