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