Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:13:22

0001 // ----------------------------------------------------------------------------
0002 // 'MakeRatioComparisonPlot.C'
0003 // Derek Anderson
0004 // 01.03.2022
0005 //
0006 // Use this quickly plot a set of numerator
0007 // distributions against a set of denominator
0008 // distributions
0009 // ----------------------------------------------------------------------------
0010 
0011 #include <iostream>
0012 #include "TH1.h"
0013 #include "TPad.h"
0014 #include "TFile.h"
0015 #include "TLine.h"
0016 #include "TError.h"
0017 #include "TString.h"
0018 #include "TLegend.h"
0019 #include "TCanvas.h"
0020 #include "TPaveText.h"
0021 
0022 using namespace std;
0023 
0024 // global constants
0025 static const UInt_t NHist(2);
0026 static const UInt_t NPlot(2);
0027 static const UInt_t NPad(2);
0028 static const UInt_t NVtx(4);
0029 static const UInt_t NTxt(3);
0030 
0031 
0032 
0033 void MakeRatioComparisonPlot() {
0034 
0035   // lower verbosity
0036   gErrorIgnoreLevel = kError;
0037   cout << "\n  Beginning plot macro..." << endl;
0038 
0039   // file parameters
0040   const TString sOutput("debugCorrelators.alexVsF4A_truthVsRecoAfterPhiFix_ptJet30.pp200run6jet40.d18m4y2023.root");
0041   const TString sInDenom[NHist] = {"reference/alex_output/pp200run6jet40.SphenixTruthUnscaled_output_trksAndChrgPars.root",
0042                                    "reference/alex_code/SphenixRecoUnscaled_output_derekTestRun.root"};
0043   const TString sInNumer[NHist] = {"./output/pp200run6jet40.forDebug_true_withPhiFix.d18m4y2023.root",
0044                                    "./output/pp200run6jet40.forDebug_reco_withPhiFix.d18m4y2023.root"};
0045 
0046   // denominator parameters
0047   const TString sHeadDenom("#bf{Reference Output}");
0048   const TString sHistDenom[NHist]  = {"EEC2",           "EEC2"};
0049   const TString sNameDenom[NHist]  = {"hReferencTruth", "hReferenceReco"};
0050   const TString sLabelDenom[NHist] = {"Truth",
0051                                       "Reco."};
0052 
0053   // numerator parameters
0054   const TString sHeadNumer("#bf{F4A Implementation}");
0055   const TString sHistNumer[NHist]  = {"hCorrelatorDrAxis_ptJet30", "hCorrelatorDrAxis_ptJet30"};
0056   const TString sNameNumer[NHist]  = {"hFun4AllTruth",             "hFun4AllReco"};
0057   const TString sLabelNumer[NHist] = {"Truth",
0058                                       "Reco."};
0059 
0060   // plot parameters
0061   const TString sTitle("");
0062   const TString sTitleX("#DeltaR");
0063   const TString sTitleY("EEC");
0064   const TString sTitleR("[f4a] / [reference]");
0065   const TString sNameRatio[NHist] = {"hRatioTruth", "hRatioReco"};
0066   const TString sOptDenom[NHist]  = {"",           "sames"};
0067   const TString sOptNumer[NHist]  = {"sames",      "sames"};
0068   const TString sOptRatio[NHist]  = {"",           "sames"};
0069   const TString sText[NTxt]       = {"#bf{#it{sPHENIX}} Simulation [Run6]", "JS 40 GeV Jet Sample", "p_{T}^{jet} #in (30, 50) GeV"};
0070   const Float_t xPlotRange[NPlot] = {0.001, 1.};
0071   const UInt_t  fColDen[NHist]    = {923, 634};
0072   const UInt_t  fColNum[NHist]    = {921, 895};
0073   const UInt_t  fMarDen[NHist]    = {20, 29};
0074   const UInt_t  fMarNum[NHist]    = {24, 30};
0075 
0076   // norm/rebin parameters
0077   const Bool_t doIntNorm(false);
0078   const Bool_t doRebinDenom[NHist] = {false, false};
0079   const Bool_t doRebinNumer[NHist] = {false, false};
0080   const UInt_t nRebinDenom[NHist]  = {2, 2};
0081   const UInt_t nRebinNumer[NHist]  = {2, 2};
0082 
0083   // open files
0084   TFile *fOutput = new TFile(sOutput.Data(), "recreate");
0085   if (!fOutput) {
0086     cerr << "PANIC: couldn't open output file!\n" << endl;
0087     return;
0088   }
0089 
0090   TFile *fDenom[NHist];
0091   TFile *fNumer[NHist];
0092   for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0093     fDenom[iHist] = new TFile(sInDenom[iHist].Data(), "read");
0094     fNumer[iHist] = new TFile(sInNumer[iHist].Data(), "read");
0095     if (!fDenom[iHist] || !fNumer[iHist]) {
0096       cerr << "PANIC: couldn't open denominator or numerator file #" << iHist << "!\n"
0097            << "       fDenom = " << fDenom[iHist] << ", fNumer = " << fNumer[iHist] << "\n"
0098            << endl;
0099       return;
0100     }
0101   }
0102   cout << "    Opened files." << endl;
0103 
0104   // grab histograms
0105   TH1D *hDenom[NHist];
0106   TH1D *hNumer[NHist];
0107   for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0108     hDenom[iHist] = (TH1D*) fDenom[iHist] -> Get(sHistDenom[iHist].Data());
0109     hNumer[iHist] = (TH1D*) fNumer[iHist] -> Get(sHistNumer[iHist].Data());
0110     if (!hDenom[iHist] || !hNumer[iHist]) {
0111       cerr << "PANIC: couldn't grab numerator or denominator histogram # " << iHist << "!\n"
0112            << "       hDenom = " << hDenom[iHist] << ", hNumer = " << hNumer[iHist] << "\n"
0113            << endl;
0114       return;
0115     }
0116     hDenom[iHist] -> SetName(sNameDenom[iHist].Data());
0117     hNumer[iHist] -> SetName(sNameNumer[iHist].Data());
0118   }
0119   cout << "    Grabbed histograms." << endl;
0120 
0121   // rebin histograms
0122   Bool_t doRebin(false);
0123   for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0124     doRebin = (doRebinDenom[iHist] || doRebinNumer[iHist]);
0125     if (doRebin) break;
0126   }
0127 
0128   if (doRebin) {
0129     for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0130       if (doRebinDenom[iHist]) hDenom[iHist] -> Rebin(nRebinDenom[iHist]);
0131       if (doRebinNumer[iHist]) hNumer[iHist] -> Rebin(nRebinNumer[iHist]);
0132     }
0133     cout << "    Rebinned histograms." << endl;
0134   }
0135 
0136   // normalize by integrals )if needed)
0137   if (doIntNorm) {
0138     for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0139       const Double_t intDenom = hDenom[iHist] -> Integral();
0140       const Double_t intNumer = hNumer[iHist] -> Integral();
0141       hDenom[iHist] -> Scale(1. / intDenom);
0142       hNumer[iHist] -> Scale(1. / intNumer);
0143     }
0144     cout << "    Normalized histograms by integral." << endl;
0145   }
0146 
0147   // calculate ratios
0148   TH1D *hRatio[NHist];
0149   for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0150     hRatio[iHist] = (TH1D*) hDenom[iHist] -> Clone();
0151     hRatio[iHist] -> Reset("ICE");
0152     hRatio[iHist] -> Divide(hNumer[iHist], hDenom[iHist], 1., 1.);
0153     hRatio[iHist] -> SetName(sNameRatio[iHist]);
0154   }
0155   cout << "    Calculated ratios." << endl;
0156 
0157   // set styles
0158   const UInt_t  fFil(0);
0159   const UInt_t  fLin(1);
0160   const UInt_t  fWid(1);
0161   const UInt_t  fTxt(42);
0162   const UInt_t  fAln(12);
0163   const UInt_t  fCnt(1);
0164   const Float_t fLab[NPad]  = {0.074, 0.04};
0165   const Float_t fTit[NPad]  = {0.074, 0.04};
0166   const Float_t fOffX[NPad] = {1.1, 1.};
0167   const Float_t fOffY[NPad] = {0.7, 1.3};
0168   for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0169     hDenom[iHist] -> SetMarkerColor(fColDen[iHist]);
0170     hDenom[iHist] -> SetMarkerStyle(fMarDen[iHist]);
0171     hDenom[iHist] -> SetFillColor(fColDen[iHist]);
0172     hDenom[iHist] -> SetFillStyle(fFil);
0173     hDenom[iHist] -> SetLineColor(fColDen[iHist]);
0174     hDenom[iHist] -> SetLineStyle(fLin);
0175     hDenom[iHist] -> SetLineWidth(fWid);
0176     hDenom[iHist] -> SetTitle(sTitle.Data());
0177     hDenom[iHist] -> SetTitleFont(fTxt);
0178     hDenom[iHist] -> GetXaxis() -> SetRangeUser(xPlotRange[0], xPlotRange[1]);
0179     hDenom[iHist] -> GetXaxis() -> SetTitle(sTitleX.Data());
0180     hDenom[iHist] -> GetXaxis() -> SetTitleFont(fTxt);
0181     hDenom[iHist] -> GetXaxis() -> SetTitleSize(fTit[1]);
0182     hDenom[iHist] -> GetXaxis() -> SetTitleOffset(fOffX[1]);
0183     hDenom[iHist] -> GetXaxis() -> SetLabelFont(fTxt);
0184     hDenom[iHist] -> GetXaxis() -> SetLabelSize(fLab[1]);
0185     hDenom[iHist] -> GetXaxis() -> CenterTitle(fCnt);
0186     hDenom[iHist] -> GetYaxis() -> SetTitle(sTitleY.Data());
0187     hDenom[iHist] -> GetYaxis() -> SetTitleFont(fTxt);
0188     hDenom[iHist] -> GetYaxis() -> SetTitleSize(fTit[1]);
0189     hDenom[iHist] -> GetYaxis() -> SetTitleOffset(fOffY[1]);
0190     hDenom[iHist] -> GetYaxis() -> SetLabelFont(fTxt);
0191     hDenom[iHist] -> GetYaxis() -> SetLabelSize(fLab[1]);
0192     hDenom[iHist] -> GetYaxis() -> CenterTitle(fCnt);
0193     hNumer[iHist] -> SetMarkerColor(fColNum[iHist]);
0194     hNumer[iHist] -> SetMarkerStyle(fMarNum[iHist]);
0195     hNumer[iHist] -> SetFillColor(fColNum[iHist]);
0196     hNumer[iHist] -> SetFillStyle(fFil);
0197     hNumer[iHist] -> SetLineColor(fColNum[iHist]);
0198     hNumer[iHist] -> SetLineStyle(fLin);
0199     hNumer[iHist] -> SetLineWidth(fWid);
0200     hNumer[iHist] -> SetTitle(sTitle.Data());
0201     hNumer[iHist] -> SetTitleFont(fTxt);
0202     hNumer[iHist] -> GetXaxis() -> SetRangeUser(xPlotRange[0], xPlotRange[1]);
0203     hNumer[iHist] -> GetXaxis() -> SetTitle(sTitleX.Data());
0204     hNumer[iHist] -> GetXaxis() -> SetTitleFont(fTxt);
0205     hNumer[iHist] -> GetXaxis() -> SetTitleSize(fTit[1]);
0206     hNumer[iHist] -> GetXaxis() -> SetTitleOffset(fOffX[1]);
0207     hNumer[iHist] -> GetXaxis() -> SetLabelFont(fTxt);
0208     hNumer[iHist] -> GetXaxis() -> SetLabelSize(fLab[1]);
0209     hNumer[iHist] -> GetXaxis() -> CenterTitle(fCnt);
0210     hNumer[iHist] -> GetYaxis() -> SetTitle(sTitleY.Data());
0211     hNumer[iHist] -> GetYaxis() -> SetTitleFont(fTxt);
0212     hNumer[iHist] -> GetYaxis() -> SetTitleSize(fTit[1]);
0213     hNumer[iHist] -> GetYaxis() -> SetTitleOffset(fOffY[1]);
0214     hNumer[iHist] -> GetYaxis() -> SetLabelFont(fTxt);
0215     hNumer[iHist] -> GetYaxis() -> SetLabelSize(fLab[1]);
0216     hNumer[iHist] -> GetYaxis() -> CenterTitle(fCnt);
0217     hRatio[iHist] -> SetMarkerColor(fColNum[iHist]);
0218     hRatio[iHist] -> SetMarkerStyle(fMarNum[iHist]);
0219     hRatio[iHist] -> SetFillColor(fColNum[iHist]);
0220     hRatio[iHist] -> SetFillStyle(fFil);
0221     hRatio[iHist] -> SetLineColor(fColNum[iHist]);
0222     hRatio[iHist] -> SetLineStyle(fLin);
0223     hRatio[iHist] -> SetLineWidth(fWid);
0224     hRatio[iHist] -> SetTitle(sTitle.Data());
0225     hRatio[iHist] -> SetTitleFont(fTxt);
0226     hRatio[iHist] -> GetXaxis() -> SetRangeUser(xPlotRange[0], xPlotRange[1]);
0227     hRatio[iHist] -> GetXaxis() -> SetTitle(sTitleX.Data());
0228     hRatio[iHist] -> GetXaxis() -> SetTitleFont(fTxt);
0229     hRatio[iHist] -> GetXaxis() -> SetTitleSize(fTit[0]);
0230     hRatio[iHist] -> GetXaxis() -> SetTitleOffset(fOffX[0]);
0231     hRatio[iHist] -> GetXaxis() -> SetLabelFont(fTxt);
0232     hRatio[iHist] -> GetXaxis() -> SetLabelSize(fLab[0]);
0233     hRatio[iHist] -> GetXaxis() -> CenterTitle(fCnt);
0234     hRatio[iHist] -> GetYaxis() -> SetTitle(sTitleR.Data());
0235     hRatio[iHist] -> GetYaxis() -> SetTitleFont(fTxt);
0236     hRatio[iHist] -> GetYaxis() -> SetTitleSize(fTit[0]);
0237     hRatio[iHist] -> GetYaxis() -> SetTitleOffset(fOffY[0]);
0238     hRatio[iHist] -> GetYaxis() -> SetLabelFont(fTxt);
0239     hRatio[iHist] -> GetYaxis() -> SetLabelSize(fLab[0]);
0240     hRatio[iHist] -> GetYaxis() -> CenterTitle(fCnt);
0241   }
0242   cout << "    Set styles." << endl;
0243 
0244   // make legend
0245   const UInt_t  fColLe(0);
0246   const UInt_t  fFilLe(0);
0247   const UInt_t  fLinLe(0);
0248   const UInt_t  nObjLe(2 * (NHist + 1));
0249   const Float_t hObjLe(nObjLe * 0.05);
0250   const Float_t yObjLe(0.1 + hObjLe);
0251   const Float_t fLegXY[NVtx] = {0.1, 0.1, 0.3, yObjLe};
0252 
0253   TLegend *leg = new TLegend(fLegXY[0], fLegXY[1], fLegXY[2], fLegXY[3]);
0254   leg -> SetFillColor(fColLe);
0255   leg -> SetFillStyle(fFilLe);
0256   leg -> SetLineColor(fColLe);
0257   leg -> SetLineStyle(fLinLe);
0258   leg -> SetTextFont(fTxt);
0259   leg -> SetTextAlign(fAln);
0260   leg -> AddEntry((TObject*)0, sHeadDenom.Data(), "");
0261   for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0262     leg -> AddEntry(hDenom[iHist], sLabelDenom[iHist], "pf");
0263   }
0264   leg -> AddEntry((TObject*)0, sHeadNumer.Data(), "");
0265   for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0266     leg -> AddEntry(hNumer[iHist], sLabelNumer[iHist], "pf");
0267   }
0268   cout << "    Made legend." << endl;
0269 
0270   // make text
0271   const UInt_t  fColTx(0);
0272   const UInt_t  fFilTx(0);
0273   const UInt_t  fLinTx(0);
0274   const UInt_t  nObjTx(NTxt);
0275   const Float_t hObjTx(nObjTx * 0.05);
0276   const Float_t yObjTx(0.1 + hObjTx);
0277   const Float_t fTxtXY[NVtx] = {0.3, 0.1, 0.5, yObjTx};
0278 
0279   TPaveText *txt = new TPaveText(fTxtXY[0], fTxtXY[1], fTxtXY[2], fTxtXY[3], "NDC NB");
0280   txt -> SetFillColor(fColTx);
0281   txt -> SetFillStyle(fFilTx);
0282   txt -> SetLineColor(fColTx);
0283   txt -> SetLineStyle(fLinTx);
0284   txt -> SetTextFont(fTxt);
0285   txt -> SetTextAlign(fAln);
0286   for (UInt_t iTxt = 0; iTxt < NTxt; iTxt++) {
0287     txt -> AddText(sText[iTxt].Data());
0288   }
0289   cout << "    Made text." << endl;
0290 
0291   // make line
0292   const UInt_t  fColLi(923);
0293   const UInt_t  fLinLi(9);
0294   const UInt_t  fWidLi(1);
0295   const Float_t fLinXY[NVtx] = {xPlotRange[0], 1., xPlotRange[1], 1.};
0296 
0297   TLine *line = new TLine(fLinXY[0], fLinXY[1], fLinXY[2], fLinXY[3]);
0298   line -> SetLineColor(fColLi);
0299   line -> SetLineStyle(fLinLi);
0300   line -> SetLineWidth(fWidLi);
0301   cout << "    Made line." << endl;
0302 
0303   // make plot
0304   const UInt_t  width(750);
0305   const UInt_t  height(950);
0306   const UInt_t  heightNR(750);
0307   const UInt_t  fMode(0);
0308   const UInt_t  fBord(2);
0309   const UInt_t  fGrid(0);
0310   const UInt_t  fTick(1);
0311   const UInt_t  fLogX(1);
0312   const UInt_t  fLogY1(1);
0313   const UInt_t  fLogY2(1);
0314   const UInt_t  fFrame(0);
0315   const Float_t fMarginL(0.15);
0316   const Float_t fMarginR(0.02);
0317   const Float_t fMarginT1(0.005);
0318   const Float_t fMarginT2(0.02);
0319   const Float_t fMarginB1(0.25);
0320   const Float_t fMarginB2(0.005);
0321   const Float_t fMarginBNR(0.15);
0322   const Float_t fPadXY1[NVtx] = {0., 0., 1., 0.35};
0323   const Float_t fPadXY2[NVtx] = {0., 0.35, 1., 1.};
0324 
0325   TCanvas *cPlot = new TCanvas("cPlot", "", width, height);
0326   TPad    *pPad1 = new TPad("pPad1", "", fPadXY1[0], fPadXY1[1], fPadXY1[2], fPadXY1[3]);
0327   TPad    *pPad2 = new TPad("pPad2", "", fPadXY2[0], fPadXY2[1], fPadXY2[2], fPadXY2[3]);
0328   cPlot     -> SetGrid(fGrid, fGrid);
0329   cPlot     -> SetTicks(fTick, fTick);
0330   cPlot     -> SetBorderMode(fMode);
0331   cPlot     -> SetBorderSize(fBord);
0332   pPad1     -> SetGrid(fGrid, fGrid);
0333   pPad1     -> SetTicks(fTick, fTick);
0334   pPad1     -> SetLogx(fLogX);
0335   pPad1     -> SetLogy(fLogY1);
0336   pPad1     -> SetBorderMode(fMode);
0337   pPad1     -> SetBorderSize(fBord);
0338   pPad1     -> SetFrameBorderMode(fFrame);
0339   pPad1     -> SetLeftMargin(fMarginL);
0340   pPad1     -> SetRightMargin(fMarginR);
0341   pPad1     -> SetTopMargin(fMarginT1);
0342   pPad1     -> SetBottomMargin(fMarginB1);
0343   pPad2     -> SetGrid(fGrid, fGrid);
0344   pPad2     -> SetTicks(fTick, fTick);
0345   pPad2     -> SetLogx(fLogX);
0346   pPad2     -> SetLogy(fLogY2);
0347   pPad2     -> SetBorderMode(fMode);
0348   pPad2     -> SetBorderSize(fBord);
0349   pPad2     -> SetFrameBorderMode(fFrame);
0350   pPad2     -> SetLeftMargin(fMarginL);
0351   pPad2     -> SetRightMargin(fMarginR);
0352   pPad2     -> SetTopMargin(fMarginT2);
0353   pPad2     -> SetBottomMargin(fMarginB2);
0354   cPlot     -> cd();
0355   pPad1     -> Draw();
0356   pPad2     -> Draw();
0357   pPad1     -> cd();
0358   hRatio[0] -> Draw(sOptRatio[0].Data());
0359   for (UInt_t iHist = 1; iHist < NHist; iHist++) {
0360     hRatio[iHist] -> Draw(sOptRatio[iHist].Data());
0361   }
0362   line      -> Draw();
0363   pPad2     -> cd();
0364   hDenom[0] -> Draw(sOptDenom[0].Data());
0365   hNumer[0] -> Draw(sOptNumer[0].Data());
0366   for(UInt_t iHist = 1; iHist < NHist; iHist++) {
0367     hDenom[iHist] -> Draw(sOptDenom[iHist].Data());
0368     hNumer[iHist] -> Draw(sOptNumer[iHist].Data());
0369   }
0370   leg     -> Draw();
0371   txt     -> Draw();
0372   fOutput -> cd();
0373   cPlot   -> Write();
0374   cPlot   -> Close();
0375 
0376   TCanvas *cPlotNR = new TCanvas("cPlot_NoRatio", "", width, heightNR);
0377   cPlotNR   -> SetGrid(fGrid, fGrid);
0378   cPlotNR   -> SetTicks(fTick, fTick);
0379   cPlotNR   -> SetLogx(fLogX);
0380   cPlotNR   -> SetLogy(fLogY2);
0381   cPlotNR   -> SetBorderMode(fMode);
0382   cPlotNR   -> SetBorderSize(fBord);
0383   cPlotNR   -> SetFrameBorderMode(fFrame);
0384   cPlotNR   -> SetLeftMargin(fMarginL);
0385   cPlotNR   -> SetRightMargin(fMarginR);
0386   cPlotNR   -> SetTopMargin(fMarginT1);
0387   cPlotNR   -> SetBottomMargin(fMarginBNR);
0388   hDenom[0] -> Draw(sOptDenom[0].Data());
0389   hNumer[0] -> Draw(sOptNumer[0].Data());
0390   for(UInt_t iHist = 1; iHist < NHist; iHist++) {
0391     hDenom[iHist] -> Draw(sOptDenom[iHist].Data());
0392     hNumer[iHist] -> Draw(sOptNumer[iHist].Data());
0393   }
0394   leg     -> Draw();
0395   txt     -> Draw();
0396   fOutput -> cd();
0397   cPlotNR -> Write();
0398   cPlotNR -> Close();
0399   cout << "    Made plot." << endl;
0400 
0401   // save histograms
0402   fOutput -> cd();
0403   for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0404     hDenom[iHist] -> Write();
0405     hNumer[iHist] -> Write();
0406     hRatio[iHist] -> Write();
0407   }
0408   cout << "    Saved histograms." << endl;
0409 
0410   // close files
0411   fOutput -> cd();
0412   fOutput -> Close();
0413   for (UInt_t iHist = 0; iHist < NHist; iHist++) {
0414     fDenom[iHist] -> cd();
0415     fDenom[iHist] -> Close();
0416     fNumer[iHist] -> cd();
0417     fNumer[iHist] -> Close();
0418   }
0419   cout << "  Finished plot!\n" << endl;
0420 
0421 }
0422 
0423 // end ------------------------------------------------------------------------