Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:14:11 

0001 #include "sPhenixStyle.h"
0002 #include "sPhenixStyle.C"
0003 #include "respFitterDef.h"
0004 const int nIter = 10;
0005 
0006 string generateFit(int order, string function);
0007 
0008 void numericInverse_2D(int save = 1)
0009 {
0010   SetsPhenixStyle();
0011   
0012   vector<float> etaBins;
0013   for(float i = etaStart; i <= etaEnd; i+=etaGap)
0014     {
0015       etaBins.push_back(i);
0016     }
0017   const int nEtaBins = etaBins.size() -1 + extraBins;
0018   std::cout << "nEtaBins: " << nEtaBins << std::endl;
0019   TFile *fin = new TFile("calibHists_hists_R04_Corr0_isLin1_2D.root");
0020 
0021   float fitStart =8;
0022   float fitEnd = 60;
0023   
0024   //draw and fit our jet energy linearity
0025   int colors[] = {1, 2, 4, kGreen +2, kViolet, kCyan, kOrange, kMagenta+2, kAzure-2};
0026   TF1 *unity = new TF1("unity","x",0,80);
0027   unity -> SetLineColor(1);
0028   unity -> SetLineStyle(8);
0029   
0030   TLegend *leg = new TLegend(0.6,0.2,0.9,0.5);
0031   leg -> SetFillStyle(0);
0032   leg -> SetBorderSize(0);
0033 
0034   TCanvas *c1 = new TCanvas("linearity","linearity");
0035   TCanvas *cUnity= new TCanvas("UnityDiff","UnityDiff");
0036   TCanvas *cDiff = new TCanvas("fitDiff","FitDiff");
0037   TGraphErrors *unityDiff[nEtaBins];
0038   TGraphErrors *chi2NDF[nEtaBins];
0039   TGraphErrors *fitDiff[nEtaBins][3];
0040   TF1 *linFit[nEtaBins];
0041   TH1F *pt_true_matched[nEtaBins];
0042   TGraphErrors *lin[nEtaBins];
0043              
0044 
0045   for(int i = 0; i < nEtaBins-1; i++)
0046     {
0047       std::cout << "etaBin: " << i << std::endl;
0048       linFit[i] = new TF1(Form("linFit_eta%d",i),"[0] + [1]*x",fitStart,fitEnd);
0049       linFit[i] -> SetLineColor(colors[i]);
0050       pt_true_matched[i] = (TH1F*)fin -> Get(Form("h_pt_true_matched_eta%d",i));
0051       lin[i] = (TGraphErrors*)fin -> Get(Form("g_jes_eta%d",i));
0052       lin[i] -> SetMarkerColor(colors[i]);
0053       lin[i] -> GetYaxis() -> SetRangeUser(0,70);
0054       lin[i] -> GetXaxis() -> SetLimits(0,80);
0055       lin[i] -> SetTitle(";p_{T,truth} [GeV/c];#LTp_{T,reco}#GT [GeV/c]");
0056       c1 -> cd();
0057       if(i == 0)lin[i] -> Draw("ap");
0058       else lin[i] -> Draw("samep");
0059       unityDiff[i] = new TGraphErrors();
0060       for(int j = 0; j < lin[i] -> GetN(); j++)
0061     {
0062       Double_t x, y;
0063       lin[i] -> GetPoint(j,x,y);
0064       
0065       unityDiff[i] -> SetPoint(j,x,y-x);
0066     }
0067       lin[i] -> Fit(linFit[i],"RQ");
0068       fitDiff[i][0] = new TGraphErrors();
0069       for(int n = 0; n < lin[i] -> GetN(); n++)
0070     {
0071       Double_t x,y; 
0072       lin[i] -> GetPoint(n,x,y);
0073               
0074       if(x < fitEnd)
0075         { 
0076           fitDiff[i][0] -> SetPoint(n, x, (y - linFit[i] -> Eval(x))/y);
0077         }
0078     }
0079 
0080       fitDiff[i][0] -> SetTitle(";p_{T,truth} [GeV/c];#frac{(Data - Fit)}{Data}");
0081       fitDiff[i][0] -> SetMarkerColor(colors[i]);
0082       fitDiff[i][0] -> GetYaxis() -> SetRangeUser(-1,1);
0083 
0084       cUnity -> cd();
0085 
0086       unityDiff[i] -> SetTitle(";p_{T,truth} [GeV/c];#LTp_{T,reco}#GT - p_{T,Truth} [GeV/c]");
0087       unityDiff[i] -> GetYaxis() -> SetRangeUser(-30,5);
0088       unityDiff[i] -> SetMarkerColor(colors[i]);
0089 
0090       if(i == 0)unityDiff[i] -> Draw("ap");
0091       else unityDiff[i] -> Draw("samep");
0092 
0093       if(i <=nEtaBins - 3)leg -> AddEntry(lin[i],Form("%g<#eta<%g",etaBins.at(i),etaBins.at(i+1)),"p");
0094       else leg -> AddEntry(lin[i],"-0.6 < #eta < 0.6","p");
0095 
0096       cDiff -> cd();
0097       if(i == 0) fitDiff[i][0] -> Draw("ap");
0098       else fitDiff[i][0] -> Draw("samep");
0099     }
0100   c1 -> cd(); leg -> Draw("same"); unity -> Draw("same");
0101   cUnity -> cd(); leg -> Draw("same"); 
0102   cDiff -> cd(); leg -> Draw("same");
0103   if(save)
0104     {
0105       c1 -> SaveAs("plots_correctionFromLin/lin_isoStack.pdf");
0106       cUnity -> SaveAs("plots_correctionFromLin/unityDiff.pdf");
0107       cDiff -> SaveAs("plots_correctionFromLin/fitDiff0.pdf");
0108     }
0109   
0110   //for the numerical inversion, first we must compute
0111   //the function R, which is the avaerage JES
0112   
0113   TGraphErrors *r_lin[nEtaBins];
0114   
0115   const int nPoints = lin[0] -> GetN();
0116   TF1 *rFit[nEtaBins][nIter]; 
0117   
0118   TCanvas *cChi2 = new TCanvas("chi2","chi2/NDF");
0119   TCanvas *c2 = new TCanvas("response","Response");
0120   int bestIter[nEtaBins];
0121   TLine *bestIterLn[nEtaBins];
0122   
0123   for(int i = 0; i < nEtaBins-1; i++)
0124     {
0125       r_lin[i] = new TGraphErrors();
0126       
0127       c2 -> cd();
0128      
0129       for(int j = 0; j < nPoints; j++)
0130     {
0131       Double_t x, y, erry, errx;
0132       lin[i] -> GetPoint(j, x, y);
0133       erry = lin[i] -> GetErrorY(j);
0134       errx = lin[i] -> GetErrorX(j);
0135       pt_true_matched[i]-> GetXaxis() -> SetRangeUser(x - errx, x + errx);
0136       x = pt_true_matched[i] -> GetMean();
0137       
0138       r_lin[i] -> SetPoint(j, x, y/x);
0139       float err = y/x * sqrt(pow(erry/y,2) /*+ pow((errx/2)/x,2)*/);
0140       r_lin[i] -> SetPointError(j, errx, err);
0141       
0142     }
0143       r_lin[i] -> SetMarkerColor(colors[i]);
0144       r_lin[i] -> SetTitle(";p_{T,truth} [GeV/c];#LTp_{T,reco}#GT/#LTp_{T,truth}#GT");
0145       
0146       string fit = "[0]";
0147       chi2NDF[i] = new TGraphErrors();
0148       chi2NDF[i] -> SetTitle(";Iterations;chi2/NDF");
0149       chi2NDF[i] -> SetMarkerColor(colors[i]);
0150       float minChi2NDF = 9999;
0151       
0152 
0153       //Double_t xstart, ystart;
0154       //r_lin[i] -> GetPoint(1,xstart,ystart);
0155       for(int l = 0; l < nIter; l++)
0156     {
0157       rFit[i][l] = new TF1("polyFitBase",fit.c_str(),fitStart,fitEnd);
0158       rFit[i][l] -> SetLineColor(colors[i]);
0159 
0160     
0161 
0162       r_lin[i] -> Fit(rFit[i][l],"Q0","0",fitStart,fitEnd);
0163      
0164       r_lin[i] -> GetYaxis() -> SetRangeUser(0.,1);
0165       float chi2, ndf;
0166       chi2 = rFit[i][l] -> GetChisquare();
0167       ndf = rFit[i][l] -> GetNDF();
0168       chi2NDF[i] -> SetPoint(l,l,chi2/ndf);
0169       
0170 
0171       
0172       if(chi2/ndf < minChi2NDF)
0173         {
0174           minChi2NDF = chi2/ndf; 
0175           bestIter[i] = l;
0176         }
0177       
0178           
0179       fit = generateFit(l+1,fit);
0180     }
0181       
0182       fitDiff[i][1] = new TGraphErrors();
0183       for(int n = 0; n < r_lin[i] -> GetN(); n++)
0184     {
0185       Double_t x,y;
0186       r_lin[i] -> GetPoint(n,x,y);
0187       if(x < fitEnd)
0188         {
0189           fitDiff[i][1] -> SetPoint(n, x, (y - rFit[i][bestIter[i]] -> Eval(x))/y);
0190         }
0191     }
0192       fitDiff[i][1] -> SetMarkerColor(colors[i]);
0193       fitDiff[i][1] -> SetTitle(";p_{T,truth} [GeV/c];#frac{(Data - Fit)}{Data}");
0194       fitDiff[i][1] -> GetYaxis() -> SetRangeUser(-1,1);
0195 
0196       bestIterLn[i] = new TLine(bestIter[i],0,bestIter[i],20);
0197       bestIterLn[i] -> SetLineColor(colors[i]);
0198       bestIterLn[i] -> SetLineStyle(8);
0199       
0200       
0201       c2 -> cd();
0202       if(i == 0)
0203     {
0204       r_lin[i] -> GetXaxis() -> SetLimits(0,80);
0205       r_lin[i] -> Draw("ap");
0206     }
0207       else r_lin[i] -> Draw("samep");
0208       rFit[i][bestIter[i]] -> Draw("same");
0209       
0210       
0211       cChi2 -> cd(); 
0212       if(i == 0)
0213     {
0214       chi2NDF[i] -> GetYaxis() -> SetRangeUser(0,20);
0215       chi2NDF[i] -> Draw("ap");
0216         
0217     }
0218       else chi2NDF[i] -> Draw("samep");
0219       bestIterLn[i] -> Draw("same");
0220 
0221       cDiff -> cd();
0222       if(i == 0)fitDiff[i][1] -> Draw("ap");
0223       else fitDiff[i][1] -> Draw("samep");
0224     }
0225   c2 -> cd();
0226   leg -> Draw("same");
0227   //gPad -> SetLogy();
0228   leg -> Draw("same");
0229   cDiff -> cd(); 
0230   leg -> Draw("same");
0231   
0232   if(save)
0233     {
0234       c2 -> SaveAs("plots_correctionFromLin/fMeanScaled.pdf");
0235       cChi2 -> SaveAs("plots_correctionFromLin/chi2NdffMean.pdf");
0236       cDiff -> SaveAs("plots_correctionFromLin/fitDiff1.pdf");
0237     }
0238   //Now, in order to get the calibration factor, we have to compute R^-1
0239   //Which is defined by R^-1(y) = R(f^-1(y))
0240   
0241   TGraphAsymmErrors *rTilde[nEtaBins];
0242   TGraphErrors *gfInv[nEtaBins];
0243   TF1 *correctionFit[nEtaBins][nIter];
0244   TF1 *fitLinearExtrap[nEtaBins];
0245   TCanvas *c3 = new TCanvas();
0246   TCanvas *c4 = new TCanvas();
0247   TCanvas *cChi22 = new TCanvas();
0248 
0249   for(int i = 0; i < nEtaBins-1; i++)
0250     {
0251       rTilde[i] = new TGraphAsymmErrors();
0252       gfInv[i] = new TGraphErrors();
0253      
0254       float lastPoint = 0;
0255       for(int j = 0; j < nPoints; j++)
0256     {
0257       Double_t x, y;
0258       lin[i] -> GetPoint(j, x, y);
0259       if(x < fitStart || x > fitEnd) continue;
0260       float fInv = linFit[i] -> GetX(y);
0261       gfInv[i] -> SetPoint(j, y, fInv);
0262       float R = rFit[i][bestIter[i]] -> Eval(fInv);
0263       rTilde[i] -> SetPoint(j, y, 1./R);
0264       lastPoint = y;      
0265     }
0266       rTilde[i] -> SetMarkerColor(colors[i]);
0267       gfInv[i] -> SetMarkerColor(colors[i]);
0268 
0269       rTilde[i] -> SetTitle(";#LTp_{T,reco}#GT [GeV/c];1/#tilde{R}(p_{T,reco})");
0270       gfInv[i] -> SetTitle(";#LTp_{T,reco}#GT [GeV/c];p_{T,truth} [GeV/c]");
0271 
0272       rTilde[i] -> GetYaxis() -> SetRangeUser(0,2);
0273       string fit = "[0]";
0274       chi2NDF[i] = new TGraphErrors();
0275       chi2NDF[i] -> SetTitle(";Iterations;chi2/NDF");
0276       chi2NDF[i] -> SetMarkerColor(colors[i]);
0277       float minChi2NDF = 9999;
0278       for(int l = 0; l < nIter; l++)
0279     {
0280       correctionFit[i][l] = new TF1("polyFitBas2",fit.c_str(),fitStart,52);
0281       
0282       rTilde[i] -> Fit(correctionFit[i][l],"Q0","",10,70);
0283       float chi2, ndf;
0284       chi2 = correctionFit[i][l] -> GetChisquare();
0285       ndf = correctionFit[i][l] -> GetNDF();
0286       chi2NDF[i] -> SetPoint(l,l,chi2/ndf);
0287       if(chi2/ndf < minChi2NDF)
0288         {
0289           minChi2NDF = chi2/ndf;
0290           bestIter[i] = l;
0291         }
0292       
0293       fit = generateFit(l+1,fit);
0294     }
0295     
0296       correctionFit[i][bestIter[i]] -> SetLineColor(colors[i]);
0297       
0298       fitDiff[i][2] = new TGraphErrors();
0299       for(int n = 0; n < rTilde[i] -> GetN(); n++)
0300     {
0301       Double_t x,y;
0302       rTilde[i] -> GetPoint(n,x,y);
0303       if(x < fitEnd)
0304         {
0305           fitDiff[i][2] -> SetPoint(n, x, (y - correctionFit[i][bestIter[i]] -> Eval(x))/y);
0306         }
0307     }
0308       
0309       fitDiff[i][2] -> SetTitle(";p_{T,reco} [GeV/c];#frac{(Data - Fit)}{Data}");
0310       fitDiff[i][2] -> SetMarkerColor(colors[i]);
0311       fitDiff[i][2] -> GetYaxis() -> SetRangeUser(-1,1);
0312 
0313       c3 -> cd();
0314       if(i == 0)rTilde[i] -> Draw("ap");
0315       else rTilde[i] -> Draw("samep");
0316       correctionFit[i][bestIter[i]] -> SetRange(fitStart,lastPoint);
0317       correctionFit[i][bestIter[i]] -> Draw("same");
0318       fitLinearExtrap[i] = new TF1(Form("corrFit_eta%d_extrap",i),"(x>[0])*[1]");
0319       std::cout << "Setting last point to: " << lastPoint << std::endl;
0320       fitLinearExtrap[i] -> FixParameter(0,lastPoint);
0321       std::cout << "Last point eval: " << correctionFit[i][bestIter[i]]->Eval(lastPoint) << std::endl;
0322       fitLinearExtrap[i] -> FixParameter(1,correctionFit[i][bestIter[i]]->Eval(lastPoint));
0323       fitLinearExtrap[i] -> SetLineColor(colors[i]);
0324       fitLinearExtrap[i] -> SetLineStyle(8);
0325       fitLinearExtrap[i] -> SetRange(lastPoint,80);
0326       fitLinearExtrap[i] -> Draw("same");
0327 
0328       c4 -> cd();
0329       if(i == 0)gfInv[i] -> Draw("ap");
0330       else gfInv[i] -> Draw("samep");
0331 
0332       bestIterLn[i] = new TLine(bestIter[i],-0.005,bestIter[i],0.03);
0333       bestIterLn[i] -> SetLineColor(colors[i]);
0334       bestIterLn[i] -> SetLineStyle(8);
0335       
0336       cChi22 -> cd(); 
0337       if(i == 0)
0338     {
0339       chi2NDF[i] -> GetYaxis() -> SetRangeUser(-0.005,0.03);
0340       chi2NDF[i] -> GetXaxis() -> SetLimits(-1,14.5);
0341 
0342       chi2NDF[i] -> Draw("ap");
0343     }
0344       else chi2NDF[i] -> Draw("samep");
0345       bestIterLn[i] -> Draw("same");
0346 
0347       cDiff -> cd();
0348       if(i == 0)fitDiff[i][2] -> Draw("ap");
0349       else fitDiff[i][2] -> Draw("samep");
0350       leg -> Draw("same");
0351     }
0352   c3 -> cd();
0353   leg -> Draw("same");
0354   c4 -> cd();
0355   leg -> Draw("same");
0356   cDiff -> cd();
0357   leg -> Draw("same");
0358   if(save)
0359     {
0360       c3 -> SaveAs("plots_correctionFromLin/jesCorrection.pdf");
0361       c4 -> SaveAs("plots_correctionFromLin/jesInv.pdf");
0362       cChi22 -> SaveAs("plots_correctionFromLin/chi2NdfCorrFit.pdf");
0363       cDiff -> SaveAs("plots_correctionFromLin/fitDiff2.pdf");
0364       
0365       TFile *corrOut = new TFile("JES_IsoCorr_NumInv.root","RECREATE");
0366       corrOut -> cd();
0367       for(int i = 0; i < nEtaBins-1; i++)
0368     {
0369       correctionFit[i][bestIter[i]] -> SetName(Form("corrFit_eta%d",i));
0370       correctionFit[i][bestIter[i]] -> Write();
0371       fitLinearExtrap[i] -> Write();
0372     }
0373       corrOut -> Close();
0374     }
0375 }
0376   
0377 string generateFit(int order, string function)
0378 {
0379   string nextTerm = Form("+ [%d]*pow(log(x),%d)", order, order);
0380   
0381   function += nextTerm;
0382   
0383   //TF1 *fit = new TF1("polyLog",function.c_str(),11,60);
0384   
0385   return function;
0386 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!