diff --git a/JetAnalyzers/BuildFile.xml b/JetAnalyzers/BuildFile.xml index a4d9509b..47f1e2ab 100644 --- a/JetAnalyzers/BuildFile.xml +++ b/JetAnalyzers/BuildFile.xml @@ -23,9 +23,9 @@ - ---> diff --git a/JetAnalyzers/bin/jet_response_fitter_x.cc b/JetAnalyzers/bin/jet_response_fitter_x.cc index b332f6e3..8f3adaa6 100644 --- a/JetAnalyzers/bin/jet_response_fitter_x.cc +++ b/JetAnalyzers/bin/jet_response_fitter_x.cc @@ -10,6 +10,7 @@ #include "JetMETAnalysis/JetUtilities/interface/CommandLine.h" #include "JetMETAnalysis/JetUtilities/interface/JetInfo.hh" +#include "JetMETAnalysis/JetAnalyzers/interface/fitCrystalBall.h" #include #include @@ -17,26 +18,27 @@ #include #include #include +#include #include #include #include +#include +#include #include #include #include #include #include -#include +#include +#include #include using namespace std; - //////////////////////////////////////////////////////////////////////////////// // define local functions //////////////////////////////////////////////////////////////////////////////// - - /// default fit with gaussian in niter iteration of mean void fit_gaussian(TH1F*& hrsp, const double nsigma, @@ -84,29 +86,35 @@ int main(int argc,char**argv) CommandLine cl; if (!cl.parse(argc,argv)) return 0; - string input = cl.getValue ("input"); - string output = cl.getValue ("output", ""); - double nsigma = cl.getValue ("nsigma", 1.5); - float jtptmin = cl.getValue ("jtptmin", 1.0); - int niter = cl.getValue ("niter", 3); - int ndfmin = cl.getValue ("ndfmin", 5); - vector algs = cl.getVector("algs", ""); - int verbose = cl.getValue ("verbose", 0); - int fittype = cl.getValue ("fittype", 0); - bool doAbsRsp= cl.getValue ("doAbsRsp", true); - bool doRelRsp= cl.getValue ("doRelRsp", true); - bool doEtaRsp= cl.getValue ("doEtaRsp", false); - bool doPhiRsp= cl.getValue ("doPhiRsp", false); - bool doFlavor= cl.getValue ("doFlavor", false); - + string input = cl.getValue ("input"); + string output = cl.getValue ("output", ""); + double nsigma = cl.getValue ("nsigma", 1.5); + float jtptmin = cl.getValue ("jtptmin", 1.0); + int niter = cl.getValue ("niter", 3); + int ndfmin = cl.getValue ("ndfmin", 5); + vector algs = cl.getVector("algs", ""); + int verbose = cl.getValue ("verbose", 0); + int fittype = cl.getValue ("fittype", 0); + bool doAbsRsp = cl.getValue ("doAbsRsp", false); + bool doRelRsp = cl.getValue ("doRelRsp", true); + bool doEtaRsp = cl.getValue ("doEtaRsp", false); + bool doPhiRsp = cl.getValue ("doPhiRsp", false); + bool doFlavor = cl.getValue ("doFlavor", false); + string histName = cl.getValue ("histName" ""); + int polDeg = cl.getValue ("polDeg", 7); + bool normalized = cl.getValue ("normalized", false); // normalized histogram to 1 + string fitDir = cl.getValue ("fitDir", ""); // save CB fit directories + if (!cl.check()) return 0; cl.print(); - if (fittype<0 || fittype>1) { - cout<<"ERROR: fittype not known, choose 0 for GAUSS, 1 for DSCB";return 0; + if (fittype==0) cout<<"*** Fitting with distributions w/ GAUSS"<mkdir(idir->GetName()); odir->cd(); - + cout<GetListOfKeys()); TKey* histKey(0); while ((histKey=(TKey*)nextHist())) { if (strcmp(histKey->GetClassName(),"TH1F")!=0) continue; - + TH1F* hrsp = (TH1F*)histKey->ReadObj(); string histname(hrsp->GetName()); @@ -193,46 +198,56 @@ int main(int argc,char**argv) continue; } - //if(verbose>0) cout << "Attempting to fit " << histname << " ... " << endl; + // + // write all the histogram used later for tau calibration study + // + if( fittype==2 && (histname.find("RelRsp_JetEta")==string::npos || histname.find("_RefPt")==string::npos) ) { + hrsp->Write(); + continue; + } + //if(verbose>0) cout << "Attempting to fit " << histname << " ... " << endl; + // //if (histname.find("RelRsp")>5&&histname.find("AbsRsp")!=0) { // hrsp->Write(); // continue; //} - - if (hrsp->Integral()>0.0) { - int fitstatus(0); - if (0==fittype) fit_gaussian(hrsp,nsigma,jtptmin,niter,verbose); - else fitstatus = fit_dscb(hrsp,nsigma,jtptmin,niter,alg,verbose); - - TF1* fitfnc = (TF1*) hrsp->GetListOfFunctions()->Last(); - if (0!=fitfnc && 0==fitstatus) fitfnc->ResetBit(TF1::kNotDraw); - - if (verbose>0&&0!=fitfnc) - cout<<"histo: "<GetName()<<"-> fnc: "<GetName()<GetNDF()0) cout<<"NDOF("<GetName()<<")=" - <GetNDF() - <<" FOR "<GetName()<GetListOfFunctions()->Delete(); - } + + // + // start fit process + // + if (hrsp->Integral()>0) { + int fitstatus(0); + if (fittype==0) fit_gaussian(hrsp,nsigma,jtptmin,niter,verbose); + else if (fittype==1) fitstatus = fit_dscb(hrsp,nsigma,jtptmin,niter,alg,verbose); + else fitCrystalBall(hrsp, alg, histName, polDeg, normalized, fitDir); + + TF1* fitfnc = (TF1*) hrsp->GetListOfFunctions()->Last(); + + if (fitfnc!=0 && fitstatus==0) fitfnc->ResetBit(TF1::kNotDraw); + + if (verbose>0 && fitfnc!=0) cout<<"histo: "<GetName()<<"-> fnc: "<GetName()<GetNDF()0) cout<<"NDOF("<GetName()<<")=" <GetNDF() <<" FOR "<GetName()<GetListOfFunctions()->Delete(); + } } else { - if (verbose>0) - cout<<"NOT ENOUGH ENTRIES FOR "<GetName()<0) cout<<"NOT ENOUGH ENTRIES FOR "<GetName()<Write(); - } + + if(fittype!=2) hrsp->Write(); + } // while nextHist loop cout<<"response fits for *"+alg+"* completed ..."<Write(); odir->DeleteAll(); delete odir; cout<<" and saved!\n"<Close(); delete ifile; cout<<" DONE."< vv; //Added by Alexx double fitrange_min(0.0); if (alg.find("pf")!=string::npos) fitrange_min = std::max(rspMax,0.3); else if (alg.find("PF")!=string::npos) fitrange_min = std::max(rspMax,0.3); else { - double first_nonzero_bin = 0.0; - for(int ibin=1; ibin<=hrsp->GetNbinsX();ibin++) { - if(first_nonzero_bin==0 && hrsp->GetBinContent(ibin)!=0) { - first_nonzero_bin = hrsp->GetBinCenter(ibin); - break; - } - } - vv.push_back(rspMax); - vv.push_back(hrsp->GetXaxis()->GetXmin()); - vv.push_back(0.2); - vv.push_back(first_nonzero_bin); - fitrange_min = *std::max_element(vv.begin(),vv.end()); + double first_nonzero_bin = 0.0; + for(int ibin=1; ibin<=hrsp->GetNbinsX();ibin++) { + if(first_nonzero_bin==0 && hrsp->GetBinContent(ibin)!=0) { + first_nonzero_bin = hrsp->GetBinCenter(ibin); + break; + } + } + vv.push_back(rspMax); + vv.push_back(hrsp->GetXaxis()->GetXmin()); + vv.push_back(0.2); + vv.push_back(first_nonzero_bin); + fitrange_min = *std::max_element(vv.begin(),vv.end()); } //else fitrange_min = std::max(rspMax,0.2); double fitrange_max = 1.7; @@ -321,7 +335,7 @@ int fit_dscb(TH1F*& hrsp, adjust_fitrange(hrsp,fitrange_min,fitrange_max); //guesstimate_fitrange(hrsp,fitrange_min,fitrange_max,alg); - + TF1* fdscb = new TF1("fdscb",fnc_dscb,fitrange_min,fitrange_max,7); // set the std values @@ -357,22 +371,22 @@ int fit_dscb(TH1F*& hrsp, fdscb->SetParLimits(6,0.,100.);//25.); if(hrsp->GetEntries()==0) { - if(verbose>0) - cout << "Warning in : Fit data is empty" << endl - << "hrsp->GetName(): " << hrsp->GetName() << endl; - return -1; + if(verbose>0) + cout << "Warning in : Fit data is empty" << endl + << "hrsp->GetName(): " << hrsp->GetName() << endl; + return -1; } else { - fitstatus = hrsp->Fit(fdscb,"RQB0+"); + fitstatus = hrsp->Fit(fdscb,"RQB0+"); } - + if (0==fitstatus) i=999; delete fdscb; fdscb = hrsp->GetFunction("fdscb"); if (0==fdscb) return -1; - + norm = fdscb->GetParameter(0); aone = fdscb->GetParameter(3); pone = fdscb->GetParameter(4); @@ -390,7 +404,7 @@ int fit_dscb(TH1F*& hrsp, if (0!=fitstatus){ cout<<"fit_fdscb() to "<GetName() - <<" failed. Fitstatus: "<GetFunction("fdscb")->Delete(); } @@ -410,7 +424,7 @@ double fnc_dscb(double*xx,double*pp) double p1 = pp[4]; double a2 = pp[5]; double p2 = pp[6]; - + double u = (x-mu)/sig; double A1 = TMath::Power(p1/TMath::Abs(a1),p1)*TMath::Exp(-a1*a1/2); double A2 = TMath::Power(p2/TMath::Abs(a2),p2)*TMath::Exp(-a2*a2/2); @@ -424,18 +438,17 @@ double fnc_dscb(double*xx,double*pp) return result; } - //______________________________________________________________________________ void fit_gaussian(TH1F*& hrsp, - const double nsigma, - const double jtptmin, - const int niter, - const int verbose) + 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(); @@ -484,31 +497,31 @@ void fit_gaussian(TH1F*& hrsp, fitfnc->SetParameter(1,peak); fitfnc->SetParameter(2,sigma); if(hrsp->GetEntries()==0) { - if(verbose>0) - cout << "Warning in : Fit data is empty" << endl - << "hrsp->GetName(): " << hrsp->GetName() << endl; - return; + if(verbose>0) + cout << "Warning in : Fit data is empty" << endl + << "hrsp->GetName(): " << hrsp->GetName() << endl; + return; } else { - fitstatus = hrsp->Fit(fitfnc,"RQ0"); + 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); + 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; - } + { + cout << "No function recorded in histogram " << hrsp->GetName() << endl; + } if (0!=fitstatus){ cout<<"fit_gaussian() to "<GetName() - <<" failed. Fitstatus: "<GetListOfFunctions()->Delete(); } } diff --git a/JetAnalyzers/interface/fitCrystalBall.h b/JetAnalyzers/interface/fitCrystalBall.h new file mode 100644 index 00000000..80341eb3 --- /dev/null +++ b/JetAnalyzers/interface/fitCrystalBall.h @@ -0,0 +1,9 @@ +#ifndef FITCRYSTALBALL_HH +#define FITCRYSTALBALL_HH + +#include +#include + +void fitCrystalBall(TH1F*& hrspOrig, std::string alg, std::string histName, int polDeg, bool normalized, std::string fitDir); + +#endif diff --git a/JetAnalyzers/src/fitCrystalBall.cc b/JetAnalyzers/src/fitCrystalBall.cc new file mode 100644 index 00000000..bd109ab4 --- /dev/null +++ b/JetAnalyzers/src/fitCrystalBall.cc @@ -0,0 +1,1043 @@ +#include "JetMETAnalysis/JetAnalyzers/interface/crystalBall.h" +#include "JetMETAnalysis/JetAnalyzers/interface/crystalBallRes.h" + +#include +#include +#include +#include +#include +#include +#include +#include "TVirtualFitter.h" +#include "TMath.h" +#include "TSpectrum.h" +#include "TRandom3.h" +#include "TCanvas.h" + +#include "TMath.h" +#include +#include +#include +#include +#include +#include +#include + +TH1F *hXaxis( TH1F* hist, double minX, double maxX); + +typedef std::vector vstring; + +using namespace std; + +void fitCrystalBall(TH1F*& hrspOrig, string alg, string histName, int polDeg, bool normalized, string fitDir) +{ + + // + // do no write all the histogram in the output file + // + TH1::AddDirectory(kFALSE); + + // + // select a particular algorithm + // + //vstring valg; + //valg.push_back("ak5tauHPSlooseCombDBcorrAll"); + //valg.push_back("ak5tauHPSlooseCombDBcorrOneProng0Pi0"); + //valg.push_back("ak5tauHPSlooseCombDBcorrOneProng1Pi0"); + //valg.push_back("ak5tauHPSlooseCombDBcorrTwoProng0Pi0"); + //valg.push_back("ak5tauHPSlooseCombDBcorrTwoProng1Pi0"); + //valg.push_back("ak5tauHPSlooseCombDBcorrThreeProng0Pi0"); + //valg.push_back("ak5tauHPSlooseCombDBcorrThreeProng1Pi0"); + //valg.push_back("ak5tauHPSlooseCombDBcorrAlll2"); + //valg.push_back("ak5tauHPSlooseCombDBcorrOneProng0Pi0l2"); + //valg.push_back("ak5tauHPSlooseCombDBcorrOneProng1Pi0l2"); + //valg.push_back("ak5tauHPSlooseCombDBcorrTwoProng0Pi0l2"); + //valg.push_back("ak5tauHPSlooseCombDBcorrTwoProng1Pi0l2"); + //valg.push_back("ak5tauHPSlooseCombDBcorrThreeProng0Pi0l2"); + //valg.push_back("ak5tauHPSlooseCombDBcorrThreeProng1Pi0l2"); + //bool findAlg{false}; + //if( find(valg.begin(), valg.end(), alg ) != valg.end() ) findAlg=true; + //if(!findAlg) return; + /////////////////////////////////////// + + // + // current histogram name + // + string histname{hrspOrig->GetName()}; + /////////////////////////////////////// + + // + // select the relative response histograms for an eta range. + // For time reason, we split the jobs in section of eta + // + if ( histname.find(histName) == string::npos) return; + /////////////////////////////////////// + + // + // select an histogram + // + //vstring vhistname; + //vhistname.push_back("RelRsp_JetEta-0.5to-0.3_RefPt20to22.5"); + //bool findHistname{false}; + //if( find(vhistname.begin(), vhistname.end(), histname) != vhistname.end() ) findHistname=true; + //if(!findHistname) return; + /////////////////////////////////////// + + // + // start fit procedure + // + cout << "\nCalpas Veeken crystalBallFit processing "<Integral("width")>0) hrspOrig->Scale(1/hrspOrig->Integral("width")); + //cout<<"hrsp integral after scaling: "<Integral("width")<Clone(); + + // to avoid pol divergence at low value, set empty bin to 0.5+/-0.5 + for(int i=1; i<=hrsp->GetNbinsX(); i++) if(hrsp->GetBinContent(i)==0){ + if(normalized){hrsp->SetBinContent(i, 0.005); hrsp->SetBinError(i, 0.005);} + else if(hrsp->Integral()<300000){hrsp->SetBinContent(i, 0.5); hrsp->SetBinError(i, 0.5);} + } + + // find the max of the hist that correspond to the tau pic btw [0, 1.05] + double maxBin{0}; + int binmax = -1; + for ( int i=0; i <= hrsp->FindBin(1.05); ++i ){ + double newMaxBin = hrsp->GetBinContent(i); + if ( newMaxBin > maxBin ) { + maxBin = newMaxBin; binmax = i; + } + } + + // pre-fit function + TF1 *testfit {0}; + TF1 *gausLeft{0}; + TF1 *expLeft {0}; + TF1 *polLeft {0}; + TF1 *sumGaus {0}; + // global fit function + TF1 *crystalBallFit{0}; + + // pre-fit definition + string gaus = "[0]*TMath::Exp(-0.5*TMath::Power((x-[1])/[2],2.))"; + string expo = "[0]*TMath::Exp([1]*(x-[2]))"; + string pol; + + // reduce chisquare and Kolmogorov-Smir test + double chi2, ks; + + // map chi2/KS->fit + map mapKsFit; // the best out of the 4 cases + map mapKsFitB; // best from Blue (case 4: Gauss1+Gauss2+Gauss3+pol+exp_left) + map mapKsFitO; // best from Orange (case 3: Gauss1+Gauss2+Gauss3+pol) + map mapKsFitM; // best from Magenta(case 2: Gauss1+Gaus2+pol) + map mapKsFitG; // best from Green (case 1: Gaus1+pol) + map mapChi2Fit; + map mapChi2FitB; + map mapChi2FitO; + map mapChi2FitM; + map mapChi2FitG; + + bool goodFit{false}, nofit{false}; + + if(polDeg==1){ + + TF1 *polRight{0}; + pol = "[0]+[1]*x"; + + int maxPar{25}; + vstring parName {"norm G1", "mean G1", "sigma G1 left", + "!norm Exp left", "slope Exp left", + "!p0 left", "p1 left", "!R1", "!R2", + "sigma G1 right", "norm G2", "!mean G2", "sigma G2", + "norm G3", "!mean G3", "sigma G3", "p0 right", "p1 right", "R3", + "drawPolLeft", "drawExp", "drawG1", "drawG2", "drawG3", "drawPolRight"}; + + double xbinmax = hrsp->GetXaxis()->GetBinUpEdge(binmax); // x position of the max bin + double ybinmax = hrsp->GetBinContent(binmax); + int bindown = binmax-4; + double minGaus = hrsp->GetXaxis()->GetBinLowEdge(bindown); // x_min position of Gauss1 (=R2) + + + + ////////////////////////////////////////////////////////////////////////// + // left side Gauss1 pre-fit + ////////////////////////////////////////////////////////////////////////// + cout<<"\nProcessing fit Gauss1 left..."<< endl; + testfit = new TF1("testfit", gaus.data(), 0., 2.); + testfit->SetParameter(0, hrsp->GetBinContent(binmax)); + testfit->FixParameter(0, hrsp->GetBinContent(binmax)); //!!!! was not fixed before!!!! + testfit->SetParameter(1, xbinmax); + testfit->FixParameter(1, xbinmax); + testfit->SetParameter(2, 0.10); + testfit->SetParLimits(2, 0., 1.e+1); + gausLeft = new TF1("gausLeft", gaus.data(), 0., 2.); + do { + hrsp->Fit(testfit, "Q0", "", minGaus, xbinmax); // Q: quiet mode; 0: do not draw default canvas + // change the current hist axis to match the fit hist axis + TH1F * hrspXaxis = hXaxis(hrsp, minGaus, xbinmax); + TF1 *fit=hrsp->GetFunction("testfit"); + // ensure that the current hist and the fit hist have the same number of point + fit->SetNpx(hrspXaxis->GetNbinsX()); + // compute Kolmogorov and reduce Chi2 test + ks = hrspXaxis->KolmogorovTest((TH1D*)fit->GetHistogram()); + chi2 = testfit->GetChisquare()/max(1, testfit->GetNDF()); + // add a new bin in the fit toward the left for the next fit iteration + --bindown; + minGaus = hrsp->GetXaxis()->GetBinLowEdge(bindown); + if ( chi2 < 5. && ks > 0.01) { // Calpas and Veeken guess + gausLeft->SetParameter(0, testfit->GetParameter(0)); + gausLeft->SetParameter(1, testfit->GetParameter(1)); + gausLeft->SetParameter(2, testfit->GetParameter(2)); + } + } while ( chi2 < 5. && minGaus > 0. && ks > 0.01 ); + // sometimes chi2 may exceed 5 in 1st interation->set default value + if(gausLeft->GetParameter(0)==0) gausLeft->SetParameter(0, testfit->GetParameter(0)); + if(gausLeft->GetParameter(1)==0) gausLeft->SetParameter(1, testfit->GetParameter(1)); + if(gausLeft->GetParameter(2)==0) gausLeft->SetParameter(2, testfit->GetParameter(2)); + //std::cout << "prefit for Gaus1: norm = " << testfit->GetParameter(0) << ", mean = " << testfit->GetParameter(1) << ", sigma = " << testfit->GetParameter(2) << std::endl; + //std::cout << " (minGaus = " << minGaus << ")" << std::endl; + + // go back on bin as it was not take into account + ++bindown; + minGaus = hrsp->GetXaxis()->GetBinCenter(bindown); + + + + ////////////////////////////////////////////////////////////////////////// + // left side expo pre-fit + ////////////////////////////////////////////////////////////////////////// + cout<<"\nProcessing fit expo left..."<< endl; + double minExpLeft = hrsp->GetXaxis()->GetBinLowEdge(bindown-3); // x_min expo position (=R1) + expLeft = new TF1("expLeft", expo.data(), 0., 2.); + if ( minGaus > 0 && hrsp->Integral(hrsp->FindBin(0), hrsp->FindBin(minGaus)) > 0. ){ + delete testfit; + testfit = new TF1("testfit", expo.data(), 0., 2.); + testfit->SetParameter(0, gausLeft->Eval(minGaus)); + testfit->SetParameter(1, 0.1); + testfit->SetParameter(2, minGaus); + testfit->FixParameter(2, minGaus); + do { + hrsp->Fit(testfit, "Q0", "", minExpLeft, minGaus); + --bindown; + minExpLeft = hrsp->GetXaxis()->GetBinCenter(bindown); + chi2 = testfit->GetChisquare()/max(1, testfit->GetNDF()); + if ( chi2 < 5. ) { + expLeft->SetParameter(0, testfit->GetParameter(0)); + expLeft->SetParameter(1, testfit->GetParameter(1)); + expLeft->SetParameter(2, testfit->GetParameter(2)); + } + } + while ( chi2 < 5. && minExpLeft > 0 ); + if(expLeft->GetParameter(0)==0) expLeft->SetParameter(0, testfit->GetParameter(0)); + if(expLeft->GetParameter(1)==0) expLeft->SetParameter(1, testfit->GetParameter(1)); + if(expLeft->GetParameter(2)==0) expLeft->SetParameter(2, testfit->GetParameter(2)); + + ++bindown; + minExpLeft = hrsp->GetXaxis()->GetBinCenter(bindown); + } + else cout<<"no point or and/or range for the expo left fit\n!!"; + + //std::cout << "prefit for Exp: norm = " << testfit->GetParameter(0) << ", slope = " << testfit->GetParameter(1) << ", offset = " << testfit->GetParameter(2) << std::endl; + + + + ////////////////////////////////////////////////////////////////////////// + // left side Pol pre-fit + ////////////////////////////////////////////////////////////////////////// + cout<<"\nProcessing fit pol left...\n"; + if( minExpLeft > 0 && hrsp->Integral(hrsp->FindBin(0), hrsp->FindBin(minExpLeft)) > 0. ){ + polLeft = new TF1("polLeft", pol.c_str(), 0, minExpLeft); + } + else cout<<"no point or and/or range for the pol left fit!!\n"; + + //std::cout << "prefit for Pol(left):" + //for(int i=0; i<2; i++) cout<<"p["<GetParameter(i) <GetXaxis()->GetBinUpEdge(binup); + delete testfit; + testfit = new TF1("sumGaus", crystalBall_1, 0., 2., maxPar); + testfit->SetParameter(0, 0.70*hrsp->GetBinContent(binmax)); + testfit->SetParameter(1, xbinmax); + testfit->SetParameter(2, 0.1); + testfit->SetParameter(3, 0.2*hrsp->GetBinContent(binmax)); + testfit->SetParameter(4, xbinmax+0.1); + testfit->SetParameter(5, 0.15); + testfit->SetParameter(6, 0.1*hrsp->GetBinContent(binmax)); + testfit->SetParameter(7, xbinmax+0.2); + testfit->SetParameter(8, 0.15); + sumGaus = new TF1("sumGaus", crystalBall_1, 0., 2., maxPar); + do { + hrsp->Fit(testfit, "Q0", "", xbinmax, maxSumGaus); + ++binup; + maxSumGaus = hrsp->GetXaxis()->GetBinCenter(binup); + chi2 = testfit->GetChisquare()/max(1, testfit->GetNDF()); + if ( chi2 < 5. ) for(int i=9; i<16; i++) sumGaus->SetParameter(i, testfit->GetParameter(i)); + } + while ( chi2 < 5. && maxSumGaus < 2. ); + for(int i=9; i<16; i++) if(sumGaus->GetParameter(i)==0) sumGaus->SetParameter(i, testfit->GetParameter(i-9)); + --binup; + maxSumGaus = hrsp->GetXaxis()->GetBinCenter(binup); + + + + ////////////////////////////////////////////////////////////////////////// + // Right side Pol pre-fit + ////////////////////////////////////////////////////////////////////////// + cout<<"\nProcessing fit pol right...\n"; + if( maxSumGaus < 2 && hrsp->Integral(hrsp->FindBin(maxSumGaus), hrsp->FindBin(2)) > 0. ){ + polRight = new TF1("polRight", pol.c_str(), maxSumGaus, 2); + } + else cout<<"no point or and/or range for the pol right fit!!\n"; + //std::cout << "prefit for Pol(right):" + //for(int i=0; i<2; i++) cout<<"p["<GetParameter(i) <GetParameter(1); + double mu2=0.15; + double mu3=0.15; + double r1=minGaus-minExpLeft; + double r2=xbinmax-minGaus; + //double r3=maxSumGaus-xbinmax; + double r3=maxSumGaus; + //cout<<"maxsumGaus: "< mrd {rd_unfix, rd_fixG1, rd_fixG2, rd_fixG3, rd_fixR1, rd_fixR2, rd_fixR3}; //7*12 + + // + // try ideal + Random value and after 100 try, select the best fit + // + TH1F* hrspClone; + + TRandom *random = new TRandom3(); + cout<<"\nStart Random...\n"; + for(int i=0; i<300; i++){ + //for(int i=1; i<2; i++){ + cout<<"\ntest number: "<SetLineWidth(2); + for(int par=0; parSetParameter(par, 0); + crystalBallFit->SetParName(par, parName[par].c_str()); + } + + + + //////////////// + // fit case 4 + //////////////// + crystalBallFit->SetLineColor(kBlue); + + if(gausLeft){ + double stepG1{0.1*random->Rndm()}; // btw [0, 1] + double stepR2{0.1*random->Rndm()}; + + if(i<7){ // "optimal value = 6 1st vector of the matrix" + if(mrd[i][6]==0.) crystalBallFit->SetParameter(1, mrd[i][0]); // mean + else crystalBallFit->FixParameter(1, mrd[i][0]); + if(mrd[i][10]==0.) crystalBallFit->SetParameter(8, mrd[i][4]); // R2 + else crystalBallFit->FixParameter(8, mrd[i][4]); + } + //else if(i>=7 && i<53){ // random up from optimal value + else if(i>=7 && i<146){ // random up from optimal value + if(mrd[0][6]==0.) crystalBallFit->SetParameter(1, mrd[0][0]+stepG1); //mean + else crystalBallFit->FixParameter(1, mrd[0][0]+stepG1); + if(mrd[0][10]==0.) crystalBallFit->SetParameter(8, mrd[0][4]+stepR2); // R2 + else crystalBallFit->FixParameter(8, mrd[0][4]+stepR2); + } + else{ // random down from optimal value + if(mrd[0][6]==0.) crystalBallFit->SetParameter(1, mrd[0][0]-stepG1); //mean + else crystalBallFit->FixParameter(1, mrd[0][0]-stepG1); + if(mrd[0][10]==0.) crystalBallFit->SetParameter(8, mrd[0][4]-stepR2); // R2 + else crystalBallFit->FixParameter(8, mrd[0][4]-stepR2); + } + crystalBallFit->SetParameter(0, gausLeft->GetParameter(0)); // norm + //crystalBallFit->SetParLimits(0, 0.5*ybinmax, ybinmax*2); // was not set!! + crystalBallFit->SetParLimits(0, 0., 1000); // was not set!! + crystalBallFit->SetParLimits(1, 0., xbinmax); + crystalBallFit->SetParameter(2, gausLeft->GetParameter(2)); // sigma Gauss1 left + //crystalBallFit->SetParLimits(2, 0.02, 10.); + crystalBallFit->SetParLimits(2, 0.02, 1.); + crystalBallFit->SetParLimits(8, 0., xbinmax); // R2 + }else{ + for(int i=0; i<9; i++) crystalBallFit->FixParameter(i, 0); + } + if(expLeft){ + double stepR1{0.1*random->Rndm()}; + if(i<7){ // "optimal value" + if(mrd[i][9]==0.) crystalBallFit->SetParameter(7, mrd[i][3]); // R1 + else crystalBallFit->FixParameter(7, mrd[i][3]); + } + else if (i>=7 && i<53){ // random up from optimal value + if(mrd[0][9]==0.) crystalBallFit->SetParameter(7, mrd[0][3]+stepR1); // R1 + else crystalBallFit->FixParameter(7, mrd[0][3]+stepR1); + } + else{ // random down from optimal value + if(mrd[0][9]==0.) crystalBallFit->SetParameter(7, mrd[0][3]-stepR1); // R1 + else crystalBallFit->FixParameter(7, mrd[0][3]-stepR1); + } + crystalBallFit->FixParameter(3, 0.); // norm + crystalBallFit->SetParameter(4, expLeft->GetParameter(1)); // slope + //crystalBallFit->SetParLimits(4, 0., 1.e+4); + crystalBallFit->SetParLimits(4, 0.01, 20); + crystalBallFit->SetParLimits(7, 0., xbinmax); + }else{ + for(int i=3; i<9; i++) crystalBallFit->FixParameter(i, 0); + } + if(polLeft){ + crystalBallFit->FixParameter(5, 0); // P0 left + crystalBallFit->SetParameter(6, polLeft ->GetParameter(1)); + }else{ + for(int i=5; i<8; i++) crystalBallFit->FixParameter(i, 0); + } + // + if(polRight){ + crystalBallFit->FixParameter(16, 0); // P0 right + crystalBallFit->SetParameter(17, polRight ->GetParameter(1)); + }else{ + for(int i=16; i<19; i++) crystalBallFit->FixParameter(i, 0); + } + + double stepG2{0.1*random->Rndm()}; + double stepG3{0.1*random->Rndm()}; + double stepR3{0.1*random->Rndm()}; + + //if(i!=269) continue;//!!!!!!!!!! to select only the problematic + + if(i<7){ // "optimal value" + if(mrd[i][7]==0.) crystalBallFit->SetParameter(11, mrd[i][1]); // Gauss2 mean + else crystalBallFit->FixParameter(11, mrd[i][1]); + if(mrd[i][8]==0.) crystalBallFit->SetParameter(14, mrd[i][2]); // Gauss3 mean + else crystalBallFit->FixParameter(14, mrd[i][2]); + if(mrd[i][11]==0.) crystalBallFit->SetParameter(18, mrd[i][5]); // R3 + else crystalBallFit->FixParameter(18, mrd[i][5]); + } + //else if(i>=7 && i<53){ // random up from optimal value + else if(i>=7 && i<146){ // random up from optimal value + if(mrd[0][7]==0.) crystalBallFit->SetParameter(11, mrd[0][1]+stepG2); + else crystalBallFit->FixParameter(11, mrd[0][1]+stepG2); + if(mrd[0][8]==0.) crystalBallFit->SetParameter(14, mrd[0][2]+stepG3); + else crystalBallFit->FixParameter(14, mrd[0][2]+stepG3); + if(mrd[0][11]==0.) crystalBallFit->SetParameter(18, mrd[0][5]+stepR3); + else crystalBallFit->FixParameter(18, mrd[0][5]+stepR3); + } + else{ // random down from optimal value + if(mrd[0][7]==0.) crystalBallFit->SetParameter(11, mrd[0][1]-stepG2); + else crystalBallFit->FixParameter(11, mrd[0][1]+stepG2); + if(mrd[0][8]==0.) crystalBallFit->SetParameter(14, mrd[0][2]-stepG3); + else crystalBallFit->FixParameter(14, mrd[0][2]+stepG3); + if(mrd[0][11]==0.) crystalBallFit->SetParameter(18, mrd[0][5]-stepR3); + else crystalBallFit->FixParameter(18, mrd[0][5]-stepR3); + } + crystalBallFit->SetParameter(9, 0.05); // sigma Gauss 1 right + //crystalBallFit->SetParLimits(9, 0.01, 10); + crystalBallFit->SetParLimits(9, 0.01, 1); + crystalBallFit->SetParameter(10, 0.5*ybinmax); // Gauss 2 norm + //crystalBallFit->SetParLimits(10, 0., ybinmax); + crystalBallFit->SetParLimits(10, 0., 1000); + crystalBallFit->SetParLimits(11, 0., 2.); + crystalBallFit->SetParameter(12, 0.2); // Gauss 2 sigma + //crystalBallFit->SetParLimits(12, 0., 10.); + crystalBallFit->SetParLimits(12, 0., 1.); + crystalBallFit->SetParameter(13, 0.25*ybinmax); // gauss 3 norm + //crystalBallFit->SetParLimits(13, 0., ybinmax); + crystalBallFit->SetParLimits(13, 0., 1000); + crystalBallFit->SetParLimits(14, 0., 2.); + crystalBallFit->SetParameter(15, 0.1); // Gauss 3 sigma + //crystalBallFit->SetParLimits(15, 0., 10.); + crystalBallFit->SetParLimits(15, 0., 1.); + // bool parameter to draw each function separatly + crystalBallFit->FixParameter(19, 1); // draw pol left + crystalBallFit->FixParameter(20, 1); // draw exp + crystalBallFit->FixParameter(21, 1); // draw g1 + crystalBallFit->FixParameter(22, 1); // draw g2 + crystalBallFit->FixParameter(23, 1); // draw g3 + crystalBallFit->FixParameter(24, 1); // draw pol right + // + hrsp->Fit("crystalBallFit", "R"); + testfit = (TF1*) hrsp->GetListOfFunctions()->Last(); + testfit->SetNpx(hrsp->GetNbinsX()); + ks=hrsp->KolmogorovTest(testfit->GetHistogram()); + chi2=testfit->GetChisquare()/max(1, testfit->GetNDF()); + cout<<"ks/chi2 case4: "<Clone(); + if( ks>1e-4 && fabs(ks-1)>0.1 && testfit->Eval(1)>0) mapKsFitB[ks]=hrspClone; // 0 and 1 are bad ks values 1e-7!! + if(chi2<50 && testfit->Eval(1)>0) mapChi2FitB[chi2]=hrspClone; + + + //////////////// + // fit case 3 + //////////////// + crystalBallFit->SetLineColor(kOrange); + crystalBallFit->FixParameter(6, 0); // pol7/range left=cst p0 + hrsp->Fit("crystalBallFit", "R"); + testfit = (TF1*) hrsp->GetListOfFunctions()->Last(); + testfit->SetNpx(hrsp->GetNbinsX()); + ks=hrsp->KolmogorovTest(testfit->GetHistogram()); + chi2=testfit->GetChisquare()/max(1, testfit->GetNDF()); + cout<<"ks/chi2 case3: "<Clone(); + if( ks>1e-4 && fabs(ks-1)>0.1 && testfit->Eval(1)>0) mapKsFitO[ks]=hrspClone; + if(chi2<50 && testfit->Eval(1)>0) mapChi2FitO[chi2]=hrspClone; + + + + //////////////// + // fit case 2 + //////////////// + crystalBallFit->SetLineColor(kMagenta); + for(int i=13; i<15; i++) crystalBallFit->FixParameter(i, 0); // gaus3 right + crystalBallFit->FixParameter(15, 1); // sigma !=0 + hrsp->Fit("crystalBallFit", "R"); + testfit = (TF1*) hrsp->GetListOfFunctions()->Last(); + testfit->SetNpx(hrsp->GetNbinsX()); + ks=hrsp->KolmogorovTest(testfit->GetHistogram()); + chi2=testfit->GetChisquare()/max(1, testfit->GetNDF()); + cout<<"ks/chi2 case2 "<Clone(); + if(ks>1e-4 && fabs(ks-1)>0.1 && testfit->Eval(1)>0) mapKsFitM[ks]=hrspClone; + if(chi2<50 && testfit->Eval(1)>0) mapChi2FitM[chi2]=hrspClone; + + + + //////////////// + // fit case 1 + //////////////// + crystalBallFit->SetLineColor(kGreen); + for(int i=10; i<12; i++) crystalBallFit->FixParameter(i, 0); // gaus2 right + crystalBallFit->FixParameter(12, 1); // sigma !=0 + hrsp->Fit("crystalBallFit", "R"); + testfit = (TF1*) hrsp->GetListOfFunctions()->Last(); + testfit->SetNpx(hrsp->GetNbinsX()); + ks=hrsp->KolmogorovTest(testfit->GetHistogram()); + chi2=testfit->GetChisquare()/max(1, testfit->GetNDF()); + cout<<"ks/chi2 case1 "<Clone(); + if( ks>1e-4 && fabs(ks-1)>0.1 && testfit->Eval(1)>0) mapKsFitG[ks]=hrspClone; + if(chi2<50 && testfit->Eval(1)>0) mapChi2FitG[chi2]=hrspClone; + } // random + cout<<"endRandom\n"; + cout<<"\nretriving best fit from ks or chi2 map..."<first]=mapKsFitB.rbegin()->second; + if(mapKsFitO.size()!=0)mapKsFit[mapKsFitO.rbegin()->first]=mapKsFitO.rbegin()->second; + if(mapKsFitM.size()!=0)mapKsFit[mapKsFitM.rbegin()->first]=mapKsFitM.rbegin()->second; + if(mapKsFitG.size()!=0)mapKsFit[mapKsFitG.rbegin()->first]=mapKsFitG.rbegin()->second; + // + if(mapChi2FitB.size()!=0)mapChi2Fit[mapChi2FitB.begin()->first]=mapChi2FitB.begin()->second; + if(mapChi2FitO.size()!=0)mapChi2Fit[mapChi2FitO.begin()->first]=mapChi2FitO.begin()->second; + if(mapChi2FitM.size()!=0)mapChi2Fit[mapChi2FitM.begin()->first]=mapChi2FitM.begin()->second; + if(mapChi2FitG.size()!=0)mapChi2Fit[mapChi2FitG.begin()->first]=mapChi2FitG.begin()->second; + + // select best fit from ks then chi2, if none best failed ks then chi2, if none keep original hist + bool selKs{false}, selChi2{false}; + + if(mapKsFit.size()!=0 || mapChi2Fit.size()!=0){ + if(mapChi2Fit.size()!=0){ + double minChi2 {mapChi2Fit.begin()->first}; + if(minChi2<5){ + goodFit=true; hrspOrig= (TH1F*) mapChi2Fit.begin()->second->Clone(); + cout<<"pass Chi2: "<first<first}; + if(maxKs>0.5){ + goodFit=true; hrspOrig= (TH1F*) mapKsFit.rbegin()->second->Clone(); + cout<<"pass Ks: "<first<second->Clone(); + cout<<"failed Chi2: "<first<second->Clone(); + cout<<"failed Ks: "<first<Write(); + } + else{ + + pol = "[0]+[1]*x+[2]*TMath::Power(x,2)+[3]*TMath::Power(x,3)+" + "[4]*TMath::Power(x,4)+[5]*TMath::Power(x,5)+" + "[6]*TMath::Power(x,6)+[7]*TMath::Power(x,7)"; + + // global fit maximun parameters and name + int maxPar{36}; + vstring parName {"norm G1", "mean G1", "sigma G1 left", "!norm Exp left", "slope Exp left", + "!p0 left", "p1 left", "p2 left", "p3 left", "p4 left", "p5 left", "p6 left", "p7 left", + "!R1", "!R2", "sigma G1 right", "norm G2", "!mean G2", "sigma G2", "norm G3", "!mean G3", "sigma G3", + "p0 right", "p1 right", "p2 right", "p3 right", "p4 right", "p5 right", "p6 right", "p7 right", + "drawPolLeft", "drawExp", "drawG1", "drawG2", "drawG3", "drawPolRight"}; + + + double xbinmax = hrsp->GetXaxis()->GetBinUpEdge(binmax); // x position of the max bin + double ybinmax = hrsp->GetBinContent(binmax); + int bindown = binmax-4; + double minGaus = hrsp->GetXaxis()->GetBinLowEdge(bindown); // x_min position of Gauss1 (=R2) + + + ////////////////////////////////////////////////////////////////////////// + // left side Gauss1 pre-fit + ////////////////////////////////////////////////////////////////////////// + cout<<"\nProcessing fit Gauss1 left..."<< endl; + testfit = new TF1("testfit", gaus.data(), 0., 2.); + testfit->SetParameter(0, hrsp->GetBinContent(binmax)); + testfit->FixParameter(0, hrsp->GetBinContent(binmax)); //!!!! was not fixed before!!!! + testfit->SetParameter(1, xbinmax); + testfit->FixParameter(1, xbinmax); + testfit->SetParameter(2, 0.10); + testfit->SetParLimits(2, 0., 1.e+1); + gausLeft = new TF1("gausLeft", gaus.data(), 0., 2.); + do { + hrsp->Fit(testfit, "Q0", "", minGaus, xbinmax); // Q: quiet mode; 0: do not draw default canvas + // change the current hist axis to match the fit hist axis + TH1F * hrspXaxis = hXaxis(hrsp, minGaus, xbinmax); + TF1 *fit=hrsp->GetFunction("testfit"); + // ensure that the current hist and the fit hist have the same number of point + fit->SetNpx(hrspXaxis->GetNbinsX()); + // compute Kolmogorov and reduce Chi2 test + ks = hrspXaxis->KolmogorovTest((TH1D*)fit->GetHistogram()); + chi2 = testfit->GetChisquare()/max(1, testfit->GetNDF()); + // add a new bin in the fit toward the left for the next fit iteration + --bindown; + minGaus = hrsp->GetXaxis()->GetBinLowEdge(bindown); + if ( chi2 < 5. && ks > 0.01) { // Calpas and Veeken guess + gausLeft->SetParameter(0, testfit->GetParameter(0)); + gausLeft->SetParameter(1, testfit->GetParameter(1)); + gausLeft->SetParameter(2, testfit->GetParameter(2)); + } + } while ( chi2 < 5. && minGaus > 0. && ks > 0.01 ); + // sometimes chi2 may exceed 5 in 1st interation->set default value + if(gausLeft->GetParameter(0)==0) gausLeft->SetParameter(0, testfit->GetParameter(0)); + if(gausLeft->GetParameter(1)==0) gausLeft->SetParameter(1, testfit->GetParameter(1)); + if(gausLeft->GetParameter(2)==0) gausLeft->SetParameter(2, testfit->GetParameter(2)); + //std::cout << "prefit for Gaus1: norm = " << testfit->GetParameter(0) << ", mean = " << testfit->GetParameter(1) << ", sigma = " << testfit->GetParameter(2) << std::endl; + //std::cout << " (minGaus = " << minGaus << ")" << std::endl; + + // go back on bin as it was not take into account + ++bindown; + minGaus = hrsp->GetXaxis()->GetBinCenter(bindown); + + + ////////////////////////////////////////////////////////////////////////// + // left side expo pre-fit + ////////////////////////////////////////////////////////////////////////// + cout<<"\nProcessing fit expo left..."<< endl; + double minExpLeft = hrsp->GetXaxis()->GetBinLowEdge(bindown-3); // x_min expo position (=R1) + expLeft = new TF1("expLeft", expo.data(), 0., 2.); + if ( minGaus > 0 && hrsp->Integral(hrsp->FindBin(0), hrsp->FindBin(minGaus)) > 0. ){ + delete testfit; + testfit = new TF1("testfit", expo.data(), 0., 2.); + testfit->SetParameter(0, gausLeft->Eval(minGaus)); + testfit->SetParameter(1, 0.1); + testfit->SetParameter(2, minGaus); + testfit->FixParameter(2, minGaus); + do { + hrsp->Fit(testfit, "Q0", "", minExpLeft, minGaus); + --bindown; + minExpLeft = hrsp->GetXaxis()->GetBinCenter(bindown); + chi2 = testfit->GetChisquare()/max(1, testfit->GetNDF()); + if ( chi2 < 5. ) { + expLeft->SetParameter(0, testfit->GetParameter(0)); + expLeft->SetParameter(1, testfit->GetParameter(1)); + expLeft->SetParameter(2, testfit->GetParameter(2)); + } + } + while ( chi2 < 5. && minExpLeft > 0 ); + if(expLeft->GetParameter(0)==0) expLeft->SetParameter(0, testfit->GetParameter(0)); + if(expLeft->GetParameter(1)==0) expLeft->SetParameter(1, testfit->GetParameter(1)); + if(expLeft->GetParameter(2)==0) expLeft->SetParameter(2, testfit->GetParameter(2)); + + ++bindown; + minExpLeft = hrsp->GetXaxis()->GetBinCenter(bindown); + } + else { + cout<<"no point or and/or range for the expo left fit\n!!"; + } + //std::cout << "prefit for Exp: norm = " << testfit->GetParameter(0) << ", slope = " << testfit->GetParameter(1) << ", offset = " << testfit->GetParameter(2) << std::endl; + + + ////////////////////////////////////////////////////////////////////////// + // left side Pol pre-fit + ////////////////////////////////////////////////////////////////////////// + cout<<"\nProcessing fit pol left...\n"; + if( minExpLeft > 0 && hrsp->Integral(hrsp->FindBin(0), hrsp->FindBin(minExpLeft)) > 0. ){ + polLeft = new TF1("polLeft", pol.c_str(), 0, minExpLeft); + } + else { + cout<<"no point or and/or range for the pol7 left fit!!\n"; + } + //std::cout << "prefit for Pol(left):" + //if(pol1) for(int i=0; i<8; i++) cout<<"p["<GetParameter(i) <SetLineWidth(2); + for(int par=0; parSetParameter(par, 0); + crystalBallFit->SetParName(par, parName[par].c_str()); + } + + //////////////// + // fit case 4 + //////////////// + crystalBallFit->SetLineColor(kBlue); + + if(gausLeft){ + double stepG1{0.1*random->Rndm()}; + double stepR2{0.1*random->Rndm()}; + + crystalBallFit->SetParameter(0, gausLeft->GetParameter(0)); // norm + //crystalBallFit->SetParLimits(0, 0., ybinmax*2); // was not set!! + if(i<6){ // "optimal value = 6 1st vector of the matrix" + if(mrd[i][5]==0.) { crystalBallFit->SetParameter(1, mrd[i][0]); crystalBallFit->SetParLimits(1, 0, xbinmax);} // mean + else crystalBallFit->FixParameter(1, mrd[i][0]); + if(mrd[i][9]==0.) crystalBallFit->SetParameter(14, mrd[i][4]); // R2 + else crystalBallFit->FixParameter(14, mrd[i][4]); + } + else if(i>=6 && i<47){ // random up from optimal value + if(mrd[0][5]==0.) { crystalBallFit->SetParameter(1, mrd[0][0]+stepG1); crystalBallFit->SetParLimits(1, 0, xbinmax);} //mean + else crystalBallFit->FixParameter(1, mrd[0][0]+stepG1); + if(mrd[0][9]==0.) crystalBallFit->SetParameter(14, mrd[0][4]+stepR2); // R2 + else crystalBallFit->FixParameter(14, mrd[0][4]+stepR2); + } + else{ // random down from optimal value + if(mrd[0][5]==0.) { crystalBallFit->SetParameter(1, mrd[0][0]-stepG1); crystalBallFit->SetParLimits(1, 0, xbinmax);} //mean + else crystalBallFit->FixParameter(1, mrd[0][0]-stepG1); + if(mrd[0][9]==0.) crystalBallFit->SetParameter(14, mrd[0][4]-stepR2); // R2 + else crystalBallFit->FixParameter(14, mrd[0][4]-stepR2); + } + crystalBallFit->SetParameter(2, gausLeft->GetParameter(2)); // sigma Gauss1 left + //crystalBallFit->SetParLimits(2, 0.02, 10.); + crystalBallFit->SetParLimits(2, 0.02, 1.); + crystalBallFit->SetParLimits(14, 0., xbinmax); + }else{ + for(int i=0; i<15; i++) crystalBallFit->FixParameter(i, 0); + } + if(expLeft){ + double stepR1{0.1*random->Rndm()}; + crystalBallFit->FixParameter(3, 0.); // norm + crystalBallFit->SetParameter(4, expLeft->GetParameter(1)); // slope + //crystalBallFit->SetParLimits(4, 0., 1.e+4); + crystalBallFit->SetParLimits(4, 0.01, 20); + if(i<6){ // "optimal value" + if(mrd[i][8]==0.) crystalBallFit->SetParameter(13, mrd[i][3]); // R1 + else crystalBallFit->FixParameter(13, mrd[i][3]); + } + else if (i>=6 && i<47){ // random up from optimal value + if(mrd[0][8]==0.) crystalBallFit->SetParameter(13, mrd[0][3]+stepR1); // R1 + else crystalBallFit->FixParameter(13, mrd[0][3]+stepR1); + } + else{ // random down from optimal value + if(mrd[0][8]==0.) crystalBallFit->SetParameter(13, mrd[0][3]-stepR1); // R1 + else crystalBallFit->FixParameter(13, mrd[0][3]-stepR1); + } + crystalBallFit->SetParLimits(13, 0., 1.); + }else{ + for(int i=3; i<15; i++) crystalBallFit->FixParameter(i, 0); + } + if(polLeft){ + for(int i=1; i<8; i++) crystalBallFit->SetParameter(i+5, polLeft ->GetParameter(i)); + crystalBallFit->FixParameter(5, 0); // P0 left + for(int i=0; i<7; i++) crystalBallFit->SetParLimits(i+6, -100, 100 ); //!!!!?? + }else{ + for(int i=5; i<=13; i++) crystalBallFit->FixParameter(i, 0); + } + // + double stepG2{0.1*random->Rndm()}; + double stepG3{0.1*random->Rndm()}; + + //if(i!=96) continue;//!!!!!!!!!! to select only the problematic + + if(i<6){ // "optimal value" + if(mrd[i][6]==0.) crystalBallFit->SetParameter(17, mrd[i][1]); // Gauss2 mean + else crystalBallFit->FixParameter(17, mrd[i][1]); + if(mrd[i][7]==0.) crystalBallFit->SetParameter(20, mrd[i][2]); // Gauss3 mean + else crystalBallFit->FixParameter(20, mrd[i][2]); + } + else if(i>=6 && i<47){ // random up from optimal value + if(mrd[0][6]==0.) crystalBallFit->SetParameter(17, mrd[0][1]+stepG2); + else crystalBallFit->FixParameter(17, mrd[0][1]+stepG2); + if(mrd[0][7]==0.) crystalBallFit->SetParameter(20, mrd[0][2]+stepG3); + else crystalBallFit->FixParameter(20, mrd[0][2]+stepG3); + } + else{ // random down from optimal value + if(mrd[0][6]==0.) crystalBallFit->SetParameter(17, mrd[0][1]-stepG2); + else crystalBallFit->FixParameter(17, mrd[0][1]+stepG2); + if(mrd[0][7]==0.) crystalBallFit->SetParameter(20, mrd[0][2]-stepG3); + else crystalBallFit->FixParameter(20, mrd[0][2]+stepG3); + } + crystalBallFit->SetParameter(15, 0.05); // sigma Gauss 1 right + //crystalBallFit->SetParLimits(15, 0.01, 10); + crystalBallFit->SetParLimits(15, 0.01, 1); + crystalBallFit->SetParameter(16, 0.01*ybinmax); // Gauss 2 norm + crystalBallFit->SetParLimits(16, 0., ybinmax); + crystalBallFit->SetParLimits(17, 0., 2.); + crystalBallFit->SetParameter(18, 0.2); // Gauss 2 sigma + //crystalBallFit->SetParLimits(18, 0., 10.); + crystalBallFit->SetParLimits(18, 0., 1.); + crystalBallFit->SetParameter(19, 0.10*ybinmax); // gauss 3 norm + crystalBallFit->SetParLimits(19, 0., ybinmax); + crystalBallFit->SetParLimits(20, 0., 2.); + crystalBallFit->SetParameter(21, 0.1); // Gauss 3 sigma + //crystalBallFit->SetParLimits(21, 0., 10.); + crystalBallFit->SetParLimits(21, 0., 1.); + crystalBallFit->SetParameter(22, 0.5); // P0 right + crystalBallFit->SetParLimits(22, 0., ybinmax); + for(int i=0; i<8; i++) crystalBallFit->SetParLimits(i+22, -100, 100 ); //!!!!?? + if(polDeg==1){ + crystalBallFit->SetParameter(23, -0.1); + crystalBallFit->SetParLimits(23, -1, 0); + for(int i=0; i<6; i++) crystalBallFit->FixParameter(i+24, 0); + } + // bool parameter to draw each function separatly + crystalBallFit->FixParameter(30, 1); // draw pol left + crystalBallFit->FixParameter(31, 1); // draw exp + crystalBallFit->FixParameter(32, 1); // draw g1 + crystalBallFit->FixParameter(33, 1); // draw g2 + crystalBallFit->FixParameter(34, 1); // draw g3 + crystalBallFit->FixParameter(35, 1); // draw pol right + // + hrsp->Fit("crystalBallFit", "R"); + testfit = (TF1*) hrsp->GetListOfFunctions()->Last(); + testfit->SetNpx(hrsp->GetNbinsX()); + ks=hrsp->KolmogorovTest(testfit->GetHistogram()); + chi2=testfit->GetChisquare()/max(1, testfit->GetNDF()); + cout<<"ks/chi2 case4: "<Clone(); + if( ks>1e-4 && fabs(ks-1)>0.1 && testfit->Eval(1)>0) mapKsFitB[ks]=hrspClone; // 0 and 1 are bad ks values 1e-7!! + if(chi2<50 && testfit->Eval(1)>0) mapChi2FitB[chi2]=hrspClone; + + + //////////////// + // fit case 3 + //////////////// + crystalBallFit->SetLineColor(kOrange); + if(polLeft)for(int i=6; i<=12; i++) crystalBallFit->FixParameter(i, 0); // pol7/range left=cst p0 + for(int i=23; i<=29; i++) crystalBallFit->FixParameter(i, 0); // pol7 right= cst p0 + hrsp->Fit("crystalBallFit", "R"); + testfit = (TF1*) hrsp->GetListOfFunctions()->Last(); + testfit->SetNpx(hrsp->GetNbinsX()); + ks=hrsp->KolmogorovTest(testfit->GetHistogram()); + chi2=testfit->GetChisquare()/max(1, testfit->GetNDF()); + cout<<"ks/chi2 case3: "<Clone(); + if( ks>1e-4 && fabs(ks-1)>0.1 && testfit->Eval(1)>0) mapKsFitO[ks]=hrspClone; + if(chi2<50 && testfit->Eval(1)>0) mapChi2FitO[chi2]=hrspClone; + + + + //////////////// + // fit case 2 + //////////////// + crystalBallFit->SetLineColor(kMagenta); + for(int i=19; i<=20; i++) crystalBallFit->FixParameter(i, 0); // gaus3 right + crystalBallFit->FixParameter(21, 1); // sigma !=0 + hrsp->Fit("crystalBallFit", "R"); + testfit = (TF1*) hrsp->GetListOfFunctions()->Last(); + testfit->SetNpx(hrsp->GetNbinsX()); + ks=hrsp->KolmogorovTest(testfit->GetHistogram()); + chi2=testfit->GetChisquare()/max(1, testfit->GetNDF()); + cout<<"ks/chi2 case2 "<Clone(); + if(ks>1e-4 && fabs(ks-1)>0.1 && testfit->Eval(1)>0) mapKsFitM[ks]=hrspClone; + if(chi2<50 && testfit->Eval(1)>0) mapChi2FitM[chi2]=hrspClone; + + + + //////////////// + // fit case 1 + //////////////// + crystalBallFit->SetLineColor(kGreen); + for(int i=16; i<=17; i++) crystalBallFit->FixParameter(i, 0); // gaus2 right + crystalBallFit->FixParameter(18, 1); // sigma !=0 + hrsp->Fit("crystalBallFit", "R"); + testfit = (TF1*) hrsp->GetListOfFunctions()->Last(); + testfit->SetNpx(hrsp->GetNbinsX()); + ks=hrsp->KolmogorovTest(testfit->GetHistogram()); + chi2=testfit->GetChisquare()/max(1, testfit->GetNDF()); + cout<<"ks/chi2 case1 "<Clone(); + if( ks>1e-4 && fabs(ks-1)>0.1 && testfit->Eval(1)>0) mapKsFitG[ks]=hrspClone; + if(chi2<50 && testfit->Eval(1)>0) mapChi2FitG[chi2]=hrspClone; + + } // random + cout<<"endRandom\n"; + cout<<"\nretriving best fit from ks or chi2 map..."<first]=mapKsFitB.rbegin()->second; + if(mapKsFitO.size()!=0)mapKsFit[mapKsFitO.rbegin()->first]=mapKsFitO.rbegin()->second; + if(mapKsFitM.size()!=0)mapKsFit[mapKsFitM.rbegin()->first]=mapKsFitM.rbegin()->second; + if(mapKsFitG.size()!=0)mapKsFit[mapKsFitG.rbegin()->first]=mapKsFitG.rbegin()->second; + // + if(mapChi2FitB.size()!=0)mapChi2Fit[mapChi2FitB.begin()->first]=mapChi2FitB.begin()->second; + if(mapChi2FitO.size()!=0)mapChi2Fit[mapChi2FitO.begin()->first]=mapChi2FitO.begin()->second; + if(mapChi2FitM.size()!=0)mapChi2Fit[mapChi2FitM.begin()->first]=mapChi2FitM.begin()->second; + if(mapChi2FitG.size()!=0)mapChi2Fit[mapChi2FitG.begin()->first]=mapChi2FitG.begin()->second; + + // select best fit from ks then chi2, if none best failed ks then chi2, if none keep original hist + bool selKs{false}, selChi2{false}; + + if(mapKsFit.size()!=0 || mapChi2Fit.size()!=0){ + if(mapKsFit.size()!=0){ + double maxKs {mapKsFit.rbegin()->first}; + if(maxKs>0.5){ + goodFit=true; hrspOrig= (TH1F*) mapKsFit.rbegin()->second->Clone(); + cout<<"pass Ks: "<first<first}; + if(minChi2<5){ + goodFit=true; hrspOrig= (TH1F*) mapChi2Fit.begin()->second->Clone(); + cout<<"pass Chi2: "<first<second->Clone(); + cout<<"failed Ks: "<first<second->Clone(); + cout<<"failed Chi2: "<first<Write(); + + } // else polDeg=7 + + // + // save fit plot + // + TCanvas *canvas = new TCanvas("canvas", "", 700, 700); + canvas->Divide(1,2); + canvas->cd(1); + TH1F *hrspOrigClone=(TH1F*)hrspOrig->Clone(); + hrspOrigClone->GetXaxis()->SetRangeUser(0.8, 1.2); + hrspOrigClone->Draw(); + canvas->cd(2); + gPad->SetLogy(); + hrspOrig->Draw(); + + // + // print CB fit result if fit exist + // + if(!nofit){ + cout<<"\nretrieving parameters from the selected fit... "< vintegral(3), vnorm(3), vrange(3), vmean(3), vsigma(4), vexp(2), vpolLeft(8), vpolRight(8); + vector voption(6); + + crystalBallRes(hrspOrig, true, vintegral, vnorm, vrange, vmean, vsigma, vexp, vpolLeft, vpolRight, voption, polDeg, canvas); // hist fit, drawfit, intg2, g3, g3, canvas(defaut) + + } + + cout<<"saving fit plot...\n"; + string savePath, dir, hrspOrigName{hrspOrig->GetName()}; + if(goodFit) dir="pass"; else dir="failed"; + savePath=fitDir+"/"+alg+"/crystalBallFit/"+dir+"/"+hrspOrigName+".png"; + + canvas->Print((savePath).c_str()); +} + + +TH1F *hXaxis( TH1F* hist, double minX, double maxX) +{ + double *xbins = new double [hist->GetNbinsX()+1]; + int nbins{0}; + for(int i=0; iGetNbinsX(); i++){ + double x {hist->GetXaxis()->GetBinCenter(i+1)}; + if(x>=minX && x<=maxX){ + xbins[nbins] = hist->GetXaxis()->GetBinLowEdge(i+1); + xbins[nbins+1]= hist->GetXaxis()->GetBinUpEdge(i+1); + ++nbins; + } + } + + TH1F *hnew = new TH1F ("hnew", "", nbins, xbins); + for(int i=0; iSetBinContent(i+1, hist->GetBinContent(hist->FindBin(hnew->GetBinCenter(i+1)))); + hnew->SetBinError (i+1, hist->GetBinError (hist->FindBin(hnew->GetBinCenter(i+1)))); + } + return hnew; +} + + + + diff --git a/JetAnalyzers/test/function.py b/JetAnalyzers/test/function.py new file mode 100644 index 00000000..12036e0a --- /dev/null +++ b/JetAnalyzers/test/function.py @@ -0,0 +1,70 @@ +from variable import * + +#-------------------------------------------------------------------------------- +# +# function to add root files +# +def buildConfigFile_hadd(haddCommand, shellFileName_full, inputFileNames, outputFileName_full): + + shellFile = open(shellFileName_full, "w") + shellFile.write("#!/bin/csh -f\n") + shellFile.write("\n") + # CV: delete output file in case it exists + shellFile.write("rm -f %s\n" % outputFileName_full) + shellFile.write("\n") + haddCommandLine = "%s %s" % (haddCommand, outputFileName_full) + for inputFileName in inputFileNames: + haddCommandLine += " %s" % inputFileName + shellFile.write("%s\n" % haddCommandLine) + shellFile.close() + + logFileName_full = shellFileName_full.replace('.csh', '.log') + + retVal = {} + retVal['shellFileName'] = shellFileName_full + retVal['outputFileName'] = outputFileName_full + retVal['logFileName'] = logFileName_full + + return retVal +#-------------------------------------------------------------------------------- + +#-------------------------------------------------------------------------------- +# +# function to make seperate string from a vector of string +# +def make_MakeFile_vstring(list_of_strings): + retVal = "" + for i, string_i in enumerate(list_of_strings): + if i > 0: + retVal += " " + retVal += string_i + return retVal +#-------------------------------------------------------------------------------- + +#-------------------------------------------------------------------------------- +# +# function to make seperate double from a vector of double +# +def make_MakeFile_vdouble(list_of_doubles): + retVal = "" + for i, double_i in enumerate(list_of_doubles): + if i > 0: + retVal += " " + retVal += "%2.1f" % double_i + return retVal +#-------------------------------------------------------------------------------- + +#-------------------------------------------------------------------------------- +# +# pt and eta binning +# +def make_jrAnalyzer_config(configFileName): + configFile = open(configFileName, "w") + configFile.write("drmax = 0.3\n") + configFile.write("etabarrelmin = -1.3\n") + configFile.write("etabarrelmax = 1.3\n") + configFile.write("binspt = %s\n" % make_MakeFile_vdouble(ptBinning)) + configFile.write("binseta = %s\n" % make_MakeFile_vdouble(etaBinning)) + configFile.close() +#-------------------------------------------------------------------------------- + diff --git a/JetAnalyzers/test/runJRAtauWorkflow.py b/JetAnalyzers/test/runJRAtauWorkflow.py index 9f688cb8..f72014be 100755 --- a/JetAnalyzers/test/runJRAtauWorkflow.py +++ b/JetAnalyzers/test/runJRAtauWorkflow.py @@ -1,526 +1,532 @@ #!/usr/bin/env python -#from TauAnalysis.TauIdEfficiency.tools.buildConfigFilesTauIdEffAnalysis import buildConfigFile_hadd +from variable import * +from function import * import os import re -version = 'v1_2enRecoveryCBa' -era = 'TauJec11V1' - -#inputFilePath = '/data2/veelken/CMSSW_4_2_x/JRAtauNtuples/Ztautau/v1_2enRecovery/' \ -# + 'user/v/veelken/CMSSW_4_2_x/JRAtauNtuples/Ztautau/v1_2enRecovery' -inputFilePath = '/hdfs/cms/store/user/calpas/DYJetsToLL_M-50_13TeV-madgraph-pythia8/DYJetsToLL_M-50_13TeV-madgraph-pythia8_TransferFunc_v1/a683f4f5fc0a3cbafdd5a17e17e4babe/' - -#outputFilePath = '/data1/veelken/tmp/JRAtau/%s' % version -outputFilePath = '/home/calpas/TransferFunction/CMSSW_7_2_5/src/JetMETAnalysis/JetAnalyzers/test/JRAtau/%s' % version - -samplesToAnalyze = [ - 'Ztautau' -] - -algorithms = [ - 'ak5tauHPSlooseCombDBcorrAll', - 'ak5tauHPSlooseCombDBcorrOneProng0Pi0', - 'ak5tauHPSlooseCombDBcorrOneProng1Pi0', - 'ak5tauHPSlooseCombDBcorrOneProng2Pi0', - 'ak5tauHPSlooseCombDBcorrThreeProng0Pi0', - 'ak5tauHPSmediumCombDBcorrAll', - 'ak5tauHPSmediumCombDBcorrOneProng0Pi0', - 'ak5tauHPSmediumCombDBcorrOneProng1Pi0', - 'ak5tauHPSmediumCombDBcorrOneProng2Pi0', - 'ak5tauHPSmediumCombDBcorrThreeProng0Pi0', - 'ak5tauHPStightCombDBcorrAll', - 'ak5tauHPStightCombDBcorrOneProng0Pi0', - 'ak5tauHPStightCombDBcorrOneProng1Pi0', - 'ak5tauHPStightCombDBcorrOneProng2Pi0', - 'ak5tauHPStightCombDBcorrThreeProng0Pi0', -] - -execDir = "%s/bin/%s/" % (os.environ['CMSSW_BASE'], os.environ['SCRAM_ARCH']) - -executable_jrAnalyzer = execDir + 'jet_response_analyzer_x' -executable_fitResponse = execDir + 'jet_response_fitter_x' -executable_fitResolution = execDir + 'jet_response_and_resolution_x' -executable_fitL3param = execDir + 'jet_l3_correction_x' -executable_fitL2param = execDir + 'jet_l2_correction_x' -executable_applyL2L3param = execDir + 'jet_apply_jec_x' -executable_showGraphs = execDir + 'jet_inspect_graphs_x' -executable_showHistos = execDir + 'jet_inspect_histos_x' -executable_showProfiles = execDir + 'jet_inspect_profiles_x' -executable_hadd = 'hadd -f' -executable_shell = '/bin/csh' -executable_python = 'python' - -# define function used for fitting tau-jet response: -# 0 = Gaussian -# 1 = Crystall-Ball function -fitOption = 1 - -jecTextFilePath = os.getcwd() - -if not os.path.exists(outputFilePath): - os.mkdir(outputFilePath) - -outputFilePath_plots = os.path.join(outputFilePath, "plots") -if not os.path.exists(outputFilePath_plots): - os.mkdir(outputFilePath_plots) + #-------------------------------------------------------------------------------- # -# build shell script for running 'hadd' in order to "harvest" histograms -# produced by FWLiteZllRecoilCorrectionAnalyzer macro +# create cristalBall fit directories # -fileNames_hadd = {} +outputFilePath = '%s/src/JetMETAnalysis/JetAnalyzers/test/JRAtau/%s' % (os.environ['CMSSW_BASE'], version) + +fitDirPol1 = outputFilePath+"/plots/pol1/" +fitDirPol7 = outputFilePath+"/plots/pol7/" +cuts = ["pass", "failed"] +for cut in cuts: + for alg in algorithms: + outputPlotPath = fitDirPol1+'%s/crystalBallFit/%s/' % (alg, cut) + if not os.path.exists(outputPlotPath): + os.makedirs(outputPlotPath) + outputPlotPath = fitDirPol7+'%s/crystalBallFit/%s/' % (alg, cut) + if not os.path.exists(outputPlotPath): + os.makedirs(outputPlotPath) + + alg+="l2" # l2 correction + outputPlotPath = fitDirPol1+'%s/crystalBallFit/%s/' % (alg, cut) + if not os.path.exists(outputPlotPath): + os.makedirs(outputPlotPath) + outputPlotPath = fitDirPol7+'%s/crystalBallFit/%s/' % (alg, cut) + if not os.path.exists(outputPlotPath): + os.makedirs(outputPlotPath) +#-------------------------------------------------------------------------------- -ntupleFileNames = os.listdir(inputFilePath) -#print(ntupleFileNames) +#-------------------------------------------------------------------------------- +# +# build shell script for running 'hadd' in order to "harvest" histograms +# +ntupleFileNames = [] +for inputFilePath in inputFilePaths: + ntupleFileNames.extend([ os.path.join(inputFilePath, file) for file in os.listdir(inputFilePath) ]) -ntupleFile_regex = r"ntupleJRAtau_(?P\w*)_(?P\d*)_(?P\d*)_(?P[a-zA-Z0-9]*).root" +ntupleFile_regex = r"[a-zA-Z0-9_/:.]*ntupleJRAtau_(?P[a-zA-Z0-9]*)_(?P\d*)_(?P\d*)_(?P[a-zA-Z0-9]*).root" ntupleFile_matcher = re.compile(ntupleFile_regex) +#-------------------------------------------------------------------------------- -#################### -def buildConfigFile_hadd(haddCommand, shellFileName_full, inputFileNames, outputFileName_full): - - """Build shell script to run 'hadd' command in order to add all histograms - in files specified by inputFileNames argument and write the sum to file outputFileName""" - - shellFile = open(shellFileName_full, "w") - shellFile.write("#!/bin/csh -f\n") - shellFile.write("\n") - # CV: delete output file in case it exists - shellFile.write("rm -f %s\n" % outputFileName_full) - shellFile.write("\n") - haddCommandLine = "%s %s" % (haddCommand, outputFileName_full) - for inputFileName in inputFileNames: - haddCommandLine += " %s" % inputFileName - shellFile.write("%s\n" % haddCommandLine) - shellFile.close() - - logFileName_full = shellFileName_full.replace('.csh', '.log') - - retVal = {} - retVal['shellFileName'] = shellFileName_full - retVal['outputFileName'] = outputFileName_full - retVal['logFileName'] = logFileName_full - - return retVal -#################### - - +#-------------------------------------------------------------------------------- +# +# add all samples files +# +haddInputFileNames = [] for sampleToAnalyze in samplesToAnalyze: - haddInputFileNames = [] + nbOfSample = [] for ntupleFileName in ntupleFileNames: if ntupleFile_matcher.match(ntupleFileName) and \ - ntupleFile_matcher.match(ntupleFileName).group('sample') == sampleToAnalyze: - haddInputFileNames.append(os.path.join(inputFilePath, ntupleFileName)) - - print "sample = %s: found %i input files." % (sampleToAnalyze, len(haddInputFileNames)) + ntupleFile_matcher.match(ntupleFileName).group('sample') == sampleToAnalyze: + haddInputFileNames.append(ntupleFileName) + nbOfSample.append(ntupleFileName) + print "sample = %s: found %i input files." % (sampleToAnalyze, len(nbOfSample)) +print "found %i input files." % (len(haddInputFileNames)) - haddShellFileName = os.path.join(outputFilePath, 'harvestJRAtauNtuples_%s.csh' % sampleToAnalyze) - haddOutputFileName = os.path.join(inputFilePath, 'ntupleJRAtau_%s_all.root' % sampleToAnalyze) +haddShellFileName = os.path.join(outputFilePath, 'harvestJRAtauNtuples.csh') +haddOutputFileName = os.path.join(outputFilePath, 'ntupleJRAtau_all.root') # copy where you can +retVal_hadd = buildConfigFile_hadd(executable_hadd, haddShellFileName, haddInputFileNames, haddOutputFileName) - retVal_hadd = \ - buildConfigFile_hadd(executable_hadd, haddShellFileName, haddInputFileNames, haddOutputFileName) - - fileNames_hadd[sampleToAnalyze] = {} - fileNames_hadd[sampleToAnalyze]['shellFileName'] = haddShellFileName - fileNames_hadd[sampleToAnalyze]['inputFileNames'] = haddInputFileNames - fileNames_hadd[sampleToAnalyze]['outputFileName'] = haddOutputFileName - fileNames_hadd[sampleToAnalyze]['logFileName'] = retVal_hadd['logFileName'] +fileNames_hadd = {} +fileNames_hadd['shellFileName'] = haddShellFileName +fileNames_hadd['inputFileNames'] = haddInputFileNames +fileNames_hadd['outputFileName'] = haddOutputFileName +fileNames_hadd['logFileName'] = retVal_hadd['logFileName'] #-------------------------------------------------------------------------------- -def make_MakeFile_vstring(list_of_strings): - retVal = "" - for i, string_i in enumerate(list_of_strings): - if i > 0: - retVal += " " - retVal += string_i - return retVal +#-------------------------------------------------------------------------------- +# +# jet response for uncalibrated tau-jets +# +fileNames_and_options_jrAnalyzer = {} +fileNames_and_options_jrAnalyzer['inputFileNames'] = [fileNames_hadd['outputFileName']] +fileNames_and_options_jrAnalyzer['outputFileName'] = os.path.join(outputFilePath, "jet_response_analyzer_uncalibrated.root") +fileNames_and_options_jrAnalyzer['configFileName'] = os.path.join(outputFilePath, "jet_response_analyzer_uncalibrated.cfg") +make_jrAnalyzer_config(fileNames_and_options_jrAnalyzer['configFileName']) +fileNames_and_options_jrAnalyzer['logFileName'] = os.path.join(outputFilePath, "jet_response_analyzer_uncalibrated.log") +fileNames_and_options_jrAnalyzer['commandLine'] = '%s -input %s -output %s -algs %s -nbinsrelrsp 250' % \ + (fileNames_and_options_jrAnalyzer['configFileName'], + make_MakeFile_vstring(fileNames_and_options_jrAnalyzer['inputFileNames']), + fileNames_and_options_jrAnalyzer['outputFileName'], + "".join([ "%s:0.3 " % algorithm for algorithm in algorithms ]) + ) +#-------------------------------------------------------------------------------- #-------------------------------------------------------------------------------- # -# initialize command-line parameters for analyzing "plain" ROOT Ntuples for uncalibrated tau-jets -# and filling tau-jet response and resolution histograms +# round 1.1: fit jet response for uncalibrated tau-jets # -fileNames_and_options_jrAnalyzer = {} +fileNames_and_options_fitResponse_uncalibrated = {} +haddFitHists = [] +for histName in histNames: + fileNames_and_options_fitResponse_uncalibrated[histName] = {} + fileNames_and_options_fitResponse_uncalibrated[histName]['inputFileNames'] = [fileNames_and_options_jrAnalyzer['outputFileName']] + fileNames_and_options_fitResponse_uncalibrated[histName]['outputFileName'] = os.path.join(outputFilePath, "responseJRAtau_%s_uncalibrated.root" % (histName)) + haddFitHists.append(fileNames_and_options_fitResponse_uncalibrated[histName]['outputFileName']) + fileNames_and_options_fitResponse_uncalibrated[histName]['logFileName'] = os.path.join(outputFilePath, "responseJRAtau_%s_uncalibrated.log" %histName) + fileNames_and_options_fitResponse_uncalibrated[histName]['commandLine'] = '-input %s -output %s -algs %s -fittype %i -histName %s -isItCalibrated uncalibrated -pol1 false -normalized true -fitDirPol7 %s' % \ + (make_MakeFile_vstring(fileNames_and_options_fitResponse_uncalibrated[histName]['inputFileNames']), + fileNames_and_options_fitResponse_uncalibrated[histName]['outputFileName'], + make_MakeFile_vstring(algorithms), + fitOption, + histName, + fitDirPol7 + ) +#-------------------------------------------------------------------------------- -ptBinning = [ - 20., 22.5, 25., 27.5, 30., 35., 40., 45., 50., 60., 80., 120., 200. -] +#-------------------------------------------------------------------------------- +# +# round 1.2: add the uncalibrated fit histograms +# +haddShellFitHistName = os.path.join(outputFilePath, 'harvestFitHist_uncalibrated.csh') +haddOutputFitHistName = os.path.join(outputFilePath, 'fitHists_uncalibrated.root') -etaBinning = [ - -2.5, -2.3, -2.1, -1.9, -1.7, -1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, - +0.1, +0.3, +0.5, +0.7, +0.9, +1.1, +1.3, +1.5, +1.7, +1.9, +2.1, +2.3, +2.5 -] - -def make_MakeFile_vdouble(list_of_doubles): - retVal = "" - for i, double_i in enumerate(list_of_doubles): - if i > 0: - retVal += " " - retVal += "%2.1f" % double_i - return retVal +retFitHist_hadd = buildConfigFile_hadd(executable_haddFitHists, haddShellFitHistName, haddFitHists, haddOutputFitHistName) -def make_jrAnalyzer_config(configFileName): - configFile = open(configFileName, "w") - configFile.write("drmax = 0.3\n") - configFile.write("etabarrelmin = -1.3\n") - configFile.write("etabarrelmax = 1.3\n") - configFile.write("binspt = %s\n" % make_MakeFile_vdouble(ptBinning)) - configFile.write("binseta = %s\n" % make_MakeFile_vdouble(etaBinning)) - configFile.close() +fitHists_hadd_uncalibrated = {} +fitHists_hadd_uncalibrated['shellFileName'] = haddShellFitHistName +fitHists_hadd_uncalibrated['inputFileNames'] = haddFitHists +fitHists_hadd_uncalibrated['outputFileName'] = haddOutputFitHistName +fitHists_hadd_uncalibrated['logFileName'] = retFitHist_hadd['logFileName'] +#-------------------------------------------------------------------------------- -for sampleToAnalyze in samplesToAnalyze: - fileNames_and_options_jrAnalyzer[sampleToAnalyze] = {} - fileNames_and_options_jrAnalyzer[sampleToAnalyze]['inputFileNames'] = \ - [ fileNames_hadd[sampleToAnalyze]['outputFileName'] ] - fileNames_and_options_jrAnalyzer[sampleToAnalyze]['outputFileName'] = \ - os.path.join(outputFilePath, "histogramsJRAtau_%s.root" % sampleToAnalyze) - fileNames_and_options_jrAnalyzer[sampleToAnalyze]['configFileName'] = \ - os.path.join(outputFilePath, "jet_response_analyzer_%s.cfg" % sampleToAnalyze) - make_jrAnalyzer_config(fileNames_and_options_jrAnalyzer[sampleToAnalyze]['configFileName']) - fileNames_and_options_jrAnalyzer[sampleToAnalyze]['logFileName'] = \ - os.path.join(outputFilePath, "jet_response_analyzer_%s.log" % sampleToAnalyze) - fileNames_and_options_jrAnalyzer[sampleToAnalyze]['commandLine'] = \ - '%s -input %s -output %s -algs %s' % \ - (fileNames_and_options_jrAnalyzer[sampleToAnalyze]['configFileName'], - make_MakeFile_vstring(fileNames_and_options_jrAnalyzer[sampleToAnalyze]['inputFileNames']), - fileNames_and_options_jrAnalyzer[sampleToAnalyze]['outputFileName'], - "".join([ "%s:0.3 " % algorithm for algorithm in algorithms ])) +#-------------------------------------------------------------------------------- +# +# round 2.1: fit jet response for uncalibrated tau-jets +# Fixed the CB parameters determined above and use pol1 instead of pol7 +# +fileNames_and_options_fitResponse_uncalibrated_pol1 = {} +haddFitHists_pol1 = [] +for histName in histNames: + fileNames_and_options_fitResponse_uncalibrated_pol1[histName] = {} + fileNames_and_options_fitResponse_uncalibrated_pol1[histName]['inputFileNames'] = [fitHists_hadd_uncalibrated['outputFileName']] + fileNames_and_options_fitResponse_uncalibrated_pol1[histName]['outputFileName'] = os.path.join(outputFilePath, "responseJRAtau_%s_uncalibrated_pol1.root" % (histName)) + haddFitHists_pol1.append(fileNames_and_options_fitResponse_uncalibrated_pol1[histName]['outputFileName']) + fileNames_and_options_fitResponse_uncalibrated_pol1[histName]['logFileName'] = os.path.join(outputFilePath, "responseJRAtau_%s_uncalibrated_plo1.log" %histName) + fileNames_and_options_fitResponse_uncalibrated_pol1[histName]['commandLine'] = '-input %s -output %s -algs %s -fittype %i -histName %s -isItCalibrated uncalibrated -pol1 true -normalized true -fitDirPol1 %s' % \ + (make_MakeFile_vstring(fileNames_and_options_fitResponse_uncalibrated_pol1[histName]['inputFileNames']), + fileNames_and_options_fitResponse_uncalibrated_pol1[histName]['outputFileName'], + make_MakeFile_vstring(algorithms), + fitOption, + histName, + fitDirPol1 + ) #-------------------------------------------------------------------------------- #-------------------------------------------------------------------------------- # -# initialize command-line parameters for fitting jet response for uncalibrated tau-jets +# round 2.2: add the uncalibrated fit histograms # -fileNames_and_options_fitResponse_uncalibrated = {} +haddShellFitHistName_pol1 = os.path.join(outputFilePath, 'harvestFitHist_uncalibrated_pol1.csh') +haddOutputFitHistName_pol1 = os.path.join(outputFilePath, 'fitHists_uncalibrated_pol1.root') -for sampleToAnalyze in samplesToAnalyze: - fileNames_and_options_fitResponse_uncalibrated[sampleToAnalyze] = {} - fileNames_and_options_fitResponse_uncalibrated[sampleToAnalyze]['inputFileNames'] = \ - [ fileNames_and_options_jrAnalyzer[sampleToAnalyze]['outputFileName'] ] - fileNames_and_options_fitResponse_uncalibrated[sampleToAnalyze]['outputFileName'] = \ - os.path.join(outputFilePath, "responseJRAtau_%s.root" % sampleToAnalyze) - fileNames_and_options_fitResponse_uncalibrated[sampleToAnalyze]['logFileName'] = \ - os.path.join(outputFilePath, "jet_response_fitter_%s.log" % sampleToAnalyze) - fileNames_and_options_fitResponse_uncalibrated[sampleToAnalyze]['commandLine'] = \ - '-input %s -output %s -algs %s -fittype %i' % \ - (make_MakeFile_vstring(fileNames_and_options_fitResponse_uncalibrated[sampleToAnalyze]['inputFileNames']), - fileNames_and_options_fitResponse_uncalibrated[sampleToAnalyze]['outputFileName'], - make_MakeFile_vstring(algorithms), - fitOption) +retFitHist_hadd_pol1 = buildConfigFile_hadd(executable_haddFitHists, haddShellFitHistName_pol1, haddFitHists_pol1, haddOutputFitHistName_pol1) + +fitHists_hadd_uncalibrated_pol1 = {} +fitHists_hadd_uncalibrated_pol1['shellFileName'] = haddShellFitHistName_pol1 +fitHists_hadd_uncalibrated_pol1['inputFileNames'] = haddFitHists_pol1 +fitHists_hadd_uncalibrated_pol1['outputFileName'] = haddOutputFitHistName_pol1 +fitHists_hadd_uncalibrated_pol1['logFileName'] = retFitHist_hadd_pol1['logFileName'] #-------------------------------------------------------------------------------- #-------------------------------------------------------------------------------- # -# initialize command-line parameters for determining L2 and L3 correction parameters +# determining L2L3 correction and fit them # -fileNames_and_options_fitL3param = {} fileNames_and_options_fitL2param = {} +fileNames_and_options_fitL2param['inputFileNames'] = [fitHists_hadd_uncalibrated_pol1['outputFileName']] +fileNames_and_options_fitL2param['outputFileName'] = os.path.join(outputFilePath, "fitL2param.root") +fileNames_and_options_fitL2param['logFileName'] = os.path.join(outputFilePath, "fitL2param.log") +fileNames_and_options_fitL2param['commandLine'] = \ + '-input %s -output %s -era %s -algs %s -formats png -batch true -l2l3 true -fitDirPol7 %s' % \ + (make_MakeFile_vstring(fileNames_and_options_fitL2param['inputFileNames']), + fileNames_and_options_fitL2param['outputFileName'], + era, + make_MakeFile_vstring(algorithms), + fitDirPol7 + ) +#-------------------------------------------------------------------------------- -for sampleToAnalyze in samplesToAnalyze: - fileNames_and_options_fitL3param[sampleToAnalyze] = {} - fileNames_and_options_fitL3param[sampleToAnalyze]['inputFileNames'] = \ - [ fileNames_and_options_fitResponse_uncalibrated[sampleToAnalyze]['outputFileName'] ] - fileNames_and_options_fitL3param[sampleToAnalyze]['outputFileName'] = \ - os.path.join(outputFilePath, "fitL3param_%s.root" % sampleToAnalyze) - fileNames_and_options_fitL3param[sampleToAnalyze]['logFileName'] = \ - os.path.join(outputFilePath, "fitL3param_%s.log" % sampleToAnalyze) - fileNames_and_options_fitL3param[sampleToAnalyze]['commandLine'] = \ - '-input %s -output %s -era %s -algs %s -formats png -batch true' % \ - (make_MakeFile_vstring(fileNames_and_options_fitL3param[sampleToAnalyze]['inputFileNames']), - fileNames_and_options_fitL3param[sampleToAnalyze]['outputFileName'], - era, - make_MakeFile_vstring(algorithms)) - - fileNames_and_options_fitL2param[sampleToAnalyze] = {} - fileNames_and_options_fitL2param[sampleToAnalyze]['inputFileNames'] = \ - [ fileNames_and_options_fitResponse_uncalibrated[sampleToAnalyze]['outputFileName'], - fileNames_and_options_fitL3param[sampleToAnalyze]['outputFileName'] ] - fileNames_and_options_fitL2param[sampleToAnalyze]['outputFileName'] = \ - os.path.join(outputFilePath, "fitL2param_%s.root" % sampleToAnalyze) - fileNames_and_options_fitL2param[sampleToAnalyze]['logFileName'] = \ - os.path.join(outputFilePath, "fitL2param_%s.log" % sampleToAnalyze) - fileNames_and_options_fitL2param[sampleToAnalyze]['commandLine'] = \ - '-input %s -l3input %s -output %s -era %s -algs %s -formats png -batch true' % \ - (fileNames_and_options_fitResponse_uncalibrated[sampleToAnalyze]['outputFileName'], - fileNames_and_options_fitL3param[sampleToAnalyze]['outputFileName'], - fileNames_and_options_fitL2param[sampleToAnalyze]['outputFileName'], - era, - make_MakeFile_vstring(algorithms)) -#-------------------------------------------------------------------------------- - -#-------------------------------------------------------------------------------- -# -# initialize command-line parameters for applying L2/L3 correction parameters -# and producing new "plain" ROOT Ntuples for calibrated tau-jets +#-------------------------------------------------------------------------------- +# +# apply L2/L3 correction (determined above) to the original tree # fileNames_and_options_applyL2L3param = {} - -for sampleToAnalyze in samplesToAnalyze: - fileNames_and_options_applyL2L3param[sampleToAnalyze] = {} - fileNames_and_options_applyL2L3param[sampleToAnalyze]['inputFileNames'] = \ - [ fileNames_and_options_jrAnalyzer[sampleToAnalyze]['outputFileName'], - fileNames_and_options_fitL3param[sampleToAnalyze]['outputFileName'], - fileNames_and_options_fitL2param[sampleToAnalyze]['outputFileName'] ] - fileNames_and_options_applyL2L3param[sampleToAnalyze]['outputFileName'] = \ - os.path.join(outputFilePath, "applyL2L3param_%s.root" % sampleToAnalyze) - fileNames_and_options_applyL2L3param[sampleToAnalyze]['logFileName'] = \ - os.path.join(outputFilePath, "applyL2L3param_%s.log" % sampleToAnalyze) - fileNames_and_options_applyL2L3param[sampleToAnalyze]['commandLine'] = \ - '-input %s -output %s -era %s -algs %s -jecpath %s -levels 2 3' % \ - (fileNames_hadd[sampleToAnalyze]['outputFileName'], - fileNames_and_options_applyL2L3param[sampleToAnalyze]['outputFileName'], - era, - make_MakeFile_vstring(algorithms), - jecTextFilePath) +fileNames_and_options_applyL2L3param['inputFileNames'] = [ fileNames_and_options_jrAnalyzer['outputFileName'], + fileNames_and_options_fitL2param['outputFileName'] ] +fileNames_and_options_applyL2L3param['outputFileName'] = os.path.join(outputFilePath, "applyL2L3param.root") +fileNames_and_options_applyL2L3param['logFileName'] = os.path.join(outputFilePath, "applyL2L3param.log") +fileNames_and_options_applyL2L3param['commandLine'] = '-input %s -output %s -era %s -algs %s -jecpath %s -levels 2 -saveitree false' % \ + (fileNames_hadd['outputFileName'], + fileNames_and_options_applyL2L3param['outputFileName'], + era, + make_MakeFile_vstring(algorithms), + jecTextFilePath) #-------------------------------------------------------------------------------- #-------------------------------------------------------------------------------- # -# initialize command-line parameters for analyzing "plain" ROOT Ntuples for calibrated tau-jets -# and filling tau-jet response and resolution histograms +# jet response for calibrated tau-jets # -fileNames_and_options_jrAnalyzer_calibrated = {} +fileNames_and_options_jrAnalyzer_calibrated = {} +fileNames_and_options_jrAnalyzer_calibrated['inputFileNames'] = [fileNames_and_options_applyL2L3param['outputFileName']] +fileNames_and_options_jrAnalyzer_calibrated['outputFileName'] = os.path.join(outputFilePath, "jet_response_analyzer_calibrated.root") +fileNames_and_options_jrAnalyzer_calibrated['configFileName'] = os.path.join(outputFilePath, "jet_response_analyzer_calibrated.cfg") +make_jrAnalyzer_config(fileNames_and_options_jrAnalyzer_calibrated['configFileName']) +fileNames_and_options_jrAnalyzer_calibrated['logFileName'] = os.path.join(outputFilePath, "jet_response_analyzer_calibrated.log") +fileNames_and_options_jrAnalyzer_calibrated['commandLine'] = '%s -input %s -output %s -algs %s -nbinsrelrsp 250' % \ + (fileNames_and_options_jrAnalyzer_calibrated['configFileName'], + make_MakeFile_vstring(fileNames_and_options_jrAnalyzer_calibrated['inputFileNames']), + fileNames_and_options_jrAnalyzer_calibrated['outputFileName'], + make_MakeFile_vstring([ "".join([ algorithm, suffix]) for algorithm in algorithms for suffix in ["l2"] ])) +#-------------------------------------------------------------------------------- -for sampleToAnalyze in samplesToAnalyze: - fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze] = {} - fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze]['inputFileNames'] = \ - [ fileNames_and_options_applyL2L3param[sampleToAnalyze]['outputFileName'] ] - fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze]['outputFileName'] = \ - os.path.join(outputFilePath, "histogramsJRAtau_%s_calibrated.root" % sampleToAnalyze) - fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze]['configFileName'] = \ - os.path.join(outputFilePath, "jet_response_analyzer_%s_calibrated.cfg" % sampleToAnalyze) - make_jrAnalyzer_config(fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze]['configFileName']) - fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze]['logFileName'] = \ - os.path.join(outputFilePath, "jet_response_analyzer_%s_calibrated.log" % sampleToAnalyze) - fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze]['commandLine'] = \ - '%s -input %s -output %s -algs %s' % \ - (fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze]['configFileName'], - make_MakeFile_vstring(fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze]['inputFileNames']), - fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze]['outputFileName'], - make_MakeFile_vstring([ "".join([ algorithm, suffix]) for algorithm in algorithms for suffix in [ "", "l2l3"] ])) + +#-------------------------------------------------------------------------------- +# +# fit jet response for calibrated tau-jets +# +fileNames_and_options_fitResponse_calibrated = {} +haddFitHists_calibrated = [] +for histName in histNames: + fileNames_and_options_fitResponse_calibrated[histName] = {} + fileNames_and_options_fitResponse_calibrated[histName]['inputFileNames'] = [fileNames_and_options_jrAnalyzer_calibrated['outputFileName']] + fileNames_and_options_fitResponse_calibrated[histName]['outputFileName'] = os.path.join(outputFilePath, "responseJRAtau_%s_calibrated.root" % (histName)) + haddFitHists_calibrated.append(fileNames_and_options_fitResponse_calibrated[histName]['outputFileName']) + fileNames_and_options_fitResponse_calibrated[histName]['logFileName'] = os.path.join(outputFilePath, "responseJRAtau_%s_calibrated.log" %histName) + fileNames_and_options_fitResponse_calibrated[histName]['commandLine'] = '-input %s -output %s -algs %s -fittype %i -histName %s -isItCalibrated calibrated -pol1 false -normalized true -fitDirPol7 %s' % \ + (make_MakeFile_vstring(fileNames_and_options_fitResponse_calibrated[histName]['inputFileNames']), + fileNames_and_options_fitResponse_calibrated[histName]['outputFileName'], + #make_MakeFile_vstring([ "".join([ algorithm, suffix]) for algorithm in algorithms for suffix in [ "", "l2"] ]), + make_MakeFile_vstring([ "".join([ algorithm, suffix]) for algorithm in algorithms for suffix in [ "l2"] ]), # now fit only the calibrate + fitOption, + histName, + fitDirPol7 + ) +#-------------------------------------------------------------------------------- + +#-------------------------------------------------------------------------------- +# +# add the calibrated fit histograms +# +haddShellFitHistName_calibrated = os.path.join(outputFilePath, 'harvestFitHist_calibrated.csh') +haddOutputFitHistName_calibrated = os.path.join(outputFilePath, 'fitHists_calibrated.root') + +retFitHist_hadd_calibrated = buildConfigFile_hadd(executable_haddFitHists, haddShellFitHistName_calibrated, haddFitHists_calibrated, haddOutputFitHistName_calibrated) + +fitHists_hadd_calibrated = {} +fitHists_hadd_calibrated['shellFileName'] = haddShellFitHistName_calibrated +fitHists_hadd_calibrated['inputFileNames'] = haddFitHists_calibrated +fitHists_hadd_calibrated['outputFileName'] = haddOutputFitHistName_calibrated +fitHists_hadd_calibrated['logFileName'] = retFitHist_hadd_calibrated['logFileName'] #-------------------------------------------------------------------------------- #-------------------------------------------------------------------------------- # -# initialize command-line parameters for fitting jet response and resolution -# for calibrated as well as uncalibrated tau-jets +# add the uncalibrated and calibrated fit histograms for the next step # -fileNames_and_options_fitResponse_calibrated = {} +haddShellFitHistName_all = os.path.join(outputFilePath, 'harvestFitHist_all.csh') +haddOutputFitHistName_all = os.path.join(outputFilePath, 'fitHists_all.root') +haddFitHists_all = [ fitHists_hadd_uncalibrated['outputFileName'], fitHists_hadd_calibrated['outputFileName']] +retFitHist_hadd_all = buildConfigFile_hadd(executable_haddFitHists, haddShellFitHistName_all, haddFitHists_all, haddOutputFitHistName_all) +fitHists_hadd_all = {} +fitHists_hadd_all['shellFileName'] = haddShellFitHistName_all +fitHists_hadd_all['inputFileNames'] = haddFitHists_all +fitHists_hadd_all['outputFileName'] = haddOutputFitHistName_all +fitHists_hadd_all['logFileName'] = retFitHist_hadd_all['logFileName'] +#-------------------------------------------------------------------------------- + +#-------------------------------------------------------------------------------- +# +# jet response and resolution graphs fileNames_and_options_fitResolution = {} +fileNames_and_options_fitResolution['inputFileNames'] = [fitHists_hadd_all['outputFileName']] +fileNames_and_options_fitResolution['outputFileName'] = os.path.join(outputFilePath, "resolutionJRAtau.root") +fileNames_and_options_fitResolution['logFileName'] = os.path.join(outputFilePath, "resolutionJRAtau.log") +fileNames_and_options_fitResolution['commandLine'] = '-input %s -output %s -algs %s -dorelrsp true' % \ + (make_MakeFile_vstring(fileNames_and_options_fitResolution['inputFileNames']), + fileNames_and_options_fitResolution['outputFileName'], + make_MakeFile_vstring([ "".join([ algorithm, suffix]) for algorithm in algorithms for suffix in [ "", "l2"] ])) +#-------------------------------------------------------------------------------- -for sampleToAnalyze in samplesToAnalyze: - fileNames_and_options_fitResponse_calibrated[sampleToAnalyze] = {} - fileNames_and_options_fitResponse_calibrated[sampleToAnalyze]['inputFileNames'] = \ - [ fileNames_and_options_jrAnalyzer_calibrated[sampleToAnalyze]['outputFileName'] ] - fileNames_and_options_fitResponse_calibrated[sampleToAnalyze]['outputFileName'] = \ - os.path.join(outputFilePath, "responseJRAtau_%s_calibrated.root" % sampleToAnalyze) - fileNames_and_options_fitResponse_calibrated[sampleToAnalyze]['logFileName'] = \ - os.path.join(outputFilePath, "jet_response_fitter_%s_calibrated.log" % sampleToAnalyze) - fileNames_and_options_fitResponse_calibrated[sampleToAnalyze]['commandLine'] = \ - '-input %s -output %s -algs %s -fittype %i' % \ - (make_MakeFile_vstring(fileNames_and_options_fitResponse_calibrated[sampleToAnalyze]['inputFileNames']), - fileNames_and_options_fitResponse_calibrated[sampleToAnalyze]['outputFileName'], - make_MakeFile_vstring([ "".join([ algorithm, suffix]) for algorithm in algorithms for suffix in [ "", "l2l3"] ]), - fitOption) - - fileNames_and_options_fitResolution[sampleToAnalyze] = {} - fileNames_and_options_fitResolution[sampleToAnalyze]['inputFileNames'] = \ - [ fileNames_and_options_fitResponse_calibrated[sampleToAnalyze]['outputFileName'] ] - fileNames_and_options_fitResolution[sampleToAnalyze]['outputFileName'] = \ - os.path.join(outputFilePath, "resolutionJRAtau_%s.root" % sampleToAnalyze) - fileNames_and_options_fitResolution[sampleToAnalyze]['logFileName'] = \ - os.path.join(outputFilePath, "jet_resolution_fitter_%s.log" % sampleToAnalyze) - fileNames_and_options_fitResolution[sampleToAnalyze]['commandLine'] = \ - '-input %s -output %s -algs %s -dorelrsp true -doetarsp true -docbfits true' % \ - (make_MakeFile_vstring(fileNames_and_options_fitResolution[sampleToAnalyze]['inputFileNames']), - fileNames_and_options_fitResolution[sampleToAnalyze]['outputFileName'], - make_MakeFile_vstring([ "".join([ algorithm, suffix]) for algorithm in algorithms for suffix in [ "", "l2l3"] ])) -#-------------------------------------------------------------------------------- - -#-------------------------------------------------------------------------------- -# -# initialize command-line parameters for making jet response plots for calibrated tau-jets -# plus resolution plots for calibrated compared to uncalibrated tau-jets + +#-------------------------------------------------------------------------------- +# +# hist: response calibrated vs uncalibrated +# graph: resolution calibrated vs uncalibrated # fileNames_and_options_showHistos = {} fileNames_and_options_showGraphs = {} -for sampleToAnalyze in samplesToAnalyze: - - fileNames_and_options_showHistos[sampleToAnalyze] = {} - fileNames_and_options_showGraphs[sampleToAnalyze] = {} +for algorithm in algorithms: + outputFilePath_plots_algorithm = outputFilePath+"/plots/%s/tauRelRsp/" %(algorithm) + if not os.path.exists(outputFilePath_plots_algorithm): + os.mkdir(outputFilePath_plots_algorithm) + + fileNames_and_options_showHistos[algorithm] = {} + fileNames_and_options_showHistos[algorithm]['inputFileNames'] = [ fitHists_hadd_all['outputFileName'] ] + fileNames_and_options_showHistos[algorithm]['outputFileName'] = "make_jet_inspect_histos_target_%s_calibrated" % (algorithm) + fileNames_and_options_showHistos[algorithm]['logFileName'] = \ + os.path.join(outputFilePath, "jet_inspect_histos_%s_calibrated.log" % (algorithm)) + fileNames_and_options_showHistos[algorithm]['commandLine'] = \ + '-inputs %s -algs %s -variables RelRsp:JetEta:RefPt %s -formats png -batch true -opath %s -colors 1 4 -peak false' % \ + (make_MakeFile_vstring(fileNames_and_options_showHistos[algorithm]['inputFileNames']), + make_MakeFile_vstring([ "".join([ algorithm, suffix]) for suffix in [ "", "l2"] ]), + "-norm true -npercanvas 1", + outputFilePath_plots_algorithm) + + + fileNames_and_options_showGraphs[algorithm] = {} + for refVariable in [ "RefPt", "JetEta" ]: + #for refVariable in [ "JetEta:RefPt" ]: + fileNames_and_options_showGraphs[algorithm][refVariable] = {} + fileNames_and_options_showGraphs[algorithm][refVariable]['inputFileNames'] = [ fileNames_and_options_fitResolution['outputFileName'] ] + fileNames_and_options_showGraphs[algorithm][refVariable]['outputFileName'] = "make_jet_inspect_graphs_%s_%s_target" % (algorithm, refVariable) + fileNames_and_options_showGraphs[algorithm][refVariable]['logFileName'] = \ + os.path.join(outputFilePath, "jet_inspect_graphs_%s_%s.log" % (algorithm, refVariable)) + if (refVariable == "RefPt"): + fileNames_and_options_showGraphs[algorithm][refVariable]['commandLine'] = \ + '-inputs %s -algs %s -variables RelRspVs%s %s -formats png -batch true -opath %s -colors 1 4 -logx true' % \ + (make_MakeFile_vstring(fileNames_and_options_showGraphs[algorithm][refVariable]['inputFileNames']), + make_MakeFile_vstring([ "".join([ algorithm, suffix]) for suffix in [ "", "l2"] ]), + refVariable, + "-ymin 0.5 -ymax 1.5 -legx 0.20 -legy 0.35 -legw 0.60", + outputFilePath_plots_algorithm) + else: + fileNames_and_options_showGraphs[algorithm][refVariable]['commandLine'] = \ + '-inputs %s -algs %s -variables RelRspVs%s %s -formats png -batch true -opath %s -colors 1 4' % \ + (make_MakeFile_vstring(fileNames_and_options_showGraphs[algorithm][refVariable]['inputFileNames']), + make_MakeFile_vstring([ "".join([ algorithm, suffix]) for suffix in [ "", "l2"] ]), + refVariable, + "-ymin 0.5 -ymax 1.5 -legx 0.20 -legy 0.35 -legw 0.60", + outputFilePath_plots_algorithm) +#-------------------------------------------------------------------------------- - for algorithm in algorithms: - outputFilePath_plots_algorithm = os.path.join(outputFilePath_plots, algorithm) - if not os.path.exists(outputFilePath_plots_algorithm): - os.mkdir(outputFilePath_plots_algorithm) - - fileNames_and_options_showHistos[sampleToAnalyze][algorithm] = {} - for refVariable in [ "RefPt", "JetEta" ]: - fileNames_and_options_showHistos[sampleToAnalyze][algorithm][refVariable] = {} - fileNames_and_options_showHistos[sampleToAnalyze][algorithm][refVariable]['inputFileNames'] = \ - [ fileNames_and_options_fitResponse_calibrated[sampleToAnalyze]['outputFileName'] ] - fileNames_and_options_showHistos[sampleToAnalyze][algorithm][refVariable]['outputFileName'] = \ - "make_jet_inspect_histos_target_%s_%s_%s_calibrated" % (sampleToAnalyze, algorithm, refVariable) - fileNames_and_options_showHistos[sampleToAnalyze][algorithm][refVariable]['logFileName'] = \ - os.path.join(outputFilePath, "jet_inspect_histos_%s_%s_%s_calibrated.log" % (sampleToAnalyze, algorithm, refVariable)) - fileNames_and_options_showHistos[sampleToAnalyze][algorithm][refVariable]['commandLine'] = \ - '-inputs %s -algs %s -variables RelRsp:%s %s -formats png -batch true -opath %s' % \ - (make_MakeFile_vstring(fileNames_and_options_showHistos[sampleToAnalyze][algorithm][refVariable]['inputFileNames']), - make_MakeFile_vstring([ "".join([ algorithm, suffix]) for suffix in [ "", "l2l3"] ]), - refVariable, - "-norm true -npercanvas 1", - outputFilePath_plots_algorithm) - - fileNames_and_options_showGraphs[sampleToAnalyze][algorithm] = {} - for refVariable in [ "RefPt", "JetEta" ]: - fileNames_and_options_showGraphs[sampleToAnalyze][algorithm][refVariable] = {} - fileNames_and_options_showGraphs[sampleToAnalyze][algorithm][refVariable]['inputFileNames'] = \ - [ fileNames_and_options_fitResolution[sampleToAnalyze]['outputFileName'] ] - fileNames_and_options_showGraphs[sampleToAnalyze][algorithm][refVariable]['outputFileName'] = \ - "make_jet_inspect_graphs_%s_%s_%s_target" % (sampleToAnalyze, algorithm, refVariable) - fileNames_and_options_showGraphs[sampleToAnalyze][algorithm][refVariable]['logFileName'] = \ - os.path.join(outputFilePath, "jet_inspect_graphs_%s_%s_%s.log" % (sampleToAnalyze, algorithm, refVariable)) - fileNames_and_options_showGraphs[sampleToAnalyze][algorithm][refVariable]['commandLine'] = \ - '-inputs %s -algs %s -variables RelRspVs%s %s -formats png -batch true -opath %s' % \ - (make_MakeFile_vstring(fileNames_and_options_showGraphs[sampleToAnalyze][algorithm][refVariable]['inputFileNames']), - make_MakeFile_vstring([ "".join([ algorithm, suffix]) for suffix in [ "", "l2l3"] ]), - refVariable, - "-ymin 0.5 -ymax 1.5 -legx 0.20 -legy 0.35 -legw 0.60", - outputFilePath_plots_algorithm) -#-------------------------------------------------------------------------------- - -# done building config files and initializing command-line parameters, now build Makefile... -makeFileName = "Makefile_runJRAtauworkflow_%s" % version +#-------------------------------------------------------------------------------- +# +# build Makefile... +# +makeFileName = "Makefile_runJRAtauworkflow_l2_%s" % version makeFile = open(makeFileName, "w") makeFile.write("\n") outputFileNames_runJRAtauworkflow = [] -for sampleName in samplesToAnalyze: +outputFileNames_runJRAtauworkflow.extend([ + fileNames_hadd['outputFileName'], + fileNames_and_options_jrAnalyzer['outputFileName'], + fitHists_hadd_uncalibrated['outputFileName'], + fileNames_and_options_fitL2param['outputFileName'], + fileNames_and_options_applyL2L3param['outputFileName'], + fileNames_and_options_jrAnalyzer_calibrated['outputFileName'], + fitHists_hadd_calibrated['outputFileName'], + fitHists_hadd_all['outputFileName'], + fileNames_and_options_fitResolution['outputFileName'] +]) + +for histName in histNames: + outputFileNames_runJRAtauworkflow.extend([ + fileNames_and_options_fitResponse_uncalibrated[histName]['outputFileName'], + fileNames_and_options_fitResponse_calibrated[histName]['outputFileName'] + ]) + +for algorithm in algorithms: outputFileNames_runJRAtauworkflow.extend([ - fileNames_and_options_jrAnalyzer[sampleName]['outputFileName'], - fileNames_and_options_fitResponse_uncalibrated[sampleName]['outputFileName'], - fileNames_and_options_fitL3param[sampleName]['outputFileName'], - fileNames_and_options_fitL2param[sampleName]['outputFileName'], - fileNames_and_options_applyL2L3param[sampleName]['outputFileName'], - fileNames_and_options_jrAnalyzer_calibrated[sampleName]['outputFileName'], - fileNames_and_options_fitResponse_calibrated[sampleName]['outputFileName'], - fileNames_and_options_fitResolution[sampleName]['outputFileName'] + fileNames_and_options_showHistos[algorithm]['outputFileName'], ]) - for algorithm in algorithms: - for refVariable in [ "RefPt", "JetEta" ]: - outputFileNames_runJRAtauworkflow.extend([ - fileNames_and_options_showHistos[sampleName][algorithm][refVariable]['outputFileName'], - fileNames_and_options_showGraphs[sampleName][algorithm][refVariable]['outputFileName'] - ]) + for refVariable in [ "RefPt", "JetEta" ]: + outputFileNames_runJRAtauworkflow.extend([ + fileNames_and_options_showGraphs[algorithm][refVariable]['outputFileName'] + ]) + + makeFile.write("all: %s\n" % make_MakeFile_vstring(outputFileNames_runJRAtauworkflow)) makeFile.write("\techo 'Finished running JRAtau Workflow.'\n") makeFile.write("\n") -for sampleName in samplesToAnalyze: - if len(fileNames_hadd[sampleName]['inputFileNames']) > 0: - makeFile.write("%s: %s\n" % - (fileNames_hadd[sampleName]['outputFileName'], - make_MakeFile_vstring(fileNames_hadd[sampleName]['inputFileNames']))) - makeFile.write("\t%s %s &> %s\n" % - (executable_shell, - fileNames_hadd[sampleName]['shellFileName'], - fileNames_hadd[sampleName]['logFileName'])) -makeFile.write("\n") -for sampleName in samplesToAnalyze: + +# add the ntuple +if len(fileNames_hadd['inputFileNames']) > 0: makeFile.write("%s: %s\n" % - (fileNames_and_options_jrAnalyzer[sampleName]['outputFileName'], - make_MakeFile_vstring(fileNames_and_options_jrAnalyzer[sampleName]['inputFileNames']))) + (fileNames_hadd['outputFileName'], + make_MakeFile_vstring(fileNames_hadd['inputFileNames']))) makeFile.write("\t%s %s &> %s\n" % - (executable_jrAnalyzer, - fileNames_and_options_jrAnalyzer[sampleName]['commandLine'], - fileNames_and_options_jrAnalyzer[sampleName]['logFileName'])) + (executable_shell, + fileNames_hadd['shellFileName'], + fileNames_hadd['logFileName'])) makeFile.write("\n") -for sampleName in samplesToAnalyze: - makeFile.write("%s: %s\n" % - (fileNames_and_options_fitResponse_uncalibrated[sampleName]['outputFileName'], - make_MakeFile_vstring(fileNames_and_options_fitResponse_uncalibrated[sampleName]['inputFileNames']))) - makeFile.write("\t%s %s &> %s\n" % - (executable_fitResponse, - fileNames_and_options_fitResponse_uncalibrated[sampleName]['commandLine'], - fileNames_and_options_fitResponse_uncalibrated[sampleName]['logFileName'])) + +# response +makeFile.write("%s: %s\n" % + (fileNames_and_options_jrAnalyzer['outputFileName'], + make_MakeFile_vstring(fileNames_and_options_jrAnalyzer['inputFileNames']))) +makeFile.write("\t%s %s &> %s\n" % + (executable_jrAnalyzer, + fileNames_and_options_jrAnalyzer['commandLine'], + fileNames_and_options_jrAnalyzer['logFileName'])) makeFile.write("\n") -for sampleName in samplesToAnalyze: - makeFile.write("%s: %s\n" % - (fileNames_and_options_fitL3param[sampleName]['outputFileName'], - make_MakeFile_vstring(fileNames_and_options_fitL3param[sampleName]['inputFileNames']))) - makeFile.write("\t%s %s &> %s\n" % - (executable_fitL3param, - fileNames_and_options_fitL3param[sampleName]['commandLine'], - fileNames_and_options_fitL3param[sampleName]['logFileName'])) + +# fit uncalibrated +for histName in histNames: + makeFile.write("%s: %s\n" % + (fileNames_and_options_fitResponse_uncalibrated[histName]['outputFileName'], + make_MakeFile_vstring(fileNames_and_options_fitResponse_uncalibrated[histName]['inputFileNames']))) + makeFile.write("\t%s %s &> %s\n" % + (executable_fitResponse, + fileNames_and_options_fitResponse_uncalibrated[histName]['commandLine'], + fileNames_and_options_fitResponse_uncalibrated[histName]['logFileName'])) + makeFile.write("\n") + +# add the fit uncalibrated hist +if len(fitHists_hadd_uncalibrated['inputFileNames']) > 0: makeFile.write("%s: %s\n" % - (fileNames_and_options_fitL2param[sampleName]['outputFileName'], - make_MakeFile_vstring(fileNames_and_options_fitL2param[sampleName]['inputFileNames']))) + (fitHists_hadd_uncalibrated['outputFileName'], + make_MakeFile_vstring(fitHists_hadd_uncalibrated['inputFileNames']))) makeFile.write("\t%s %s &> %s\n" % - (executable_fitL2param, - fileNames_and_options_fitL2param[sampleName]['commandLine'], - fileNames_and_options_fitL2param[sampleName]['logFileName'])) + (executable_shell, + fitHists_hadd_uncalibrated['shellFileName'], + fitHists_hadd_uncalibrated['logFileName'])) makeFile.write("\n") -for sampleName in samplesToAnalyze: - makeFile.write("%s: %s\n" % - (fileNames_and_options_applyL2L3param[sampleName]['outputFileName'], - make_MakeFile_vstring(fileNames_and_options_applyL2L3param[sampleName]['inputFileNames']))) - makeFile.write("\t%s %s &> %s\n" % - (executable_applyL2L3param, - fileNames_and_options_applyL2L3param[sampleName]['commandLine'], - fileNames_and_options_applyL2L3param[sampleName]['logFileName'])) + +# determine l2l3 +makeFile.write("%s: %s\n" % + (fileNames_and_options_fitL2param['outputFileName'], + make_MakeFile_vstring(fileNames_and_options_fitL2param['inputFileNames']))) +makeFile.write("\t%s %s &> %s\n" % + (executable_fitL2param, + fileNames_and_options_fitL2param['commandLine'], + fileNames_and_options_fitL2param['logFileName'])) +makeFile.write("\n") + +# apply l2l3 +makeFile.write("%s: %s\n" % + (fileNames_and_options_applyL2L3param['outputFileName'], + make_MakeFile_vstring(fileNames_and_options_applyL2L3param['inputFileNames']))) +makeFile.write("\t%s %s &> %s\n" % + (executable_applyL2L3param, + fileNames_and_options_applyL2L3param['commandLine'], + fileNames_and_options_applyL2L3param['logFileName'])) +makeFile.write("\n") + +# calibrated response +makeFile.write("%s: %s\n" % + (fileNames_and_options_jrAnalyzer_calibrated['outputFileName'], + make_MakeFile_vstring(fileNames_and_options_jrAnalyzer_calibrated['inputFileNames']))) +makeFile.write("\t%s %s &> %s\n" % + (executable_jrAnalyzer, + fileNames_and_options_jrAnalyzer_calibrated['commandLine'], + fileNames_and_options_jrAnalyzer_calibrated['logFileName'])) makeFile.write("\n") -for sampleName in samplesToAnalyze: + +# fit calibrated response +for histName in histNames: + makeFile.write("%s: %s\n" % + (fileNames_and_options_fitResponse_calibrated[histName]['outputFileName'], + make_MakeFile_vstring(fileNames_and_options_fitResponse_calibrated[histName]['inputFileNames']))) + makeFile.write("\t%s %s &> %s\n" % + (executable_fitResponse, + fileNames_and_options_fitResponse_calibrated[histName]['commandLine'], + fileNames_and_options_fitResponse_calibrated[histName]['logFileName'])) + makeFile.write("\n") + +# add the calibrated fit hist +if len(fitHists_hadd_calibrated['inputFileNames']) > 0: makeFile.write("%s: %s\n" % - (fileNames_and_options_jrAnalyzer_calibrated[sampleName]['outputFileName'], - make_MakeFile_vstring(fileNames_and_options_jrAnalyzer_calibrated[sampleName]['inputFileNames']))) + (fitHists_hadd_calibrated['outputFileName'], + make_MakeFile_vstring(fitHists_hadd_calibrated['inputFileNames']))) makeFile.write("\t%s %s &> %s\n" % - (executable_jrAnalyzer, - fileNames_and_options_jrAnalyzer_calibrated[sampleName]['commandLine'], - fileNames_and_options_jrAnalyzer_calibrated[sampleName]['logFileName'])) + (executable_shell, + fitHists_hadd_calibrated['shellFileName'], + fitHists_hadd_calibrated['logFileName'])) makeFile.write("\n") -for sampleName in samplesToAnalyze: + +#add uncalibrated and calibrated fit +if len(fitHists_hadd_all['inputFileNames']) > 0: makeFile.write("%s: %s\n" % - (fileNames_and_options_fitResponse_calibrated[sampleName]['outputFileName'], - make_MakeFile_vstring(fileNames_and_options_fitResponse_calibrated[sampleName]['inputFileNames']))) + (fitHists_hadd_all['outputFileName'], + make_MakeFile_vstring(fitHists_hadd_all['inputFileNames']))) makeFile.write("\t%s %s &> %s\n" % - (executable_fitResponse, - fileNames_and_options_fitResponse_calibrated[sampleName]['commandLine'], - fileNames_and_options_fitResponse_calibrated[sampleName]['logFileName'])) + (executable_shell, + fitHists_hadd_all['shellFileName'], + fitHists_hadd_all['logFileName'])) +makeFile.write("\n") + +# resolution +makeFile.write("%s: %s\n" % + (fileNames_and_options_fitResolution['outputFileName'], + #make_MakeFile_vstring(fitHists_hadd_calibrated['outputFileName']))) + fitHists_hadd_all['outputFileName'])) +makeFile.write("\t%s %s &> %s\n" % + (executable_fitResolution, + fileNames_and_options_fitResolution['commandLine'], + fileNames_and_options_fitResolution['logFileName'])) +makeFile.write("\n") + +# histo graph +for algorithm in algorithms: makeFile.write("%s: %s\n" % - (fileNames_and_options_fitResolution[sampleName]['outputFileName'], - make_MakeFile_vstring(fileNames_and_options_fitResolution[sampleName]['inputFileNames']))) + (fileNames_and_options_showHistos[algorithm]['outputFileName'], + fitHists_hadd_all['outputFileName'])) makeFile.write("\t%s %s &> %s\n" % - (executable_fitResolution, - fileNames_and_options_fitResolution[sampleName]['commandLine'], - fileNames_and_options_fitResolution[sampleName]['logFileName'])) -makeFile.write("\n") -for sampleName in samplesToAnalyze: - for algorithm in algorithms: - for refVariable in [ "RefPt", "JetEta" ]: - makeFile.write("%s: %s\n" % - (fileNames_and_options_showHistos[sampleName][algorithm][refVariable]['outputFileName'], - make_MakeFile_vstring(fileNames_and_options_showHistos[sampleName][algorithm][refVariable]['inputFileNames']))) - makeFile.write("\t%s %s &> %s\n" % - (executable_showHistos, - fileNames_and_options_showHistos[sampleName][algorithm][refVariable]['commandLine'], - fileNames_and_options_showHistos[sampleName][algorithm][refVariable]['logFileName'])) - makeFile.write("%s: %s\n" % - (fileNames_and_options_showGraphs[sampleName][algorithm][refVariable]['outputFileName'], - make_MakeFile_vstring(fileNames_and_options_showGraphs[sampleName][algorithm][refVariable]['inputFileNames']))) - makeFile.write("\t%s %s &> %s\n" % - (executable_showGraphs, - fileNames_and_options_showGraphs[sampleName][algorithm][refVariable]['commandLine'], - fileNames_and_options_showGraphs[sampleName][algorithm][refVariable]['logFileName'])) -makeFile.write("\n") + (executable_showHistos, + fileNames_and_options_showHistos[algorithm]['commandLine'], + fileNames_and_options_showHistos[algorithm]['logFileName'])) + for refVariable in [ "RefPt", "JetEta" ]: + makeFile.write("%s: %s\n" % + (fileNames_and_options_showGraphs[algorithm][refVariable]['outputFileName'], + fileNames_and_options_fitResolution['outputFileName'])) + makeFile.write("\t%s %s &> %s\n" % + (executable_showGraphs, + fileNames_and_options_showGraphs[algorithm][refVariable]['commandLine'], + fileNames_and_options_showGraphs[algorithm][refVariable]['logFileName'])) +makeFile.write("\n") + makeFile.write(".PHONY: clean\n") makeFile.write("clean:\n") makeFile.write("\trm -f %s\n" % make_MakeFile_vstring(outputFileNames_runJRAtauworkflow)) makeFile.write("\techo 'Finished deleting old files.'\n") makeFile.write("\n") makeFile.close() +#-------------------------------------------------------------------------------- + +#print("Finished building Makefile. Now execute 'make -j 16 -f %s'." % makeFileName) +print("Finished building Makefile. Now execute 'make -j -f %s'." % makeFileName) -print("Finished building Makefile. Now execute 'make -j 8 -f %s'." % makeFileName) diff --git a/JetAnalyzers/test/runSVfitMEMWorkflow.py b/JetAnalyzers/test/runSVfitMEMWorkflow.py new file mode 100755 index 00000000..4fdd1dbe --- /dev/null +++ b/JetAnalyzers/test/runSVfitMEMWorkflow.py @@ -0,0 +1,247 @@ +#!/usr/bin/env python + +from variable import * +from function import * + +import os +import re +import subprocess + +inputFilePath = '%s/src/JetMETAnalysis/JetAnalyzers/test/JRAtau' % (os.environ['CMSSW_BASE']) +outputParPath = '/home/calpas/svfitMEM/CMSSW_8_0_0/src/TauAnalysis/SVfitTF/data/L2L3Corr' +print "plots path: %s" %(inputFilePath) +print "parameters path: %s " %(outputParPath) + +algorithms = [ + 'ak5tauHPSlooseCombDBcorrAll', + 'ak5tauHPSlooseCombDBcorrOneProng0Pi0', + #'ak5tauHPSlooseCombDBcorrOneProng1Pi0', + #'ak5tauHPSlooseCombDBcorrTwoProng0Pi0', + #'ak5tauHPSlooseCombDBcorrTwoProng1Pi0', + #'ak5tauHPSlooseCombDBcorrThreeProng0Pi0', + #'ak5tauHPSlooseCombDBcorrThreeProng1Pi0', + ## + #'ak5tauHPSlooseCombDBcorrAlll2', + #'ak5tauHPSlooseCombDBcorrOneProng0Pi0l2', + #'ak5tauHPSlooseCombDBcorrOneProng1Pi0l2', + #'ak5tauHPSlooseCombDBcorrTwoProng0Pi0l2', + #'ak5tauHPSlooseCombDBcorrTwoProng1Pi0l2', + #'ak5tauHPSlooseCombDBcorrThreeProng0Pi0l2', + #'ak5tauHPSlooseCombDBcorrThreeProng1Pi0l2', +] + +#-------------------------------------------------------------------------------- +# +# create cristalBall fit directories +# +outputFilePath = '%s/src/JetMETAnalysis/JetAnalyzers/test/JRAtau/%s' % (os.environ['CMSSW_BASE'], version) +fitDir = outputFilePath+"/plots" +cuts = ["pass", "failed"] +for cut in cuts: + for alg in algorithms: + outputPlotPath = fitDir+'/%s/crystalBallFit/%s/' % (alg, cut) + if not os.path.exists(outputPlotPath): + os.makedirs(outputPlotPath) +#-------------------------------------------------------------------------------- + +#-------------------------------------------------------------------------------- +# +# build Makefile... +# +makeFileName = "Makefile_runSVfitMEMworkflow_l2_%s" % version +makeFile = open(makeFileName, "w") +makeFile.write("\n") +outputFileNames_runJRAtauworkflow = [] +#-------------------------------------------------------------------------------- + +#-------------------------------------------------------------------------------- +# +# add all samples files +# +haddInputFileNames = [] +#for sampleToAnalyze in samplesToAnalyze: +# nbOfSample = [] +# for ntupleFileName in ntupleFileNames: +# if ntupleFile_matcher.match(ntupleFileName) and \ +# ntupleFile_matcher.match(ntupleFileName).group('sample') == sampleToAnalyze: +# haddInputFileNames.append(ntupleFileName) +# nbOfSample.append(ntupleFileName) +# print "sample = %s: found %i input files." % (sampleToAnalyze, len(nbOfSample)) +#print "found %i input files." % (len(haddInputFileNames)) + +haddShellFileName = os.path.join(outputFilePath, 'harvestJRAtauNtuples.csh') +haddOutputFileName = os.path.join(outputFilePath, 'ntupleJRAtau_all.root') # copy where you can +retVal_hadd = buildConfigFile_hadd(executable_hadd, haddShellFileName, haddInputFileNames, haddOutputFileName) + +fileNames_hadd = {} +fileNames_hadd['shellFileName'] = haddShellFileName +fileNames_hadd['inputFileNames'] = haddInputFileNames +fileNames_hadd['outputFileName'] = haddOutputFileName +fileNames_hadd['logFileName'] = retVal_hadd['logFileName'] +#outputFileNames_runJRAtauworkflow.extend([ fileNames_hadd['outputFileName'],]) + + +# makefile +#if len(fileNames_hadd['inputFileNames']) > 0: +# makeFile.write("%s: %s\n" % +# (fileNames_hadd['outputFileName'], +# make_MakeFile_vstring(fileNames_hadd['inputFileNames']))) +# makeFile.write("\t%s %s &> %s\n" % +# (executable_shell, +# fileNames_hadd['shellFileName'], +# fileNames_hadd['logFileName'])) +#makeFile.write("\n") +#-------------------------------------------------------------------------------- + +#-------------------------------------------------------------------------------- +# +# jet response for uncalibrated tau-jets +# +fileNames_and_options_jrAnalyzer = {} +fileNames_and_options_jrAnalyzer['inputFileNames'] = [fileNames_hadd['outputFileName']] +fileNames_and_options_jrAnalyzer['outputFileName'] = os.path.join(outputFilePath, "jet_response_analyzer_uncalibrated.root") +fileNames_and_options_jrAnalyzer['configFileName'] = os.path.join(outputFilePath, "jet_response_analyzer_uncalibrated.cfg") +make_jrAnalyzer_config(fileNames_and_options_jrAnalyzer['configFileName']) +fileNames_and_options_jrAnalyzer['logFileName'] = os.path.join(outputFilePath, "jet_response_analyzer_uncalibrated.log") +fileNames_and_options_jrAnalyzer['commandLine'] = '%s -input %s -output %s -algs %s -nbinsrelrsp 250' % \ + (fileNames_and_options_jrAnalyzer['configFileName'], + make_MakeFile_vstring(fileNames_and_options_jrAnalyzer['inputFileNames']), + fileNames_and_options_jrAnalyzer['outputFileName'], + "".join([ "%s:0.3 " % algorithm for algorithm in algorithms ]) + ) +#outputFileNames_runJRAtauworkflow.extend([ fileNames_and_options_jrAnalyzer['outputFileName'],]) + +# makefile +#makeFile.write("%s: %s\n" % +# (fileNames_and_options_jrAnalyzer['outputFileName'], +# make_MakeFile_vstring(fileNames_and_options_jrAnalyzer['inputFileNames']))) +#makeFile.write("\t%s %s &> %s\n" % +# (executable_jrAnalyzer, +# fileNames_and_options_jrAnalyzer['commandLine'], +# fileNames_and_options_jrAnalyzer['logFileName'])) +#makeFile.write("\n") +#-------------------------------------------------------------------------------- + +#-------------------------------------------------------------------------------- +# +# round 1.1: fit uncalibrated jet response +# +fileNames_and_options_fitResponse_uncalibrated = {} +for alg in algorithms: + haddFitHists = [] + for histName in histNames: + fileNames_and_options_fitResponse_uncalibrated[histName] = {} + fileNames_and_options_fitResponse_uncalibrated[histName]['inputFileNames'] = [fileNames_and_options_jrAnalyzer['outputFileName']] + fileNames_and_options_fitResponse_uncalibrated[histName]['outputFileName'] = os.path.join(fitDir, "%s/responseJRAtau_%s.root" % (alg, histName)) + haddFitHists.append(fileNames_and_options_fitResponse_uncalibrated[histName]['outputFileName']) + fileNames_and_options_fitResponse_uncalibrated[histName]['logFileName'] = os.path.join(fitDir, "%s/responseJRAtau_%s.log" % (alg, histName)) + fileNames_and_options_fitResponse_uncalibrated[histName]['commandLine'] = '-input %s -output %s -algs %s -fittype %i -histName %s -polDeg 1 -normalized true -fitDir %s' % \ + (make_MakeFile_vstring(fileNames_and_options_fitResponse_uncalibrated[histName]['inputFileNames']), + fileNames_and_options_fitResponse_uncalibrated[histName]['outputFileName'], + alg, + fitOption, + histName, + fitDir + ) + outputFileNames_runJRAtauworkflow.extend([ fileNames_and_options_fitResponse_uncalibrated[histName]['outputFileName'],]) + + # makefile + makeFile.write("%s: %s\n" % + (fileNames_and_options_fitResponse_uncalibrated[histName]['outputFileName'], + make_MakeFile_vstring(fileNames_and_options_fitResponse_uncalibrated[histName]['inputFileNames']))) + makeFile.write("\t%s %s &> %s\n" % + (executable_fitResponse, + fileNames_and_options_fitResponse_uncalibrated[histName]['commandLine'], + fileNames_and_options_fitResponse_uncalibrated[histName]['logFileName'])) + makeFile.write("\n") + #-------------------------------------------------------------------------------- + + #-------------------------------------------------------------------------------- + # + # round 1.2: add the uncalibrated fit histograms + # + haddShellFitHistName = os.path.join(fitDir, '%s/harvestFitHist.csh' %(alg)) + haddOutputFitHistName = os.path.join(fitDir, '%s/fitHists.root' %(alg)) + + retFitHist_hadd = buildConfigFile_hadd(executable_haddFitHists, haddShellFitHistName, haddFitHists, haddOutputFitHistName) + + fitHists_hadd_uncalibrated = {} + fitHists_hadd_uncalibrated['shellFileName'] = haddShellFitHistName + fitHists_hadd_uncalibrated['inputFileNames'] = haddFitHists + fitHists_hadd_uncalibrated['outputFileName'] = haddOutputFitHistName + fitHists_hadd_uncalibrated['logFileName'] = retFitHist_hadd['logFileName'] + outputFileNames_runJRAtauworkflow.extend([ fitHists_hadd_uncalibrated['outputFileName'],]) + + # makefile + if len(fitHists_hadd_uncalibrated['inputFileNames']) > 0: + makeFile.write("%s: %s\n" % + (fitHists_hadd_uncalibrated['outputFileName'], + make_MakeFile_vstring(fitHists_hadd_uncalibrated['inputFileNames']))) + makeFile.write("\t%s %s &> %s\n" % + (executable_shell, + fitHists_hadd_uncalibrated['shellFileName'], + fitHists_hadd_uncalibrated['logFileName'])) + makeFile.write("\n") + #-------------------------------------------------------------------------------- + + + #-------------------------------------------------------------------------------- + # + # determine L2L3 corrections and fit them + # + #pars=[6, 17] + #for par in pars: + for par in range(0, 19): # from 0 to 18!! + fileNames_and_options_fitL2param = {} + fileNames_and_options_fitL2param['inputFileNames'] = [fitHists_hadd_uncalibrated['outputFileName']] + + outdir = outputParPath+"/"+alg+"/"+str(par) + subprocess.call(['rm', '-r', outdir]) + subprocess.check_output(['mkdir', '-p', outdir]) + + fileNames_and_options_fitL2param['outputFileName'] = os.path.join(outdir, "fitL2param.root") + outputFileNames_runJRAtauworkflow.extend([ fileNames_and_options_fitL2param['outputFileName'],]) + fileNames_and_options_fitL2param['logFileName'] = os.path.join(outdir, "fitL2param.log") + fileNames_and_options_fitL2param['commandLine'] = \ + '-input %s -output %s -era %s -algs %s -formats png -batch true -l2l3 true -cbPar %i -outputDir %s' % \ + (make_MakeFile_vstring(fileNames_and_options_fitL2param['inputFileNames']), + fileNames_and_options_fitL2param['outputFileName'], + era, + alg, + par, + outdir) + + makeFile.write("%s: %s\n" % + (fileNames_and_options_fitL2param['outputFileName'], + #fileNames_and_options_fitResponse_uncalibrated_pol1['outputFileName'], + fitHists_hadd_uncalibrated['outputFileName'])) + makeFile.write("\t%s %s &> %s\n" % + (executable_fitL2param, + fileNames_and_options_fitL2param['commandLine'], + fileNames_and_options_fitL2param['logFileName'])) + makeFile.write("\n") + +#makeFile.write(".PHONY: clean\n") +#makeFile.write("clean:\n") +#makeFile.write("\trm -f %s\n" % make_MakeFile_vstring(outputFileNames_runJRAtauworkflow)) +#makeFile.write("\techo 'Finished deleting old files.'\n") +#makeFile.write("\n") + +makeFile.close() + + +#-------------------------------------------------------------------------------- +# +# write the all target on top of the make file +# +append_copy = open(makeFileName, "r") +original_text = append_copy.read() +append_copy.close() +append_copy = open(makeFileName, "w") +append_copy.write("all: %s\n" % make_MakeFile_vstring(outputFileNames_runJRAtauworkflow)) +append_copy.write(original_text) +append_copy.close() + +#print("Finished building Makefile. Now execute 'make -j 16 -f %s'." % makeFileName) +print("Finished building Makefile. Now execute 'make -j 23 -f %s'." % makeFileName) + diff --git a/JetAnalyzers/test/variable.py b/JetAnalyzers/test/variable.py new file mode 100644 index 00000000..77843e34 --- /dev/null +++ b/JetAnalyzers/test/variable.py @@ -0,0 +1,296 @@ + +import os +import re +import subprocess + + +#-------------------------------------------------------------------------------- +# +# executable +# +execDir = "%s/bin/%s/" % (os.environ['CMSSW_BASE'], os.environ['SCRAM_ARCH']) +executable_jrAnalyzer = execDir + 'jet_response_analyzer_x' +executable_fitResponse = execDir + 'jet_response_fitter_x' +executable_fitL2param = execDir + 'jet_l2_correction_x' # using option here to process L3 correction in one go!! +executable_applyL2L3param = execDir + 'jet_apply_jec_x' +executable_fitResolution = execDir + 'jet_response_and_resolution_x' +executable_showGraphs = execDir + 'jet_inspect_graphs_x' +executable_showHistos = execDir + 'jet_inspect_histos_x' +#executable_showProfiles = execDir + 'jet_inspect_profiles_x' +executable_hadd = 'hadd -f -k' # -f(target != source), -k(skip corrupt files) +executable_haddFitHists = 'hadd -f' +executable_shell = '/bin/csh' +executable_python = 'python' +#-------------------------------------------------------------------------------- + +toCal = [ + 'uncalibrated', + #'calibrated', +] + + +#-------------------------------------------------------------------------------- +# correction file directory +# +jecTextFilePath = os.getcwd() +#-------------------------------------------------------------------------------- + + +#-------------------------------------------------------------------------------- +# +# define function used for fitting tau-jet response: +# 0 = Gaussian +# 1 = Crystal-Ball function +# 2 = Calpas Veelken CrystalBall +# +fitOption = 2 +#-------------------------------------------------------------------------------- + + + +ptBinning = [ 15., 20., 22.5, 25., 27.5, 30., 35., 40., + 45., 50., 60., 80., 120., 200., 500., 3000. ] #15 + +etaBinning = [ -2.3, -2.1, -1.9, -1.7, -1.5, -1.3, + -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, + +0.1, +0.3, +0.5, +0.7, +0.9, +1.1, + +1.3, +1.5, +1.7, +1.9, +2.1, +2.3 ] #23 + + +version = 'v1_2enRecoveryCBa' + +era = 'TauJec11V1' + + +inputFilePaths = [ + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/GluGluHToTauTauM125/', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/VBFHToTauTauM125', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/DYJetsToLLM50', + # + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/ZprimeToTauTauM500', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/ZprimeToTauTauM1000', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/ZprimeToTauTauM1500', + #'/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/ZprimeToTauTauM2000', # all jobs failed!! + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/ZprimeToTauTauM2500', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/ZprimeToTauTauM3500', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/ZprimeToTauTauM4000', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/ZprimeToTauTauM4500', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/ZprimeToTauTauM5000', + # + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM1000', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM1200', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM1400', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM1600', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM1800', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM2000', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM2200', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM2400', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM2600', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM2800', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM3200', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM3400', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM3600', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM3800', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM4000', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM4200', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM4400', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM4600', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM4800', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM5000', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM5200', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM5400', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM5600', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/WprimeToTauNuM5800', + # + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM80', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM100', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM110', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM120', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM130', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM140', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM160', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM180', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM200', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM250', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM300', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM350', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM400', + #'/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM450', #crab: WARNING: No sites are hosting any part of data for block: + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM500', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM800', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM900', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM1000', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM1200', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM1500', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM1600', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM1800', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM2000', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM2300', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM2600', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM2900', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToHToTauTauM3200', + # + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM80', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM90', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM120', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM130', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM140', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM180', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM200', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM250', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM300', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM350', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM400', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM450', + #'/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM500', # job failed!! + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM700', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM800', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM900', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM1000', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM1200', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM1400', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM1500', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM1800', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM2300', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM2600', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM2900', + '/hdfs/cms/store/user/calpas/JRAtauNtuple/v4/SUSYGluGluToBBHToTauTauM3200', +] + + +#-------------------------------------------------------------------------------- +# +# build shell script for running 'hadd' in order to "harvest" histograms +# +ntupleFileNames = [] +for inputFilePath in inputFilePaths: + ntupleFileNames.extend([ os.path.join(inputFilePath, file) for file in os.listdir(inputFilePath) ]) + +ntupleFile_regex = r"[a-zA-Z0-9_/:.]*ntupleJRAtau_(?P[a-zA-Z0-9]*)_(?P\d*)_(?P\d*)_(?P[a-zA-Z0-9]*).root" +ntupleFile_matcher = re.compile(ntupleFile_regex) +#-------------------------------------------------------------------------------- + + +samplesToAnalyze = [ + 'VBFHToTauTauM125', + 'GluGluHToTauTauM125', + 'DYJetsToLLM50', + 'ZprimeToTauTauM500', + 'ZprimeToTauTauM1000', + 'ZprimeToTauTauM1500', + #'ZprimeToTauTauM2000', + 'ZprimeToTauTauM2500', + 'ZprimeToTauTauM3500', + 'ZprimeToTauTauM4000', + 'ZprimeToTauTauM4500', + 'ZprimeToTauTauM5000', + 'WprimeToTauNuM1000', + 'WprimeToTauNuM1200', + 'WprimeToTauNuM1400', + 'WprimeToTauNuM1600', + 'WprimeToTauNuM1800', + 'WprimeToTauNuM2000', + 'WprimeToTauNuM2200', + 'WprimeToTauNuM2400', + 'WprimeToTauNuM2600', + 'WprimeToTauNuM2800', + 'WprimeToTauNuM3200', + 'WprimeToTauNuM3400', + 'WprimeToTauNuM3600', + 'WprimeToTauNuM3800', + 'WprimeToTauNuM4000', + 'WprimeToTauNuM4200', + 'WprimeToTauNuM4400', + 'WprimeToTauNuM4600', + 'WprimeToTauNuM4800', + 'WprimeToTauNuM5000', + 'WprimeToTauNuM5200', + 'WprimeToTauNuM5400', + 'WprimeToTauNuM5600', + 'WprimeToTauNuM5800', + # + 'SUSYGluGluToHToTauTauM80', + 'SUSYGluGluToHToTauTauM100', + 'SUSYGluGluToHToTauTauM110', + 'SUSYGluGluToHToTauTauM120', + 'SUSYGluGluToHToTauTauM130', + 'SUSYGluGluToHToTauTauM140', + 'SUSYGluGluToHToTauTauM160', + 'SUSYGluGluToHToTauTauM180', + 'SUSYGluGluToHToTauTauM200', + 'SUSYGluGluToHToTauTauM250', + 'SUSYGluGluToHToTauTauM300', + 'SUSYGluGluToHToTauTauM350', + 'SUSYGluGluToHToTauTauM400', + #'SUSYGluGluToHToTauTauM450', + 'SUSYGluGluToHToTauTauM500', + 'SUSYGluGluToHToTauTauM800', + 'SUSYGluGluToHToTauTauM900', + 'SUSYGluGluToHToTauTauM1000', + 'SUSYGluGluToHToTauTauM1200', + 'SUSYGluGluToHToTauTauM1500', + 'SUSYGluGluToHToTauTauM1600', + 'SUSYGluGluToHToTauTauM1800', + 'SUSYGluGluToHToTauTauM2000', + 'SUSYGluGluToHToTauTauM2300', + 'SUSYGluGluToHToTauTauM2600', + 'SUSYGluGluToHToTauTauM2900', + 'SUSYGluGluToHToTauTauM3200', + # + 'SUSYGluGluToBBHToTauTauM80', + 'SUSYGluGluToBBHToTauTauM90', + 'SUSYGluGluToBBHToTauTauM120', + 'SUSYGluGluToBBHToTauTauM130', + 'SUSYGluGluToBBHToTauTauM140', + 'SUSYGluGluToBBHToTauTauM180', + 'SUSYGluGluToBBHToTauTauM200', + 'SUSYGluGluToBBHToTauTauM250', + 'SUSYGluGluToBBHToTauTauM300', + 'SUSYGluGluToBBHToTauTauM350', + 'SUSYGluGluToBBHToTauTauM400', + 'SUSYGluGluToBBHToTauTauM450', + #'SUSYGluGluToBBHToTauTauM500', + 'SUSYGluGluToBBHToTauTauM700', + 'SUSYGluGluToBBHToTauTauM800', + 'SUSYGluGluToBBHToTauTauM900', + 'SUSYGluGluToBBHToTauTauM1000', + 'SUSYGluGluToBBHToTauTauM1200', + 'SUSYGluGluToBBHToTauTauM1400', + 'SUSYGluGluToBBHToTauTauM1500', + 'SUSYGluGluToBBHToTauTauM1800', + 'SUSYGluGluToBBHToTauTauM2300', + 'SUSYGluGluToBBHToTauTauM2600', + 'SUSYGluGluToBBHToTauTauM2900', + 'SUSYGluGluToBBHToTauTauM3200', +] + + + +# for time reason rn parallel jobs for different eta for the crystalball fit +histNames = [ + "RelRsp_JetEta-2.3to-2.1_RefPt", + "RelRsp_JetEta-2.1to-1.9_RefPt", + "RelRsp_JetEta-1.9to-1.7_RefPt", + "RelRsp_JetEta-1.7to-1.5_RefPt", + "RelRsp_JetEta-1.5to-1.3_RefPt", + "RelRsp_JetEta-1.3to-1.1_RefPt", + "RelRsp_JetEta-1.1to-0.9_RefPt", + "RelRsp_JetEta-0.9to-0.7_RefPt", + "RelRsp_JetEta-0.7to-0.5_RefPt", + "RelRsp_JetEta-0.5to-0.3_RefPt", + "RelRsp_JetEta-0.3to-0.1_RefPt", + "RelRsp_JetEta-0.1to0.1_RefPt", + "RelRsp_JetEta0.1to0.3_RefPt", + "RelRsp_JetEta0.3to0.5_RefPt", + "RelRsp_JetEta0.5to0.7_RefPt", + "RelRsp_JetEta0.7to0.9_RefPt", + "RelRsp_JetEta0.9to1.1_RefPt", + "RelRsp_JetEta1.1to1.3_RefPt", + "RelRsp_JetEta1.3to1.5_RefPt", + "RelRsp_JetEta1.5to1.7_RefPt", + "RelRsp_JetEta1.7to1.9_RefPt", + "RelRsp_JetEta1.9to2.1_RefPt", + "RelRsp_JetEta2.1to2.3_RefPt", + ] + + diff --git a/JetUtilities/interface/L2Creator.hh b/JetUtilities/interface/L2Creator.hh index a77573dc..61638e06 100644 --- a/JetUtilities/interface/L2Creator.hh +++ b/JetUtilities/interface/L2Creator.hh @@ -14,6 +14,7 @@ #include "JetMETAnalysis/JetUtilities/interface/ObjectLoader.h" #include "JetMETAnalysis/JetUtilities/interface/RootStyle.h" #include "JetMETAnalysis/JetUtilities/interface/JetInfo.hh" +#include "JetMETAnalysis/JetUtilities/interface/crystalBallRes.h" #include "TApplication.h" #include "TFile.h" @@ -30,6 +31,7 @@ #include "TLatex.h" #include "TError.h" #include "TSpectrum.h" +#include "TMath.h" #include #include @@ -38,6 +40,8 @@ #include #include #include +#include +#include #include using namespace std; @@ -76,14 +80,14 @@ private: TString output, outputDir, l2calofit, l2pffit; vector formats, algs; bool l2l3, mpv, delphes; - int maxFitIter; + int maxFitIter, cbPar; TFile* ofile; TFile* ifile; TFile* l3file; TDirectoryFile* l3dir; TDirectoryFile* odir; TGraphErrors* gl3rsp; - TF1* fl3rsp; + TF1* fl3rsp; JetInfo* ji; ObjectLoader hl_rsp; ObjectLoader hl_refpt; @@ -93,4 +97,4 @@ private: vector vrelcor_eta; }; -#endif \ No newline at end of file +#endif diff --git a/JetUtilities/interface/crystalBall.h b/JetUtilities/interface/crystalBall.h new file mode 100644 index 00000000..552877f8 --- /dev/null +++ b/JetUtilities/interface/crystalBall.h @@ -0,0 +1,19 @@ +#ifndef CRYSTALBALL_HH +#define CRYSTALBALL_HH + +// sum of 3gaus+pol7 +double sum_gaus_fnc(double *x, double *par); +// sum of 3gaus+pol7 +double sum_gaus_pol_fnc(double *x, double *par); +// exp +double exp_fcn(double x, double norm, double slope); +// gauss +double gauss_fcn(double x, double norm, double mean, double sigma); +// pol7 +double pol_fcn(double x, double p0, double p1=0, double p2=0, double p3=0, double p4=0, double p5=0, double p6=0, double p7=0); +// global fits +double crystalBall(double *xx, double *pp); +double crystalBall_1(double *xx, double *pp); + + +#endif diff --git a/JetUtilities/interface/crystalBallRes.h b/JetUtilities/interface/crystalBallRes.h new file mode 100644 index 00000000..bba927a6 --- /dev/null +++ b/JetUtilities/interface/crystalBallRes.h @@ -0,0 +1,18 @@ +#ifndef CRYSTALBALLRES_HH +#define CRYSTALBALLRES_HH + +#include +#include +#include "TCanvas.h" +#include + +typedef std::vector vdouble; +typedef std::vector vbool; +typedef std::vector vint; + +void crystalBallRes( TH1 *hrsp, bool drawCrystalBall, vdouble &vintegral, vdouble &vnorm, + vdouble &vrange, vdouble &vmean, vdouble &vsigma, + vdouble &vexp, vdouble &vpolLeft, vdouble &vpolRight, + vint &voption, int polDeg=1, TCanvas *canvas=0); + +#endif diff --git a/JetUtilities/src/L2Creator.cc b/JetUtilities/src/L2Creator.cc index 27474a23..050bb434 100644 --- a/JetUtilities/src/L2Creator.cc +++ b/JetUtilities/src/L2Creator.cc @@ -23,6 +23,7 @@ L2Creator::L2Creator() { mpv = false; delphes = false; maxFitIter = 30; + cbPar = -1; } //______________________________________________________________________________ @@ -34,7 +35,7 @@ L2Creator::L2Creator(CommandLine& cl) { era = cl.getValue ("era"); l3input = cl.getValue ("l3input", "l3.root"); output = cl.getValue ("output", "l2.root"); - outputDir = cl.getValue ("outputDir", "./"); + outputDir = cl.getValue ("outputDir", ""); formats = cl.getVector ("formats", ""); algs = cl.getVector ("algs", ""); l2l3 = cl.getValue ("l2l3", true); @@ -43,6 +44,7 @@ L2Creator::L2Creator(CommandLine& cl) { mpv = cl.getValue ("mpv", false); delphes = cl.getValue ("delphes", false); maxFitIter = cl.getValue ("maxFitIter", 30); + cbPar = cl.getValue ("cbPar", -1); if (!cl.partialCheck()) return; cl.print(); @@ -50,32 +52,31 @@ L2Creator::L2Creator(CommandLine& cl) { void L2Creator::resetForNextAlgorithm() { if(ji!=0) - delete ji; + delete ji; if(odir!=0) - delete odir; + delete odir; if(l3dir!=0) - delete l3dir; - if(gl3rsp!=0) + delete l3dir; + if(gl3rsp!=0) delete gl3rsp; - if(fl3rsp!=0) + if(fl3rsp!=0) delete fl3rsp; hl_rsp.reset(); hl_refpt.reset(); - hl_jetpt.reset(); + hl_jetpt.reset(); - vabsrsp_eta.erase(vabsrsp_eta.begin(),vabsrsp_eta.end()); + vabsrsp_eta.erase(vabsrsp_eta.begin(),vabsrsp_eta.end()); vabscor_eta.erase(vabscor_eta.begin(),vabscor_eta.end()); vrelcor_eta.erase(vrelcor_eta.begin(),vrelcor_eta.end()); } //______________________________________________________________________________ void L2Creator::openOutputFile() { - if(!outputDir.EndsWith("/")) outputDir+="/"; - ofile = TFile::Open(outputDir+output,"RECREATE"); + ofile = TFile::Open(output,"RECREATE"); if (!ofile->IsOpen()) { - cout<<"Can't create "<IsOpen()); + cout<<"Can't create "<IsOpen()); } } @@ -83,8 +84,8 @@ void L2Creator::openOutputFile() { void L2Creator::openInputFile() { ifile = TFile::Open(input.c_str(),"READ"); if (!ifile->IsOpen()) { - cout<<"Can't open "<IsOpen()); + cout<<"Can't open "<IsOpen()); } } @@ -121,7 +122,7 @@ void L2Creator::loopOverAlgorithms() { if (strcmp(dirKey->GetClassName(),"TDirectoryFile")!=0) continue; TDirectoryFile* idir = (TDirectoryFile*)dirKey->ReadObj(); string alg(idir->GetName()); if (!contains(algs,alg)) continue; - + algIndex++; cout << alg << " ... " << endl; @@ -146,21 +147,19 @@ void L2Creator::loopOverAlgorithms() { bool gotL3Rsp = getL3Rsp(); if(!gotL3Rsp) continue; }//if (!l2l3) - // // Make output directory in the file l2.root (or other name if reset by user) // odir = (TDirectoryFile*)ofile->mkdir(alg.c_str()); odir->cd(); - + // // Load the input histograms from jra.root or jra_f.root (or other name if reset by user) // hl_rsp.load_objects(idir,"RelRsp:JetEta:RefPt"); hl_refpt.load_objects(idir,"RefPt:JetEta:RefPt"); hl_jetpt.load_objects(idir,"JetPt:JetEta:RefPt"); - - + // // Absolute response/correction as a function of pT for each eta bin // Needed for both L2 only corrections and L2L3 corrections @@ -168,218 +167,316 @@ void L2Creator::loopOverAlgorithms() { loopOverEtaBins(); if(!l2l3) { - // - // Relative (L2) response/correction as a function of pT for each eta bin - // Needed for L2 only corrections, but not the L2L3 corrections - // - doRelCorFits(); - } + // + // Relative (L2) response/correction as a function of pT for each eta bin + // Needed for L2 only corrections, but not the L2L3 corrections + // + doRelCorFits(); + } // // write the L2 correction text file for the current algorithm // writeTextFileForCurrentAlgorithm(); - + cout<getAbbreviation(); - vector indices; - TH1F* hrsp(0); - hl_rsp.begin_loop(); + string alg = ji->getAbbreviation(); + vector indices; + TH1F* hrsp(0); + hl_rsp.begin_loop(); - while ((hrsp=hl_rsp.next_object(indices))) { + while ((hrsp=hl_rsp.next_object(indices))) { - unsigned int ieta=indices[0]; - unsigned int ipt =indices[1]; + unsigned int ieta=indices[0]; + unsigned int ipt =indices[1]; // // create new graphs if a new eta bin comes around // - if (ipt==0) { - vabsrsp_eta.push_back(new TGraphErrors()); - vabscor_eta.push_back(new TGraphErrors()); - stringstream ss; - ss<SetName(("AbsRspVsRefPt_JetEta"+ss.str()).c_str()); - vabscor_eta.back()->SetName(("AbsCorVsJetPt_JetEta"+ss.str()).c_str()); - } + if (ipt==0) { + vabsrsp_eta.push_back(new TGraphErrors()); + vabscor_eta.push_back(new TGraphErrors()); + stringstream ss; + ss<SetName(("AbsRspVsRefPt_JetEta"+ss.str()).c_str()); + vabscor_eta.back()->SetName(("AbsCorVsJetPt_JetEta"+ss.str()).c_str()); + } + // // only add points to the graphs if the current histo is not empty // the current setting might be a little high // - if (hrsp->GetEntries() > 10) {//hrsp->Integral()!=0) { - - TF1* frsp = (TF1*)hrsp->GetListOfFunctions()->Last(); - //std::cout << "hrspName = " << hrsp->GetName() << ": frsp = " << frsp << std::endl; - TH1F* hrefpt = hl_refpt.object(indices); - TH1F* hjetpt = hl_jetpt.object(indices); - - assert(hrefpt->GetEntries()>0&&hjetpt->GetEntries()>0); - - double refpt =hrefpt->GetMean(); - double erefpt =hrefpt->GetMeanError(); - double jetpt =hjetpt->GetMean(); - double ejetpt =hjetpt->GetMeanError(); - - double peak; - double epeak; - if(alg.find("calo")!=string::npos) { - peak = (frsp==0 || !mpv)?hrsp->GetMean():frsp->GetParameter(1); - epeak = (frsp==0 || !mpv)?hrsp->GetMeanError():frsp->GetParError(1); - } - else if(alg.find("pf")!=string::npos) { - peak = (frsp==0 || !mpv)?hrsp->GetMean():frsp->GetParameter(1); - epeak = (frsp==0 || !mpv)?hrsp->GetMeanError():frsp->GetParError(1); - } - else { - peak = (frsp==0 || !mpv)?hrsp->GetMean():frsp->GetParameter(1); - epeak = (frsp==0 || !mpv)?hrsp->GetMeanError():frsp->GetParError(1); - } - - double absrsp = peak; - double eabsrsp = epeak; - double abscor = 0.0; - double eabscor = 0.0; - - if (absrsp > 0) { - abscor =1.0/absrsp; - eabscor = abscor*abscor*epeak; - } - if ((abscor>0) && (absrsp>0) && (eabscor>1e-5) && (eabscor/abscor<0.5) && (eabsrsp>1e-4) && (eabsrsp/absrsp<0.5)) { - int n = vabsrsp_eta.back()->GetN(); - vabsrsp_eta.back()->SetPoint (n,refpt, absrsp); - vabsrsp_eta.back()->SetPointError(n,erefpt,eabsrsp); - vabscor_eta.back()->SetPoint (n,jetpt, abscor); - vabscor_eta.back()->SetPointError(n,ejetpt,eabscor); - } - else cout << "absrsp " << absrsp << " and eabsrsp " << eabsrsp << " and abscor " << abscor << " and eabscor " << eabscor << endl; - } + if (hrsp->GetEntries()<1) continue; + + TF1* frsp {0}; + frsp = (TF1*)hrsp->GetListOfFunctions()->Last(); + //std::cout << "hrspName = " << hrsp->GetName() << ": frsp = " << frsp << std::endl; + TH1F* hrefpt = hl_refpt.object(indices); + TH1F* hjetpt = hl_jetpt.object(indices); + + assert(hrefpt->GetEntries()>0&&hjetpt->GetEntries()>0); + + double refpt {hrefpt->GetMean()}; + double erefpt{hrefpt->GetMeanError()}; + double jetpt {hjetpt->GetMean()}; + double ejetpt{hjetpt->GetMeanError()}; + + double peak{0}, epeak{0}; + + if(alg.find("calo")!=string::npos) { + peak = (frsp==0 || !mpv)?hrsp->GetMean():frsp->GetParameter(1); + epeak = (frsp==0 || !mpv)?hrsp->GetMeanError():frsp->GetParError(1); + } + else if(alg.find("pf")!=string::npos) { + peak = (frsp==0 || !mpv)?hrsp->GetMean():frsp->GetParameter(1); + epeak = (frsp==0 || !mpv)?hrsp->GetMeanError():frsp->GetParError(1); + } + else if(alg.find("tau")!=string::npos) { + + if( frsp!=0 && frsp->Eval(1)!=0 ) { + + // retrieve cb norm integral + vector vintegral(3), vnorm(3), vrange(3), vmean(3), vsigma(4), vexp(2), vpolLeft(8), vpolRight(8); + vector voption(6); + + crystalBallRes(hrsp, false, vintegral, vnorm, vrange, vmean, vsigma, vexp, vpolLeft, vpolRight, voption); // hist fit, drawfit, intg2, g3, g3, canvas(defaut) + + // integral of the 3 gauss + double sumIntG{0}; for( auto &i : vintegral) sumIntG+=i; + + // norm of the 3 gauss + double sumNormG{0}; for( auto &i : vnorm) sumNormG+=i; + + if(cbPar!=-1){ + // peak + peak = frsp->GetParameter(cbPar); + + // norm Gaus + if(cbPar==0||cbPar==10||cbPar==13) { + if(frsp->GetParameter(cbPar)<0.25*sumIntG || frsp->GetParameter(cbPar)<0.01) epeak = max(7., frsp->GetParError(cbPar)); + else epeak = max(0.5*frsp->GetParameter(cbPar), frsp->GetParError(cbPar)); + } + // mean Gaus + else if(cbPar==1||cbPar==11||cbPar==14) { + if (cbPar==1 && vintegral[0]<0.25*sumIntG) epeak = max(0.15, frsp->GetParError(cbPar)); + else if(cbPar==11 && vintegral[1]<0.25*sumIntG) epeak = max(0.15, frsp->GetParError(cbPar)); + else if(cbPar==14 && vintegral[2]<0.25*sumIntG) epeak = max(0.15, frsp->GetParError(cbPar)); + else epeak = max(0.1, frsp->GetParError(cbPar)); + } + // sigma Gaus + else if(cbPar==2||cbPar==9||cbPar==12||cbPar==15) { + if (cbPar==2 && vintegral[0]<0.25*sumIntG) epeak = max(0.15, frsp->GetParError(cbPar)); + else if(cbPar==9 && vintegral[0]<0.25*sumIntG) epeak = max(0.15, frsp->GetParError(cbPar)); + else if(cbPar==12 && vintegral[1]<0.25*sumIntG) epeak = max(0.15, frsp->GetParError(cbPar)); + else if(cbPar==15 && vintegral[2]<0.25*sumIntG) epeak = max(0.15, frsp->GetParError(cbPar)); + else epeak = max(0.1, frsp->GetParError(cbPar)); + } + //slope exp + else if(cbPar==4) { + epeak = max(5., frsp->GetParError(cbPar)); + } + // range + else if( cbPar==7 || cbPar==8 || cbPar==18) { + if (frsp->GetParameter(cbPar)<0) { peak=0; epeak = 1;} + if (frsp->GetParameter(cbPar)>2) { peak=2; epeak = 1;} + else epeak = max(0.1, frsp->GetParError(cbPar)); + } + // pol p1*x + else if(cbPar==6 || cbPar==17){ + if(frsp->GetParameter(cbPar)>20) { peak=0; epeak=max(2., frsp->GetParError(cbPar));} + else if(frsp->GetParameter(cbPar)<-20){ peak=0; epeak=max(2., frsp->GetParError(cbPar));} + else epeak = max(0.01, frsp->GetParError(cbPar)); + } + } + else { + double histPeak{hrsp->GetXaxis()->GetBinCenter(hrsp->GetMaximumBin())}; + + // select the closest mean to 1 only if the corresponding norm and integral is>25% of tot gauss + //double closestToOne{100}; + //for (unsigned i=0; i0.25*sumIntG || vNormG[i]>0.25*sumNormG) && vmean[i]>=0 && vmean[i]<=2) { + // peak=vmean[i]; closestToOne=dist; + // } + //} + //if(closestToOne==100) cout<<"NO GAUSS 1||2||3 mean selected!! using histogram peak..."<25% of tot gauss + double closestToHistPeak{100}; + for (unsigned i=0; i0.25*sumIntG || vnorm[i]>0.25*sumNormG) && vmean[i]>=0 && vmean[i]<=2) { + peak=vmean[i]; closestToHistPeak=dist; + } + } + if(closestToHistPeak==100) cout<<"NO GAUSS 1||2||3 mean selected!! using histogram peak..."<GetMean():frsp->GetParameter(1); + epeak = (frsp==0 || !mpv)?hrsp->GetMeanError():frsp->GetParError(1); + } + + double absrsp = peak; + double eabsrsp = epeak; + double abscor = 0.0; + double eabscor = 0.0; + + if(cbPar!=-1) { cout<<"cbpar!=-1"<GetN(); + vabsrsp_eta.back()->SetPoint (n, refpt, absrsp); + vabsrsp_eta.back()->SetPointError(n, erefpt, eabsrsp); + vabscor_eta.back()->SetPoint (n, jetpt, abscor); + vabscor_eta.back()->SetPointError(n, ejetpt, eabscor); + + cout<<"\njetPt/RefPt: "<GetN()!=0 && (vabscor_eta.back())->GetN()!=0) { - cout << "Doing fits for " << vabscor_eta.back()->GetName() << " ... " << endl; - TGraphErrors* gabsrsp = vabsrsp_eta.back(); - TGraphErrors* gabscor = vabscor_eta.back(); - TF1* fabscor(0); - int npoints = gabscor->GetN(); - double xmin(1.0),xmax(100.0); - if (npoints > 0) { - // - // we don't want to fit for pt less than 5 GeV as even a corrected calo jet of 10 - // will be at least 5 Gev in raw energy. - // - //xmin = gabscor->GetX()[0]; - //xmin = max(gabscor->GetX()[0],3.0); - xmin = max(gabscor->GetX()[0],6.0); - xmax = gabscor->GetX()[gabscor->GetN()-1]; - } - - if (npoints<3 && !delphes) { - gabscor->SetPoint (0, 10.0,1.0); - gabscor->SetPointError(0, 0.0,0.0); - gabscor->SetPoint (1,100.0,1.0); - gabscor->SetPointError(1, 0.0,0.0); - fabscor = new TF1("fit","[0]",10.0,100.0); - fabscor->FixParameter(0,1.0); - } - else if (!l2l3 && npoints > 2 && gabscor->GetN()<10) { - fabscor=new TF1("fit","[0]+[1]*log10(x)+[2]*pow(log10(x),2)",xmin,xmax); - fabscor->SetParameter(0,1.0); - fabscor->SetParameter(1,0.0); - fabscor->SetParameter(2,0.0); - } - else { - if (alg.find("pf")!=string::npos || alg.find("puppi")!=string::npos) { - // - // Delphes - // - if(delphes) { - //fcn = "[0]+[1]*log10(x)+[2]*pow(log10(x),2)+[3]*pow(log10(x),3)+[4]*pow(log10(x),4)+[5]*pow(log10(x),5)+[6]*pow(log10(x),6)+[7]*pow(log10(x),7)+[8]*pow(log10(x),8)+[9]*pow(log10(x),9)"; - TString fcn = "[0]+[1]*log10(x)+[2]*pow(log10(x),2)+([3]/pow(log10(x),3))+([4]/pow(log10(x),4))+([5]/pow(log10(x),5))"; - fabscor=new TF1("fit",fcn.Data(),xmin,xmax); - } - - // - // online (HLT) - // - if(alg.find("HLT")!=string::npos){ - fabscor=new TF1("fit","(x>=[6])*([0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5])))+(x<[6])*[7]",xmin,xmax); - fabscor->FixParameter(6,xmin); - fabscor->FixParameter(7,0.0); - } - // - // offline - // - else { - if (xmin<10) xmin=10; - TString fcn = getOfflinePFFunction(); - fabscor=new TF1("fit",fcn.Data(),xmin,xmax); - setOfflinePFParameters(gabscor, fabscor,xmin); - } - } - else if (alg.find("trk")!=string::npos) { - fabscor=new TF1("fit","[0]+[1]*pow(x/500.0,[2])+[3]/log10(x)+[4]*log10(x)",xmin,xmax); - fabscor->SetParameter(0,1.7); - fabscor->SetParameter(1,0.7); - fabscor->SetParameter(2,3.0); - fabscor->SetParLimits(2,1,10); - fabscor->SetParameter(3,0.0); - fabscor->SetParameter(4,0.0); - } - else if (alg.find("jpt")!=string::npos || alg.find("tau")!=string::npos) { - fabscor=new TF1("fit","[0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5]))",xmin,xmax); - // - // INITIAL PARAMS: These are the fitted parameters that work for JetEta0.261to0.348 - // - fabscor->SetParameter(0,-8.9); - fabscor->SetParameter(1,530); - fabscor->SetParameter(2,16); - fabscor->SetParameter(3,-22); - fabscor->SetParameter(4,0.06); - fabscor->SetParameter(5,-0.28); - } - else if (alg.find("calo")!=string::npos) { - if(l2calofit.EqualTo("standard",TString::kIgnoreCase)) { - if (xmin<6) xmin=6; - fabscor=new TF1("fit","[0]+[1]/(pow(log10(x),[2])+[3])",xmin,xmax); - fabscor->SetParameter(0,1.0); - fabscor->SetParameter(1,5.0); - fabscor->SetParameter(2,3.0); - fabscor->SetParameter(3,3.0); - - fabscor->SetParLimits(3,0,100); - } - else if(l2calofit.EqualTo("DynamicMin",TString::kIgnoreCase)) { - xmin=findPeak(gabscor,0,3,3,false); - TString fcn = "[0] + ([1]/(pow(log10(x),2)+[2])) + ([3]*exp(-([4]*((log10(x)-[5])*(log10(x)-[5]))))) + ([6]*exp(-([7]*((log10(x)-[8])*(log10(x)-[8])))))"; - fabscor=new TF1("fit",fcn.Data(),xmin,xmax); - fabscor->SetParameter(0,-0.0221278); - fabscor->SetParameter(1,119.265); - fabscor->SetParameter(2,100); - fabscor->SetParameter(3,-0.0679365); - fabscor->SetParameter(4,2.82597); - fabscor->SetParameter(5,1.8277); - fabscor->SetParameter(6,-0.0679365); - fabscor->SetParameter(7,3.82597); - fabscor->SetParameter(8,1.8277); - fabscor->SetParLimits(6,-20,10); - fabscor->SetParLimits(7,0,100); - fabscor->SetParLimits(3,-15,15); - fabscor->SetParLimits(4,0,500); - fabscor->SetParLimits(0,-2,25); - fabscor->SetParLimits(1,0,250); - } - } - else { - cout << "WARNING::Cannot determine fit function for " << alg << "." << endl; - } + if (ipt==hl_jetpt.nobjects(1)-1 && (vabsrsp_eta.back())->GetN()!=0 && (vabscor_eta.back())->GetN()!=0) { + cout << "Doing fits for " << vabscor_eta.back()->GetName() << " ... " << endl; + TGraphErrors* gabsrsp = vabsrsp_eta.back(); + TGraphErrors* gabscor = vabscor_eta.back(); + TF1* fabscor(0); + int npoints = gabscor->GetN(); + double xmin(1.0),xmax(100.0); + if (npoints > 0) { + // + // we don't want to fit for pt less than 5 GeV as even a corrected calo jet of 10 + // will be at least 5 Gev in raw energy. + // + //xmin = gabscor->GetX()[0]; + //xmin = max(gabscor->GetX()[0],3.0); + xmin = max(gabscor->GetX()[0],6.0); + xmax = gabscor->GetX()[gabscor->GetN()-1]; + } + + if (npoints<3 && !delphes) { + gabscor->SetPoint (0, 10.0,1.0); + gabscor->SetPointError(0, 0.0,0.0); + gabscor->SetPoint (1,100.0,1.0); + gabscor->SetPointError(1, 0.0,0.0); + fabscor = new TF1("fit","[0]",10.0,100.0); + fabscor->FixParameter(0,1.0); + } + else if (!l2l3 && npoints > 2 && gabscor->GetN()<10) { + fabscor=new TF1("fit","[0]+[1]*log10(x)+[2]*pow(log10(x),2)",xmin,xmax); + fabscor->SetParameter(0,1.0); + fabscor->SetParameter(1,0.0); + fabscor->SetParameter(2,0.0); + } + else { + if (alg.find("pf")!=string::npos || alg.find("puppi")!=string::npos) { + // + // Delphes + // + if(delphes) { + //fcn = "[0]+[1]*log10(x)+[2]*pow(log10(x),2)+[3]*pow(log10(x),3)+[4]*pow(log10(x),4)+[5]*pow(log10(x),5)+[6]*pow(log10(x),6)+[7]*pow(log10(x),7)+[8]*pow(log10(x),8)+[9]*pow(log10(x),9)"; + TString fcn = "[0]+[1]*log10(x)+[2]*pow(log10(x),2)+([3]/pow(log10(x),3))+([4]/pow(log10(x),4))+([5]/pow(log10(x),5))"; + fabscor=new TF1("fit",fcn.Data(),xmin,xmax); + } + + // + // online (HLT) + // + if(alg.find("HLT")!=string::npos){ + fabscor=new TF1("fit","(x>=[6])*([0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5])))+(x<[6])*[7]",xmin,xmax); + fabscor->FixParameter(6,xmin); + fabscor->FixParameter(7,0.0); + } + // + // offline + // + else { + if (xmin<10) xmin=10; + TString fcn = getOfflinePFFunction(); + fabscor=new TF1("fit",fcn.Data(),xmin,xmax); + setOfflinePFParameters(gabscor, fabscor,xmin); + } + } + else if (alg.find("trk")!=string::npos) { + fabscor=new TF1("fit","[0]+[1]*pow(x/500.0,[2])+[3]/log10(x)+[4]*log10(x)",xmin,xmax); + fabscor->SetParameter(0,1.7); + fabscor->SetParameter(1,0.7); + fabscor->SetParameter(2,3.0); + fabscor->SetParLimits(2,1,10); + fabscor->SetParameter(3,0.0); + fabscor->SetParameter(4,0.0); + } + else if (alg.find("jpt")!=string::npos || alg.find("tau")!=string::npos) { + fabscor=new TF1("fit","[0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5]))",xmin,xmax); + // + // INITIAL PARAMS: These are the fitted parameters that work for JetEta0.261to0.348 + // + fabscor->SetParameter(0,-8.9); + fabscor->SetParameter(1,530); + fabscor->SetParameter(2,16); + fabscor->SetParameter(3,-22); + fabscor->SetParameter(4,0.06); + fabscor->SetParameter(5,-0.28); + } + else if (alg.find("calo")!=string::npos) { + if(l2calofit.EqualTo("standard",TString::kIgnoreCase)) { + if (xmin<6) xmin=6; + fabscor=new TF1("fit","[0]+[1]/(pow(log10(x),[2])+[3])",xmin,xmax); + fabscor->SetParameter(0,1.0); + fabscor->SetParameter(1,5.0); + fabscor->SetParameter(2,3.0); + fabscor->SetParameter(3,3.0); + + fabscor->SetParLimits(3,0,100); + } + else if(l2calofit.EqualTo("DynamicMin",TString::kIgnoreCase)) { + xmin=findPeak(gabscor,0,3,3,false); + TString fcn = "[0] + ([1]/(pow(log10(x),2)+[2])) + ([3]*exp(-([4]*((log10(x)-[5])*(log10(x)-[5]))))) + ([6]*exp(-([7]*((log10(x)-[8])*(log10(x)-[8])))))"; + fabscor=new TF1("fit",fcn.Data(),xmin,xmax); + fabscor->SetParameter(0,-0.0221278); + fabscor->SetParameter(1,119.265); + fabscor->SetParameter(2,100); + fabscor->SetParameter(3,-0.0679365); + fabscor->SetParameter(4,2.82597); + fabscor->SetParameter(5,1.8277); + fabscor->SetParameter(6,-0.0679365); + fabscor->SetParameter(7,3.82597); + fabscor->SetParameter(8,1.8277); + fabscor->SetParLimits(6,-20,10); + fabscor->SetParLimits(7,0,100); + fabscor->SetParLimits(3,-15,15); + fabscor->SetParLimits(4,0,500); + fabscor->SetParLimits(0,-2,25); + fabscor->SetParLimits(1,0,250); + } + } + else { + cout << "WARNING::Cannot determine fit function for " << alg << "." << endl; + } } // @@ -390,21 +487,33 @@ void L2Creator::loopOverEtaBins() { perform_smart_fit(gabscor,fabscor,maxFitIter); gErrorIgnoreLevel = origIgnoreLevel; if (alg.find("pf")!=string::npos) { - if (alg.find("HLT")!=string::npos) { - ((TF1*)gabscor->GetListOfFunctions()->First())->FixParameter(7,fabscor->Eval(fabscor->GetParameter(6))); - fabscor->FixParameter(7,fabscor->Eval(fabscor->GetParameter(6))); - } + if (alg.find("HLT")!=string::npos) { + ((TF1*)gabscor->GetListOfFunctions()->First())->FixParameter(7,fabscor->Eval(fabscor->GetParameter(6))); + fabscor->FixParameter(7,fabscor->Eval(fabscor->GetParameter(6))); + } } // // format the graphs // - gabscor->GetListOfFunctions()->First()->ResetBit(TF1::kNotDraw); - gabsrsp->SetMarkerStyle(20); - gabscor->SetMarkerStyle(20); + gabscor->GetListOfFunctions()->First()->ResetBit(TF1::kNotDraw); + gabsrsp->SetMarkerStyle(20); + gabscor->SetMarkerStyle(20); + + // + // save graph + // + TCanvas *canvas = new TCanvas("", "", 700, 700); + canvas->cd(); + gPad->SetLogx(); + gabscor->Draw(); + + string savePath {outputDir+"/"+gabscor->GetName()+".png"}; + canvas->Print((savePath).c_str()); + odir->cd(); - gabsrsp->Write(); - gabscor->Write(); + gabsrsp->Write(); + gabscor->Write(); ofile->Write(); } } @@ -412,13 +521,13 @@ void L2Creator::loopOverEtaBins() { //______________________________________________________________________________ bool L2Creator::getL3Rsp() { - string alg = ji->getAbbreviation(); + string alg = ji->getAbbreviation(); l3dir = (TDirectoryFile*)l3file->Get(alg.c_str()); if (l3dir==0) { cout<<"Failed to rerieve L3 correction for "<Get("L3RspVsRefPt"); fl3rsp = (TF1*)gl3rsp->GetListOfFunctions()->First(); if (0==fl3rsp) { @@ -431,372 +540,374 @@ bool L2Creator::getL3Rsp() { //______________________________________________________________________________ void L2Creator::setAndFitFLogAndFGaus(TGraphErrors* gabscor, TF1* flog, TF1* fgaus, double xmin) { - flog = new TF1("flog", "[3]+TMath::LogNormal(TMath::Log10(x), [0], [1], [2])", 3, 20); - flog->SetLineColor(kBlue); - flog->SetRange(xmin,20.0); - flog->SetParameters(0.5, 0, 1); - flog->SetParLimits(0,0.00000001,1000); - flog->SetParLimits(1,0,1000); - flog->SetParLimits(2,0.00000001,1000); - flog->FixParameter(1,0.0); + flog = new TF1("flog", "[3]+TMath::LogNormal(TMath::Log10(x), [0], [1], [2])", 3, 20); + flog->SetLineColor(kBlue); + flog->SetRange(xmin,20.0); + flog->SetParameters(0.5, 0, 1); + flog->SetParLimits(0,0.00000001,1000); + flog->SetParLimits(1,0,1000); + flog->SetParLimits(2,0.00000001,1000); + flog->FixParameter(1,0.0); gabscor->Fit("flog","RBQS"); - fgaus = new TF1("fgaus","[0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5]))",20,4000); - fgaus->SetLineColor(kOrange); - fgaus->SetRange(20,4000); - fgaus->SetParameter(0,0.5); - fgaus->SetParameter(1,9.0); - fgaus->SetParameter(2,8.0); - fgaus->SetParameter(3,-0.3); - fgaus->SetParameter(4,0.6); - fgaus->SetParameter(5,1.0); - fgaus->SetParLimits(2,0.1,100); - fgaus->SetParLimits(3,-100,0); - fgaus->SetParLimits(4,0,100); + fgaus = new TF1("fgaus","[0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5]))",20,4000); + fgaus->SetLineColor(kOrange); + fgaus->SetRange(20,4000); + fgaus->SetParameter(0,0.5); + fgaus->SetParameter(1,9.0); + fgaus->SetParameter(2,8.0); + fgaus->SetParameter(3,-0.3); + fgaus->SetParameter(4,0.6); + fgaus->SetParameter(5,1.0); + fgaus->SetParLimits(2,0.1,100); + fgaus->SetParLimits(3,-100,0); + fgaus->SetParLimits(4,0,100); gabscor->Fit("fgaus","RQS"); } //______________________________________________________________________________ TString L2Creator::getOfflinePFFunction() { - if(l2pffit.EqualTo("standard",TString::kIgnoreCase)) { - return "[0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5]))"; - } - else if(l2pffit.EqualTo("standard+Gaussian",TString::kIgnoreCase)) { - return "[0] + ([1]/(pow(log10(x),2)+[2])) + ([3]*exp(-([4]*((log10(x)-[5])*(log10(x)-[5]))))) + ([6]*exp(-([7]*((log10(x)-[8])*(log10(x)-[8])))))"; - } - else if(l2pffit.EqualTo("LogNormal+Gaussian+fixed",TString::kIgnoreCase)) { - return "([0]+TMath::LogNormal(TMath::Log10(x),[1],[2],[3]))+([4]+[5]/(pow(log10(x),2)+[6])+[7]*exp(-[8]*(log10(x)-[9])*(log10(x)-[9])))"; - } - else if(l2pffit.EqualTo("LogNormal+Gaussian",TString::kIgnoreCase)) { - return "(x<=[10])*([0]+TMath::LogNormal(TMath::Log10(x), [1], [2], [3]))+(x>[10])*([4]+[5]/(pow(log10(x),2)+[6])+[7]*exp(-[8]*(log10(x)-[9])*(log10(x)-[9])))"; - //return "(x<=[10])*([0]+((1.0/(TMath::Log10(x)*TMath::Sqrt(2.0*TMath::Pi()*TMath::Power([2],2))))*TMath::Exp(-TMath::Power(TMath::Log(TMath::Log10(x))-[1],2)/(2.0*TMath::Power([2],2)))))+(x>[10])*([4]+[5]/(pow(log10(x),2)+[6])+[7]*exp(-[8]*(log10(x)-[9])*(log10(x)-[9])))"; - } - else { - cout << "ERROR::getOfflinePFFunction::Unknown PF function choice." << endl; - return ""; - } + if(l2pffit.EqualTo("standard",TString::kIgnoreCase)) { + return "[0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5]))"; + } + else if(l2pffit.EqualTo("standard+Gaussian",TString::kIgnoreCase)) { + return "[0] + ([1]/(pow(log10(x),2)+[2])) + ([3]*exp(-([4]*((log10(x)-[5])*(log10(x)-[5]))))) + ([6]*exp(-([7]*((log10(x)-[8])*(log10(x)-[8])))))"; + } + else if(l2pffit.EqualTo("LogNormal+Gaussian+fixed",TString::kIgnoreCase)) { + return "([0]+TMath::LogNormal(TMath::Log10(x),[1],[2],[3]))+([4]+[5]/(pow(log10(x),2)+[6])+[7]*exp(-[8]*(log10(x)-[9])*(log10(x)-[9])))"; + } + else if(l2pffit.EqualTo("LogNormal+Gaussian",TString::kIgnoreCase)) { + return "(x<=[10])*([0]+TMath::LogNormal(TMath::Log10(x), [1], [2], [3]))+(x>[10])*([4]+[5]/(pow(log10(x),2)+[6])+[7]*exp(-[8]*(log10(x)-[9])*(log10(x)-[9])))"; + //return "(x<=[10])*([0]+((1.0/(TMath::Log10(x)*TMath::Sqrt(2.0*TMath::Pi()*TMath::Power([2],2))))*TMath::Exp(-TMath::Power(TMath::Log(TMath::Log10(x))-[1],2)/(2.0*TMath::Power([2],2)))))+(x>[10])*([4]+[5]/(pow(log10(x),2)+[6])+[7]*exp(-[8]*(log10(x)-[9])*(log10(x)-[9])))"; + } + else { + cout << "ERROR::getOfflinePFFunction::Unknown PF function choice." << endl; + return ""; + } } //______________________________________________________________________________ void L2Creator::setOfflinePFParameters(TGraphErrors* gabscor, TF1* fabscor, double xmin) { - TF1* flog(0); - TF1* fgaus(0); - if(l2pffit.Contains("LogNormal+Gaussian",TString::kIgnoreCase)) { - //int origIgnoreLevel = gErrorIgnoreLevel; - //gErrorIgnoreLevel = kBreak; - setAndFitFLogAndFGaus(gabscor, flog, fgaus, xmin); - //gErrorIgnoreLevel = origIgnoreLevel; - //Original function - } + TF1* flog(0); + TF1* fgaus(0); + if(l2pffit.Contains("LogNormal+Gaussian",TString::kIgnoreCase)) { + //int origIgnoreLevel = gErrorIgnoreLevel; + //gErrorIgnoreLevel = kBreak; + setAndFitFLogAndFGaus(gabscor, flog, fgaus, xmin); + //gErrorIgnoreLevel = origIgnoreLevel; + //Original function + } - if(l2pffit.EqualTo("standard",TString::kIgnoreCase)) { - fabscor->SetParameter(0,0.5); - fabscor->SetParameter(1,9.0); - fabscor->SetParameter(2,8.0); - fabscor->SetParameter(3,-0.3); - fabscor->SetParameter(4,0.6); - fabscor->SetParameter(5,1.0); - fabscor->SetParLimits(2,0.1,100); - fabscor->SetParLimits(3,-100,0); - fabscor->SetParLimits(4,0,100); - } - else if(l2pffit.EqualTo("standard+Gaussian",TString::kIgnoreCase)) { - fabscor->SetParameter(0,-0.0221278); - fabscor->SetParameter(1,119.265); - fabscor->SetParameter(2,100); - fabscor->SetParameter(3,-0.0679365); - fabscor->SetParameter(4,2.82597); - fabscor->SetParameter(5,1.8277); - fabscor->SetParameter(6,-0.0679365); - fabscor->SetParameter(7,3.82597); - fabscor->SetParameter(8,1.8277); - fabscor->SetParLimits(6,-20,10); - fabscor->SetParLimits(7,0,100); - fabscor->SetParLimits(3,-15,15); - fabscor->SetParLimits(4,0,500); - fabscor->SetParLimits(0,-2,25); - fabscor->SetParLimits(1,0,250); - } - else if(l2pffit.EqualTo("LogNormal+Gaussian+fixed",TString::kIgnoreCase)) { - fabscor->SetRange(3,2000); - fabscor->SetParameter(0,flog->GetParameter(3)); - fabscor->SetParameter(1,flog->GetParameter(0)); - fabscor->SetParameter(2,flog->GetParameter(1)); - fabscor->SetParameter(3,flog->GetParameter(2)); - fabscor->SetParLimits(1,0,5); - fabscor->SetParLimits(2,-20,20); - fabscor->SetParLimits(3,0,5); - fabscor->SetParameter(4,fgaus->GetParameter(0)); - fabscor->SetParameter(5,fgaus->GetParameter(1)); - fabscor->SetParameter(6,fgaus->GetParameter(2)); - fabscor->SetParameter(7,fgaus->GetParameter(3)); - fabscor->SetParameter(8,fgaus->GetParameter(4)); - fabscor->SetParameter(9,fgaus->GetParameter(5)); - fabscor->SetParLimits(6,0.1,100); - fabscor->SetParLimits(7,-100,0); - fabscor->SetParLimits(8,0,100); - } - else if(l2pffit.EqualTo("LogNormal+Gaussian",TString::kIgnoreCase)) { - fabscor->SetRange(xmin,4000); - fabscor->SetParameter(0,flog->GetParameter(3)); - fabscor->SetParameter(1,flog->GetParameter(0)); - fabscor->SetParameter(2,flog->GetParameter(1)); - fabscor->SetParameter(3,flog->GetParameter(2)); - fabscor->SetParLimits(1,0.00000001,1000); - fabscor->SetParLimits(2,0,1000); - fabscor->SetParLimits(3,0.00000001,1000); - fabscor->FixParameter(2,0.0); - fabscor->SetParameter(4,fgaus->GetParameter(0)); - fabscor->SetParameter(5,fgaus->GetParameter(1)); - fabscor->SetParameter(6,fgaus->GetParameter(2)); - fabscor->SetParameter(7,fgaus->GetParameter(3)); - fabscor->SetParameter(8,fgaus->GetParameter(4)); - fabscor->SetParameter(9,fgaus->GetParameter(5)); - fabscor->SetParLimits(6,0.1,100); - fabscor->SetParLimits(7,-100,0); - fabscor->SetParLimits(8,0,100); - fabscor->SetParameter(10,20); - } + if(l2pffit.EqualTo("standard",TString::kIgnoreCase)) { + fabscor->SetParameter(0,0.5); + fabscor->SetParameter(1,9.0); + fabscor->SetParameter(2,8.0); + fabscor->SetParameter(3,-0.3); + fabscor->SetParameter(4,0.6); + fabscor->SetParameter(5,1.0); + fabscor->SetParLimits(2,0.1,100); + fabscor->SetParLimits(3,-100,0); + fabscor->SetParLimits(4,0,100); + } + else if(l2pffit.EqualTo("standard+Gaussian",TString::kIgnoreCase)) { + fabscor->SetParameter(0,-0.0221278); + fabscor->SetParameter(1,119.265); + fabscor->SetParameter(2,100); + fabscor->SetParameter(3,-0.0679365); + fabscor->SetParameter(4,2.82597); + fabscor->SetParameter(5,1.8277); + fabscor->SetParameter(6,-0.0679365); + fabscor->SetParameter(7,3.82597); + fabscor->SetParameter(8,1.8277); + fabscor->SetParLimits(6,-20,10); + fabscor->SetParLimits(7,0,100); + fabscor->SetParLimits(3,-15,15); + fabscor->SetParLimits(4,0,500); + fabscor->SetParLimits(0,-2,25); + fabscor->SetParLimits(1,0,250); + } + else if(l2pffit.EqualTo("LogNormal+Gaussian+fixed",TString::kIgnoreCase)) { + fabscor->SetRange(3,2000); + fabscor->SetParameter(0,flog->GetParameter(3)); + fabscor->SetParameter(1,flog->GetParameter(0)); + fabscor->SetParameter(2,flog->GetParameter(1)); + fabscor->SetParameter(3,flog->GetParameter(2)); + fabscor->SetParLimits(1,0,5); + fabscor->SetParLimits(2,-20,20); + fabscor->SetParLimits(3,0,5); + fabscor->SetParameter(4,fgaus->GetParameter(0)); + fabscor->SetParameter(5,fgaus->GetParameter(1)); + fabscor->SetParameter(6,fgaus->GetParameter(2)); + fabscor->SetParameter(7,fgaus->GetParameter(3)); + fabscor->SetParameter(8,fgaus->GetParameter(4)); + fabscor->SetParameter(9,fgaus->GetParameter(5)); + fabscor->SetParLimits(6,0.1,100); + fabscor->SetParLimits(7,-100,0); + fabscor->SetParLimits(8,0,100); + } + else if(l2pffit.EqualTo("LogNormal+Gaussian",TString::kIgnoreCase)) { + fabscor->SetRange(xmin,4000); + fabscor->SetParameter(0,flog->GetParameter(3)); + fabscor->SetParameter(1,flog->GetParameter(0)); + fabscor->SetParameter(2,flog->GetParameter(1)); + fabscor->SetParameter(3,flog->GetParameter(2)); + fabscor->SetParLimits(1,0.00000001,1000); + fabscor->SetParLimits(2,0,1000); + fabscor->SetParLimits(3,0.00000001,1000); + fabscor->FixParameter(2,0.0); + fabscor->SetParameter(4,fgaus->GetParameter(0)); + fabscor->SetParameter(5,fgaus->GetParameter(1)); + fabscor->SetParameter(6,fgaus->GetParameter(2)); + fabscor->SetParameter(7,fgaus->GetParameter(3)); + fabscor->SetParameter(8,fgaus->GetParameter(4)); + fabscor->SetParameter(9,fgaus->GetParameter(5)); + fabscor->SetParLimits(6,0.1,100); + fabscor->SetParLimits(7,-100,0); + fabscor->SetParLimits(8,0,100); + fabscor->SetParameter(10,20); + } } //______________________________________________________________________________ Double_t L2Creator::findPeak(TGraphErrors* gabscor, int ipeak, int npeaks, int res, bool verbose) { - //Make a TH1F to find the peak - int size = gabscor->GetN(); - map > points; - Float_t* xlow = new Float_t[size]; - Float_t* y = new Float_t[size]; - Float_t* ye = new Float_t[size]; - for(int ipoint=0; ipointGetX()[ipoint]] = make_pair(gabscor->GetY()[ipoint],gabscor->GetEY()[ipoint]); - } - int ipoint=0; - for (std::map >::iterator it=points.begin(); it!=points.end(); ++it) { - xlow[ipoint] = it->first; - y[ipoint] = it->second.first; - ye[ipoint] = it->second.second; - ipoint++; - } - TH1F* hpeak = new TH1F("hpeak","hpeak",gabscor->GetN()-1,xlow); - for(int ibin=0; ibinGetNbinsX(); ibin++) { - hpeak->SetBinContent(ibin,y[ibin]); - hpeak->SetBinError(ibin,ye[ibin]); - } - TSpectrum *spec = new TSpectrum(npeaks); - spec->Search(hpeak,res,"nobackground nodraw goff"); - Double_t* xpos = spec->GetPositionX(); - Double_t* ypos = spec->GetPositionY(); - if(verbose) { - std::cout << "peak: " << xpos[ipeak] << std::endl; - std::cout << "peak height: " << ypos[ipeak] << std::endl; - } - delete hpeak; - return xpos[ipeak]; + //Make a TH1F to find the peak + int size = gabscor->GetN(); + map > points; + Float_t* xlow = new Float_t[size]; + Float_t* y = new Float_t[size]; + Float_t* ye = new Float_t[size]; + for(int ipoint=0; ipointGetX()[ipoint]] = make_pair(gabscor->GetY()[ipoint],gabscor->GetEY()[ipoint]); + } + int ipoint=0; + for (std::map >::iterator it=points.begin(); it!=points.end(); ++it) { + xlow[ipoint] = it->first; + y[ipoint] = it->second.first; + ye[ipoint] = it->second.second; + ipoint++; + } + TH1F* hpeak = new TH1F("hpeak","hpeak",gabscor->GetN()-1,xlow); + for(int ibin=0; ibinGetNbinsX(); ibin++) { + hpeak->SetBinContent(ibin,y[ibin]); + hpeak->SetBinError(ibin,ye[ibin]); + } + TSpectrum *spec = new TSpectrum(npeaks); + spec->Search(hpeak,res,"nobackground nodraw goff"); + Double_t* xpos = spec->GetPositionX(); + Double_t* ypos = spec->GetPositionY(); + if(verbose) { + std::cout << "peak: " << xpos[ipeak] << std::endl; + std::cout << "peak height: " << ypos[ipeak] << std::endl; + } + delete hpeak; + return xpos[ipeak]; } //______________________________________________________________________________ void L2Creator::doRelCorFits() { - string fnc_as_str = (ji->getAbbreviation().find("trk")!=string::npos) ? - "[0]+[1]*log10(x)+[2]*pow(log10(x),2)+[3]*pow(log10(x),3)+[4]*pow(x/500.0,3)" : - "[0]+[1]*log10(x)+[2]*pow(log10(x),2)+[3]*pow(log10(x),3)+[4]*pow(log10(x),4)"; - - vector indices; - TH1F* hjetpt(0); - hl_jetpt.begin_loop(); - while ((hjetpt=hl_jetpt.next_object(indices))) { - - unsigned int ieta = indices[0]; - unsigned int ipt = indices[1]; - - // create a new graph if a new eta bin comes around - if (ipt==0) { - vrelcor_eta.push_back(new TGraphErrors()); - stringstream ss; - ss<SetName(("RelCorVsJetPt_JetEta"+ss.str()).c_str()); - } - - // only add a point to the graph if the current histo is not empty - if (hjetpt->Integral()!=0) { - TF1* fabscor =vabscor_eta[ieta]->GetFunction("fit"); - double jetpt =hjetpt->GetMean(); - if(!fabscor) continue; - double refpt =jetpt*fabscor->Eval(jetpt); - double l3cor = 1; - if (!l2l3) l3cor = fl3rsp->Eval(refpt); - double controlpt=refpt*l3cor; - double relcor =controlpt/jetpt; - if (relcor > 5) - cout<<"WARNING !!! suspicious point: "<GetName() - <<", jet pt = "<GetN(); - vrelcor_eta.back()->SetPoint(n,jetpt,relcor); - } - } - - // fit the graph if the last pt of the current eta bin comes around - if (ipt==hl_jetpt.nobjects(1)-1 && (vrelcor_eta.back())->GetN()!=0) { - TGraph* grelcor = vrelcor_eta.back(); - double xmin = grelcor->GetX()[0]; - double xmax = grelcor->GetX()[grelcor->GetN()-1]; - TF1* frelcor = new TF1("fit",fnc_as_str.c_str(),xmin,xmax); - - frelcor->SetParameter(0,0.0); - frelcor->SetParameter(1,0.0); - frelcor->SetParameter(2,0.0); - frelcor->SetParameter(3,0.0); - frelcor->SetParameter(4,0.0); - frelcor->SetParameter(5,0.0); - - if (grelcor->GetN()<2) { - grelcor->SetPoint(0,10,1.0); - grelcor->SetPoint(1,100,1.0); - frelcor->FixParameter(1,0.0); - frelcor->FixParameter(2,0.0); - frelcor->FixParameter(3,0.0); - frelcor->FixParameter(4,0.0); - frelcor->FixParameter(5,0.0); - } - else if (grelcor->GetN()==2) { - frelcor->FixParameter(2,0.0); - frelcor->FixParameter(3,0.0); - frelcor->FixParameter(4,0.0); - frelcor->FixParameter(5,0.0); - } - - grelcor->Fit(frelcor,"QRB0"); - grelcor->GetListOfFunctions()->First()->ResetBit(TF1::kNotDraw); - grelcor->SetMarkerStyle(20); - grelcor->Write(); - } + string fnc_as_str = (ji->getAbbreviation().find("trk")!=string::npos) ? + "[0]+[1]*log10(x)+[2]*pow(log10(x),2)+[3]*pow(log10(x),3)+[4]*pow(x/500.0,3)" : + "[0]+[1]*log10(x)+[2]*pow(log10(x),2)+[3]*pow(log10(x),3)+[4]*pow(log10(x),4)"; + + vector indices; + TH1F* hjetpt(0); + hl_jetpt.begin_loop(); + while ((hjetpt=hl_jetpt.next_object(indices))) { + + unsigned int ieta = indices[0]; + unsigned int ipt = indices[1]; + + // create a new graph if a new eta bin comes around + if (ipt==0) { + vrelcor_eta.push_back(new TGraphErrors()); + stringstream ss; + ss<SetName(("RelCorVsJetPt_JetEta"+ss.str()).c_str()); } + + // only add a point to the graph if the current histo is not empty + if (hjetpt->Integral()!=0) { + TF1* fabscor =vabscor_eta[ieta]->GetFunction("fit"); + double jetpt =hjetpt->GetMean(); + if(!fabscor) continue; + double refpt =jetpt*fabscor->Eval(jetpt); + double l3cor = 1; + if (!l2l3) l3cor = fl3rsp->Eval(refpt); + double controlpt=refpt*l3cor; + double relcor =controlpt/jetpt; + if (relcor > 5) + cout<<"WARNING !!! suspicious point: "<GetName() + <<", jet pt = "<GetN(); + vrelcor_eta.back()->SetPoint(n,jetpt,relcor); + } + } + + // fit the graph if the last pt of the current eta bin comes around + if (ipt==hl_jetpt.nobjects(1)-1 && (vrelcor_eta.back())->GetN()!=0) { + TGraph* grelcor = vrelcor_eta.back(); + double xmin = grelcor->GetX()[0]; + double xmax = grelcor->GetX()[grelcor->GetN()-1]; + TF1* frelcor = new TF1("fit",fnc_as_str.c_str(),xmin,xmax); + + frelcor->SetParameter(0,0.0); + frelcor->SetParameter(1,0.0); + frelcor->SetParameter(2,0.0); + frelcor->SetParameter(3,0.0); + frelcor->SetParameter(4,0.0); + frelcor->SetParameter(5,0.0); + + if (grelcor->GetN()<2) { + grelcor->SetPoint(0,10,1.0); + grelcor->SetPoint(1,100,1.0); + frelcor->FixParameter(1,0.0); + frelcor->FixParameter(2,0.0); + frelcor->FixParameter(3,0.0); + frelcor->FixParameter(4,0.0); + frelcor->FixParameter(5,0.0); + } + else if (grelcor->GetN()==2) { + frelcor->FixParameter(2,0.0); + frelcor->FixParameter(3,0.0); + frelcor->FixParameter(4,0.0); + frelcor->FixParameter(5,0.0); + } + + grelcor->Fit(frelcor,"QRB0"); + grelcor->GetListOfFunctions()->First()->ResetBit(TF1::kNotDraw); + grelcor->SetMarkerStyle(20); + grelcor->Write(); + } + } } //______________________________________________________________________________ bool L2Creator::contains(const vector& collection,const string& element) { - vector::const_iterator it; - for (it=collection.begin();it!=collection.end();++it) - if ((*it)==element) return true; - return false; + vector::const_iterator it; + for (it=collection.begin();it!=collection.end();++it) + if ((*it)==element) return true; + return false; } //______________________________________________________________________________ void L2Creator::perform_smart_fit(TGraphErrors * gabscor, TF1 * fabscor, int maxFitIter) { - int fitIter = 0; - vector bestPars; - double bestRChi2 = 0; - do { - // - // do the fit, get the results and the parameters of the fitted function - // - TFitResultPtr fitResPtr = gabscor->Fit(fabscor,"RQS+"); - vector auxPars = fitResPtr.Get()->Parameters(); - - // - // compute the reduced chi2 of this fit and if it is the best fit so far - // then save the parameters - // - double rchi2 = fitResPtr.Get()->Chi2()/ fitResPtr.Get()->Ndf(); - if (fitResPtr.Get()->Ndf() == 0) rchi2 = 0; - if (rchi2 > 0 && (rchi2 2 || bestRChi2 == 0 ) && fitIter < maxFitIter); - - // - // set the best parameters and chi2 to the fit function - // - TF1 * ffh = gabscor->GetFunction("fit"); - for (unsigned int np=0;np < bestPars.size() ; np++){ - ffh->SetParameter(np,bestPars[np]); - fabscor->SetParameter(np,bestPars[np]); - } - fabscor->SetChisquare(bestRChi2 * fabscor->GetNDF()); - ffh->SetChisquare(bestRChi2 * fabscor->GetNDF()); - - // - // warn if the fit diverges at low pt - // - if (fabscor->Integral(0,10) > 60) - cout << "\t***ERROR***, fit for histo " << gabscor->GetName() << " diverges at low pt" << endl; - - // - // check for failed fits - // a chi2 of zero is symptomatic of a failed fit. - // - if (bestRChi2 < 0.001) { - cout<<"\t***ERROR***, FIT HAS FAILED for histo "<GetName() - <<" which has a reduced chi2="< 5){ - if (bestRChi2 > 10) - cout<<"\t***ERROR***,"; - else - cout<<"\tWARNING,"; - - cout<<" fit for histo "<GetName() - <<" has a reduced chi2="< bestPars; + double bestRChi2 = 0; + do { + // + // do the fit, get the results and the parameters of the fitted function + // + //TFitResultPtr fitResPtr = gabscor->Fit(fabscor,"RQS+"); + TFitResultPtr fitResPtr = gabscor->Fit(fabscor,"RQS"); + vector auxPars = fitResPtr.Get()->Parameters(); + + // + // compute the reduced chi2 of this fit and if it is the best fit so far + // then save the parameters + // + double rchi2 = fitResPtr.Get()->Chi2()/ fitResPtr.Get()->Ndf(); + if (fitResPtr.Get()->Ndf() == 0) rchi2 = 0; + if (rchi2 > 0 && (rchi2 2 || bestRChi2 == 0 ) && fitIter < maxFitIter); + + // + // set the best parameters and chi2 to the fit function + // + TF1 * ffh = gabscor->GetFunction("fit"); + for (unsigned int np=0;np < bestPars.size() ; np++){ + ffh->SetParameter(np,bestPars[np]); + fabscor->SetParameter(np,bestPars[np]); + } + fabscor->SetChisquare(bestRChi2 * fabscor->GetNDF()); + ffh->SetChisquare(bestRChi2 * fabscor->GetNDF()); + + // + // warn if the fit diverges at low pt + // + if (fabscor->Integral(0,10) > 60) + cout << "\t***ERROR***, fit for histo " << gabscor->GetName() << " diverges at low pt" << endl; + + // + // check for failed fits + // a chi2 of zero is symptomatic of a failed fit. + // + if (bestRChi2 < 0.001) { + cout<<"\t***ERROR***, FIT HAS FAILED for histo "<GetName() + <<" which has a reduced chi2="< 5){ + if (bestRChi2 > 10) + cout<<"\t***ERROR***,"; + else + cout<<"\tWARNING,"; + + cout<<" fit for histo "<GetName() + <<" has a reduced chi2="<getAlias() + ".txt"; - ofstream fout(txtfilename); - fout.setf(ios::right); - - unsigned int vector_size = 0; - if(l2l3) vector_size = vabscor_eta.size(); //For L2L3 Corrections Together - else vector_size = vrelcor_eta.size(); //For L2 & L3 Corrections Separate - for (unsigned int ieta=0;ietaGetListOfFunctions()->Last(); - if(frelcor!=0) { - if(ieta==0 || (ieta==1 && delphes)) - fout<<"{1 JetEta 1 JetPt max(0.0001,"<GetTitle()<<") Correction L2Relative}"<GetExpFormula()<<") Correction L2Relative}"<GetX()[0]; - double ptmax = grelcor->GetX()[grelcor->GetN()-1]; - fout<GetNpar()+2) //Number of parameters + 2 - <GetNpar(); p++) { - fout<GetParameter(p); - } - fout<getAlias() + ".txt"; + ofstream fout(txtfilename); + fout.setf(ios::right); + + unsigned int vector_size = 0; + if(l2l3) vector_size = vabscor_eta.size(); //For L2L3 Corrections Together + else vector_size = vrelcor_eta.size(); //For L2 & L3 Corrections Separate + for (unsigned int ieta=0;ietaGetListOfFunctions()->Last(); + if(frelcor!=0) { + if(ieta==0 || (ieta==1 && delphes)){ + if(cbPar==17) fout<<"{1 JetEta 1 JetPt max(-1.,"<GetTitle()<<") Correction L3Absolute}"<GetTitle()<<") Correction L3Absolute}"<GetX()[0]; + double ptmax = grelcor->GetX()[grelcor->GetN()-1]; + fout<GetNpar()+2) //Number of parameters + 2 + <GetNpar(); p++) { + fout<GetParameter(p); + } + fout<Write(); - ofile->Close(); - delete ofile; - ifile->Close(); - delete ifile; - cout<<"DONE"<Write(); + ofile->Close(); + delete ofile; + ifile->Close(); + delete ifile; + cout<<"DONE"<= 0. && x < R1 ) result += pol_fcn(x, p0_left, p1_left, p2_left, p3_left, p4_left, p5_left, p6_left, p7_left); + if (drawExp && x >= R1 && x < R2 ) result += exp_fcn(x, normExp_left, slopeExp_left); + if (drawGaus1 && x >= R2 && x < meanGaus1) result += gauss_fcn(x, normGaus1, meanGaus1, sigGaus1_left); + if (drawAnyRight && x >= meanGaus1 && x < 2. ) result += sum_gaus_pol_fnc(xx, pp); + return result; +} // crystalBall + + +double crystalBall_1(double *xx, double *pp) +{ + double x = xx[0]; + + // left + double normGaus1 = pp[0]; + double meanGaus1 = pp[1]; + double sigGaus1_left = pp[2]; + + double slopeExp_left = pp[4]; + + double p1_left = pp[6]; + + double p1_right = pp[17]; + + double R2 = meanGaus1 - 0.5*sigGaus1_left - pp[8]; + double R1 = R2 - pp[7]; // min exp + double R3 = pp[18]; + + bool drawPolLeft = pp[19]; + bool drawExp = pp[20]; + bool drawGaus1 = pp[21]; + bool drawGaus2 = pp[22]; + bool drawGaus3 = pp[23]; + bool drawPolRight = pp[24]; + bool drawAnyRight = drawGaus1 || drawGaus2 || drawGaus3; + + // remove discontinuity at x=R2: solve exp=gauss for norm exp and replace R2 by (meanGaus1-R2*R2) + double normExp_left = gauss_fcn(R2, normGaus1, meanGaus1, sigGaus1_left)/TMath::Exp(slopeExp_left*R2); + + // remove discontinuity at x=R1: solve pol1=exp for p0 left + double p0_left = exp_fcn(R1, normExp_left, slopeExp_left) - p1_left*R1 ; + + // remove discontinuity at x=R3: solve pol1=sumGaus for p0 right + double r3[]{R3}; + double p0_right = sum_gaus_fnc(r3, pp)-p1_right*R3; + + // compute probability + double result = 0.; + if (drawPolLeft && x >= 0. && x < R1 ) result += pol_fcn (x, p0_left, p1_left); + if (drawExp && x >= R1 && x < R2 ) result += exp_fcn (x, normExp_left, slopeExp_left); + if (drawGaus1 && x >= R2 && x < meanGaus1) result += gauss_fcn(x, normGaus1, meanGaus1, sigGaus1_left); + if (drawAnyRight && x >= meanGaus1 && x < R3 ) result += sum_gaus_fnc(xx, pp); + if (drawPolRight && x >= R3 && x < 2. ) result += pol_fcn (x, p0_right, p1_right); + + return result; +} // crystalBall_1 + + + +////////////////////// +double exp_fcn(double x, double norm, double slope){ + return norm*TMath::Exp(slope*x); +} + +double gauss_fcn(double x, double norm, double mean, double sigma){ + return norm*TMath::Exp(-0.5*TMath::Power((x - mean)/sigma, 2.)); +} + +double pol_fcn(double x, double p0, double p1, double p2, double p3, double p4, double p5, double p6, double p7){ + return p0 + p1*x + p2*TMath::Power(x,2) + + p3*TMath::Power(x,3) + p4*TMath::Power(x,4) + + p5*TMath::Power(x,5) + p6*TMath::Power(x,6) + + p7*TMath::Power(x,7); +} + +double sum_gaus_pol_fnc(double *xx, double *pp) +{ + double x = xx[0]; + + double normGaus1Left = pp[0]; + double meanGaus1 = pp[1]; + double valGaus1_at_meanGaus1 = normGaus1Left; + + double normGaus2 = pp[16]; + double meanGaus2 = meanGaus1 + pp[17]; + double sigGaus2 = pp[18]; + double valGaus2_at_meanGaus1 = gauss_fcn(meanGaus1, normGaus2, meanGaus2, sigGaus2); + double valGaus2 = gauss_fcn(x, normGaus2, meanGaus2, sigGaus2); + + double normGaus3 = pp[19]; + double meanGaus3 = meanGaus2 + pp[20]; + double sigGaus3 = pp[21]; + double valGaus3_at_meanGaus1 = gauss_fcn(meanGaus1, normGaus3, meanGaus3, sigGaus3); + double valGaus3 = gauss_fcn(x, normGaus3, meanGaus3, sigGaus3); + + double p0 = pp[22]; + double p1 = pp[23]; + double p2 = pp[24]; + double p3 = pp[25]; + double p4 = pp[26]; + double p5 = pp[27]; + double p6 = pp[28]; + double p7 = pp[29]; + double valPol_at_meanGaus1 = pol_fcn(meanGaus1, p0, p1, p2, p3, p4, p5, p6, p7); + double valPol = pol_fcn(x, p0, p1, p2, p3, p4, p5, p6, p7); + + double normGaus1Right = valGaus1_at_meanGaus1 - (valGaus2_at_meanGaus1 + valGaus3_at_meanGaus1 + valPol_at_meanGaus1); + double sigGaus1_right = pp[15]; + double valGaus1_right = gauss_fcn(x, normGaus1Right, meanGaus1, sigGaus1_right); + + bool drawGaus1 = pp[32]; + bool drawGaus2 = pp[33]; + bool drawGaus3 = pp[34]; + bool drawPolRight = pp[35]; + + //if ( drawGaus1 && drawGaus2 && drawGaus3 && drawPolRight && x > 0.99 && x < 1.01 ) { + // std::cout << ":" << std::endl; + // std::cout << " x = " << x << ": g1 = " << valGaus1_right << ", g2 = " << valGaus2 << ", g3 = " << valGaus3 << ", pol = " << valPol << std::endl; + //} + + double sum = 0.; + if ( drawGaus1 ) sum += valGaus1_right; + if ( drawGaus2 ) sum += valGaus2; + if ( drawGaus3 ) sum += valGaus3; + if ( drawPolRight ) sum += valPol; + if ( normGaus1Right < 0. ) sum /= (1. + normGaus1Right*normGaus1Right); + return sum; +} + + +double sum_gaus_fnc(double *xx, double *pp) +{ + double x = xx[0]; + + double normGaus1Left = pp[0]; + double meanGaus1 = pp[1]; + double valGaus1_at_meanGaus1 = normGaus1Left; + + double normGaus2 = pp[10]; + double meanGaus2 = meanGaus1 + pp[11]; + double sigGaus2 = pp[12]; + double valGaus2_at_meanGaus1 = gauss_fcn(meanGaus1, normGaus2, meanGaus2, sigGaus2); + double valGaus2 = gauss_fcn(x, normGaus2, meanGaus2, sigGaus2); + + double normGaus3 = pp[13]; + double meanGaus3 = meanGaus2 + pp[14]; + double sigGaus3 = pp[15]; + double valGaus3_at_meanGaus1 = gauss_fcn(meanGaus1, normGaus3, meanGaus3, sigGaus3); + double valGaus3 = gauss_fcn(x, normGaus3, meanGaus3, sigGaus3); + + double normGaus1Right = valGaus1_at_meanGaus1 - (valGaus2_at_meanGaus1 + valGaus3_at_meanGaus1); + double sigGaus1_right = pp[9]; + double valGaus1_right = gauss_fcn(x, normGaus1Right, meanGaus1, sigGaus1_right); + + bool drawGaus1 = pp[21]; + bool drawGaus2 = pp[22]; + bool drawGaus3 = pp[23]; + + //if ( drawGaus1 && drawGaus2 && drawGaus3 && x > 0.99 && x < 1.01 ) { + // std::cout << ":" << std::endl; + // std::cout << " x = " << x << ": g1 = " << valGaus1_right << ", g2 = " << valGaus2 << ", g3 = " << valGaus3 << std::endl; + //} + + double sum = 0.; + if ( drawGaus1 ) sum += valGaus1_right; + if ( drawGaus2 ) sum += valGaus2; + if ( drawGaus3 ) sum += valGaus3; + if ( normGaus1Right < 0. ) sum /= (1. + normGaus1Right*normGaus1Right); + return sum; +} + + + + + diff --git a/JetUtilities/src/crystalBallRes.cc b/JetUtilities/src/crystalBallRes.cc new file mode 100644 index 00000000..fe08c3ea --- /dev/null +++ b/JetUtilities/src/crystalBallRes.cc @@ -0,0 +1,213 @@ +#include "JetMETAnalysis/JetUtilities/interface/crystalBall.h" +#include "JetMETAnalysis/JetUtilities/interface/crystalBallRes.h" + +#include +#include +#include + +#include +#include + + +using namespace std; + +void crystalBallRes( TH1 *hrsp, bool drawCrystalBall, vdouble &vintegral, vdouble &vnorm, + vdouble &vrange, vdouble &vmean, vdouble &vsigma, + vdouble &vexp, vdouble &vpolLeft, vdouble &vpolRight, + vint &voption, int polDeg, TCanvas *canvas) +{ + // + // retrieve the fit + // + TF1 *crystalBallFit=(TF1*)hrsp->GetListOfFunctions()->Last(); + if(!crystalBallFit) {cout<<"crystalBallRes: pointor to CB fit failed!!"; return;} + + // + // print fit parameters + // + cout<<"\n"<GetNpar()<<" fit parameters:\n"; + for(int par=0; parGetNpar(); par++) cout<GetParName(par)<<": "<GetParameter(par)<GetParameter(1); // mean gauss1 + vrange[1] = vmean[0]-0.5*crystalBallFit->GetParameter(2)-crystalBallFit->GetParameter(14); // min gauss1 + vrange[0] = vrange[1]-crystalBallFit->GetParameter(13); // min exp + vmean[1] = vmean[0]+crystalBallFit->GetParameter(17); // mean gauss2 + vmean[2] = vmean[1]+crystalBallFit->GetParameter(20); // mean gauss3 + // + vsigma[0] = crystalBallFit->GetParameter(2); // sigma gauss1 right + vsigma[1] = crystalBallFit->GetParameter(15); // sigma gauss1 left + vsigma[2] = crystalBallFit->GetParameter(18); // sigma gauss2 + vsigma[3] = crystalBallFit->GetParameter(21); // sigma gauss3 + // + vnorm[0] = crystalBallFit->GetParameter(0); // norm gauss1 + vnorm[1] = crystalBallFit->GetParameter(16); // norm gauss2 + vnorm[2] = crystalBallFit->GetParameter(19); // norm gauss3 + // + vexp[0] = crystalBallFit->GetParameter(3); // norm exp + vexp[1] = crystalBallFit->GetParameter(4); // slope exp + // + for(int i=0; i<8; i++) vpolLeft[i] = crystalBallFit->GetParameter(i+5); + // + for(int i=0; i<8; i++) vpolRight[i] = crystalBallFit->GetParameter(i+22); + // + for(int i=0; i<6; i++) voption[i] = crystalBallFit->GetParameter(i+30); + + // will draw fit in different color for each range if option set + vbool sum { 1, 1, 1, 1, 1, 1 }; + vbool polL { 1, 0, 0, 0, 0, 0 }; + vbool exp { 0, 1, 0, 0, 0, 0 }; + vbool g1 { 0, 0, 1, 0, 0, 0 }; + vbool g2 { 0, 0, 0, 1, 0, 0 }; + vbool g3 { 0, 0, 0, 0, 1, 0 }; + vbool polR { 0, 0, 0, 0, 0, 1 }; + vector drawOpt{sum, polL, exp, g1, g2, g3, polR}; + + unsigned long nfnc{drawOpt.size()}; + + TF1 *fnc[nfnc]; + + for(unsigned i=0; iGetNpar()); + + for( int par=0; parGetNpar(); par++ ){ + fnc[i]->SetParName (par, crystalBallFit->GetParName(par)); + fnc[i]->SetParameter(par, crystalBallFit->GetParameter(par)); + } + for(unsigned j=0; jSetParameter(j+30, drawOpt[i][j]); + + //if(i==0) fnc[i]->SetLineColor(kBlack); // sum + if(i==1) fnc[i]->SetLineColor(kGray+2); // polL + if(i==2) fnc[i]->SetLineColor(kCyan); // exp + if(i==3){fnc[i]->SetLineColor(kRed); vintegral[0]=fnc[i]->Integral(0,2)/hrsp->GetBinWidth(1);} // g1 + if(i==4){fnc[i]->SetLineColor(kMagenta); vintegral[1]=fnc[i]->Integral(0,2)/hrsp->GetBinWidth(1);} // g2 + if(i==5){fnc[i]->SetLineColor(kOrange+4); vintegral[2]=fnc[i]->Integral(0,2)/hrsp->GetBinWidth(1);} // g3 + if(i==6) fnc[i]->SetLineColor(kGray+2); // polR + if(drawCrystalBall && canvas!=0){ canvas->cd(2); fnc[i]->Draw("same"); } + } + + + // + // print fit info + // + cout<<"\norder after fit: "<GetXaxis()->GetBinCenter(hrsp->GetMaximumBin())}; + cout<<"|mu1-peak|: "<GetParameter(1); // mean gauss1 + vrange[1] = vmean[0]-0.5*crystalBallFit->GetParameter(2)-crystalBallFit->GetParameter(8); // min gauss1 + vrange[0] = vrange[1]-crystalBallFit->GetParameter(7); // min exp + vrange[2] = crystalBallFit->GetParameter(18); // maxSumGaus + vmean[1] = vmean[0]+crystalBallFit->GetParameter(11); // mean gauss2 + vmean[2] = vmean[1]+crystalBallFit->GetParameter(14); // mean gauss3 + // + vsigma[0] = crystalBallFit->GetParameter(2); // sigma gauss1 left + vsigma[1] = crystalBallFit->GetParameter(9); // sigma gauss1 right + vsigma[2] = crystalBallFit->GetParameter(12); // sigma gauss2 + vsigma[3] = crystalBallFit->GetParameter(15); // sigma gauss3 + // + vnorm[0] = crystalBallFit->GetParameter(0); // norm gauss1 + vnorm[1] = crystalBallFit->GetParameter(10); // norm gauss2 + vnorm[2] = crystalBallFit->GetParameter(13); // norm gauss3 + // + vexp[0] = crystalBallFit->GetParameter(3); // norm exp + vexp[1] = crystalBallFit->GetParameter(4); // slope exp + // + for(int i=0; i<2; i++) vpolLeft[i] = crystalBallFit->GetParameter(i+5); + // + for(int i=0; i<2; i++) vpolRight[i] = crystalBallFit->GetParameter(i+16); + // + for(int i=0; i<6; i++) voption[i] = crystalBallFit->GetParameter(i+19); + + cout<<"\nfit integral [0, 2] : "<Integral(0., 2.) <Integral(0., vrange[0]) <Integral(vrange[0], vrange[1])<Integral(vrange[1], vmean[0]) <Integral(vmean[0], vrange[2]) <Integral(vrange[2], 2.) < drawOpt{sum, polL, exp, g1, g2, g3, polR}; + + unsigned long nfnc{drawOpt.size()}; + + TF1 *fnc[nfnc]; + + for(unsigned i=0; iGetNpar()); + + for( int par=0; parGetNpar(); par++ ){ + fnc[i]->SetParName (par, crystalBallFit->GetParName(par)); + fnc[i]->SetParameter(par, crystalBallFit->GetParameter(par)); + } + for(unsigned j=0; jSetParameter(j+19, drawOpt[i][j]); + + if(i==0) fnc[i]->SetLineColor(kBlack); // sum + if(i==1) fnc[i]->SetLineColor(kGray+2); // polL + if(i==2) fnc[i]->SetLineColor(kCyan); // exp + if(i==3){fnc[i]->SetLineColor(kRed); vintegral[0]=fnc[i]->Integral(0,2)/hrsp->GetBinWidth(1);} // g1 + if(i==4){fnc[i]->SetLineColor(kMagenta); vintegral[1]=fnc[i]->Integral(0,2)/hrsp->GetBinWidth(1);} // g2 + if(i==5){fnc[i]->SetLineColor(kOrange+4); vintegral[2]=fnc[i]->Integral(0,2)/hrsp->GetBinWidth(1);} // g3 + if(i==6) fnc[i]->SetLineColor(kGray+2); // polR + if(drawCrystalBall && canvas!=0){ canvas->cd(2); fnc[i]->Draw("same"); } + } + + + // + // print fit info + // + cout<<"\norder after fit: "<GetXaxis()->GetBinCenter(hrsp->GetMaximumBin())}; + cout<<"|mu1-peak|: "<