diff --git a/JetUtilities/src/L2Creator.cc b/JetUtilities/src/L2Creator.cc index de67892a..83fc9a86 100644 --- a/JetUtilities/src/L2Creator.cc +++ b/JetUtilities/src/L2Creator.cc @@ -24,6 +24,7 @@ L2Creator::L2Creator() { histogramMetric = HistUtil::getHistogramMetricType(histMet); delphes = false; maxFitIter = 30; + ptclipfit=false; } //______________________________________________________________________________ @@ -46,6 +47,10 @@ L2Creator::L2Creator(CommandLine& cl) { histMet = cl.getValue ("histMet", "mu_h"); histogramMetric = HistUtil::getHistogramMetricType(histMet); + ptclip = cl.getValue ("ptclip", 0.); + statTh = cl.getValue ("statTh", 4); + ptclipfit = cl.getValue ("ptclipfit", false); + if (!cl.partialCheck()) return; cl.print(); } @@ -211,6 +216,14 @@ void L2Creator::loopOverEtaBins() { TH1F* hrsp(0); hl_rsp.begin_loop(); + //print fit results for all eta bins in a txt file + TString txtFitResultsFilename = outputDir + "Fit_Results.txt"; + ofstream FitResults(txtFitResultsFilename); + FitResults.setf(ios::right); + + //edw gia draw prob,chi2 vs eta + double etabin[82], Chi2NDF[82], Prob[82], Prob_0p009[82]; + while ((hrsp=hl_rsp.next_object(indices))) { unsigned int ieta=indices[0]; @@ -232,7 +245,7 @@ void L2Creator::loopOverEtaBins() { // only add points to the graphs if the current histo is not empty // the current setting might be a little high // - if (hrsp->GetEntries() > 4) {//hrsp->Integral()!=0) { + if (hrsp->GetEntries() > 30) {//hrsp->Integral()!=0) {//EDW 4 //TF1* frsp = (TF1*)hrsp->GetListOfFunctions()->Last(); //std::cout << "hrspName = " << hrsp->GetName() << ": frsp = " << frsp << std::endl; @@ -277,19 +290,23 @@ void L2Creator::loopOverEtaBins() { double eabsrsp = epeak; double abscor = 0.0; double eabscor = 0.0; - - if (absrsp > 0) { +//cout << "EDW 1 " << "eta" << vabscor_eta.back()->GetName() << "absrsp " << absrsp << " and eabsrsp " << eabsrsp << " and abscor " << abscor << " and eabscor " << eabscor << endl; + if (absrsp > 0 ) {//edw test abscor =1.0/absrsp; eabscor = abscor*abscor*epeak; + + +//cout << "EDW 2 " << "eta" << vabscor_eta.back()->GetName() << endl; } - if ((abscor>0) && (absrsp>0) && (eabscor>1e-5) && (eabscor/abscor<0.5) && (eabsrsp>1e-4) && (eabsrsp/absrsp<0.5)) { + 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); + // cout << "EDW 3" << "eta" << vabscor_eta.back()->GetName() << "refpt" << refpt << "jetpt" << jetpt << "absrsp " << absrsp << " and eabsrsp " << eabsrsp << " and abscor " << abscor << " and eabscor " << eabscor << endl; } - else cout << "absrsp " << absrsp << " and eabsrsp " << eabsrsp << " and abscor " << abscor << " and eabscor " << eabscor << endl; + //else cout << "EDW 4" << "absrsp " << absrsp << " and eabsrsp " << eabsrsp << " and abscor " << abscor << " and eabscor " << eabscor << endl; } // @@ -376,6 +393,25 @@ void L2Creator::loopOverEtaBins() { fabscor=new TF1("fit",fcn.Data(),xmin,xmax); setOfflinePFParameters(gabscor, fabscor,xmin,xmax); + std::cout<<"ieta = "<77){ //for PUPPI + fabscor->FixParameter(3,0.); + fabscor->FixParameter(6,0.); + } + } + if(l2pffit.Contains("spline",TString::kIgnoreCase)) { vabscor_eta_spline.back()->setPartialFunction(fabscor); } @@ -537,6 +573,54 @@ void L2Creator::loopOverEtaBins() { } } + + //edw ptclipfit + if (ptclipfit) + { + if (xmin > 0.0001) + { + int nPar = fabscor->GetNpar(); + int clipPar = nPar; + + TString clip = TString::Format("((x<10)*([%d]))+((x>=10)*("+(TString)fabscor->GetTitle()+"))", clipPar); + TF1 * fabscornew = new TF1(fabscor->GetName(), clip, 0.001, 6500); + for (int ip=0; ipSetParameter(ip, fabscor->GetParameter(ip)); + } + + fabscornew->SetParameter(clipPar, fabscor->Eval(xmin)); + + fabscornew->SetChisquare(fabscor->GetChisquare()); + fabscornew->SetNDF(fabscor->GetNDF()); + + fabscor = fabscornew; + gabscor->GetListOfFunctions()->Clear(); + gabscor->GetListOfFunctions()->AddLast(fabscor); + } + } + + //EDW print chi2 and prob for each fit in every eta bin + std::cout<<"Chi2/NDF = "<GetChisquare()/fabscor->GetNDF()<GetName()<<"\n"; + FitResults<<"Chi2/NDF = "<GetChisquare()<<"/"<GetNDF()<<" = "<GetChisquare()/fabscor->GetNDF()<<"\n"; + FitResults<<"Prob = "<GetProb()<<"\n\n\n"; + FitResults<<"Number of points = "<GetChisquare()/fabscor->GetNDF()>1000) {Chi2NDF[ieta] =0; cout<GetChisquare()/fabscor->GetNDF(); + + Prob[ieta] = fabscor->GetProb(); + if(fabscor->GetProb()<0.009) Prob_0p009[ieta] = fabscor->GetProb(); + else Prob_0p009[ieta] = 0.01; + + // // format the graphs // @@ -553,6 +637,36 @@ void L2Creator::loopOverEtaBins() { vabscor_eta_spline.push_back(nullptr); } } + + FitResults.close(); + + + //===== Draw graph Chi2/NDF vs Eta ===== + TGraph *g_chi2ndf_vs_eta = new TGraph(82, etabin, Chi2NDF); + TCanvas *pad1 = new TCanvas("pad1", "",1); + g_chi2ndf_vs_eta->GetYaxis()->SetTitle("Chi2/ndf"); + g_chi2ndf_vs_eta->GetXaxis()->SetTitle("eta"); + g_chi2ndf_vs_eta->SetMarkerStyle(20); + g_chi2ndf_vs_eta->Draw("AP"); + pad1->SaveAs(Form("%s/Chi2NDF_vs_eta.png",outputDir.Data())); + + TGraph *g_prob_vs_eta = new TGraph(82, etabin, Prob); + TCanvas *pad2 = new TCanvas("pad2", "",1); + g_prob_vs_eta->GetYaxis()->SetTitle("Prob."); + g_prob_vs_eta->GetXaxis()->SetTitle("eta"); + g_prob_vs_eta->SetMarkerStyle(20); + g_prob_vs_eta->Draw("AP"); + pad2->SaveAs(Form("%s/Prob_vs_eta.png",outputDir.Data())); + + TGraph *g_prob_0p009_vs_eta = new TGraph(82, etabin, Prob_0p009); + TCanvas *pad3 = new TCanvas("pad3", "",1); + g_prob_0p009_vs_eta->GetYaxis()->SetTitle("Prob"); + g_prob_0p009_vs_eta->GetXaxis()->SetTitle("eta"); + g_prob_0p009_vs_eta->SetMarkerStyle(20); + g_prob_0p009_vs_eta->Draw("AP"); + pad3->SaveAs(Form("%s/prob_0p009_vs_eta.png",outputDir.Data())); + + } //______________________________________________________________________________ @@ -804,7 +918,7 @@ void L2Creator::makeCanvas(string makeCanvasVariable) { spline_func->SetRange(bounds.first,bounds.second); spline_func->SetLineColor(igraph%nperpad+1); spline_func->SetLineStyle(kDotted); - spline_func->Draw("same"); +// spline_func->Draw("same"); } } @@ -874,7 +988,8 @@ TString L2Creator::getOfflinePFFunction() { 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])))))"; + 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])))))"; +// return "[0]+[1]*pow(x,[2])+[3]*exp(-[4]*max(0.,(log10(x)-[5]))*(log10(x)-[5]))"; //Modification of 0,1,2 with powerlaw first term + single-sided log-gaus by Mikko } 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])))"; @@ -939,7 +1054,7 @@ void L2Creator::setOfflinePFParameters(TGraphErrors* gabscor, TF1* fabscor, doub fabscor->SetParLimits(4,0,100); } else if(l2pffit.EqualTo("standard+Gaussian",TString::kIgnoreCase)) { - fabscor->SetParameter(0,-0.0221278); +/* fabscor->SetParameter(0,-0.0221278); fabscor->SetParameter(1,119.265); fabscor->SetParameter(2,100); fabscor->SetParameter(3,-0.0679365); @@ -954,6 +1069,34 @@ void L2Creator::setOfflinePFParameters(TGraphErrors* gabscor, TF1* fabscor, doub fabscor->SetParLimits(4,0,500); fabscor->SetParLimits(0,-2,25); fabscor->SetParLimits(1,0,250); +*/ + fabscor->SetParameter(0,0.0221278); //v2 + fabscor->SetParameter(1,14.265); + fabscor->SetParameter(2,10); + 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,5); + fabscor->SetParLimits(5,0,50); + fabscor->SetParLimits(0,-2,50); + fabscor->SetParLimits(2,0,200); + fabscor->SetParLimits(1,0,250); + + +/* fabscor->SetRange(8.,3500.); + fabscor->SetParameter(0,1.09892); //Modification of 0,1,2 with powerlaw first term + single-sided log-gaus by Mikko + fabscor->SetParameter(1,-0.46266); + fabscor->SetParameter(2,-0.152203); + fabscor->SetParameter(3,0.0951505); + fabscor->SetParameter(4,1.16479); + fabscor->SetParameter(5,1.53048); +*/ } else if(l2pffit.EqualTo("LogNormal+Gaussian+fixed",TString::kIgnoreCase)) { fabscor->SetRange(3,2000); @@ -1313,6 +1456,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //For eta-dependent spline clipping int pt_limit = 70; + float pt_clip = ptclip; unsigned int vector_size = 0; vector_size = vabscor_eta.size(); @@ -1345,6 +1489,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { bool abovePtLimit = false; bool lastLine = false; + bool firstline = true; for(int isection=0; isectiongetNSections(); isection++) { if(lastLine) continue; @@ -1364,7 +1509,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //if(isection==spline->getNSections()-1) { // bounds.second = 6500; //} - + if(bounds.second < pt_clip) continue; if(isection==spline->getNSections()-1) lastLine = true; if(bounds.second >= pt_limit) { abovePtLimit = true; @@ -1373,16 +1518,17 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() { //For expediency of Summer16_25nsV5_MC do eta-dependent clipping fout<getNpar()+2) //Number of parameters + 2 - <setParameters(isection); for(int p=0; pgetNpar(); p++) { fout<GetParameter(p); } fout<