Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MakeRatioComparisonPlot.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MakeRatioComparisonPlot.C
1 // ----------------------------------------------------------------------------
2 // 'MakeRatioComparisonPlot.C'
3 // Derek Anderson
4 // 01.03.2022
5 //
6 // Use this quickly plot a set of numerator
7 // distributions against a set of denominator
8 // distributions
9 // ----------------------------------------------------------------------------
10 
11 #include <iostream>
12 #include "TH1.h"
13 #include "TPad.h"
14 #include "TFile.h"
15 #include "TLine.h"
16 #include "TError.h"
17 #include "TString.h"
18 #include "TLegend.h"
19 #include "TCanvas.h"
20 #include "TPaveText.h"
21 
22 using namespace std;
23 
24 // global constants
25 static const UInt_t NHist(2);
26 static const UInt_t NPlot(2);
27 static const UInt_t NPad(2);
28 static const UInt_t NVtx(4);
29 static const UInt_t NTxt(3);
30 
31 
32 
33 void MakeRatioComparisonPlot() {
34 
35  // lower verbosity
36  gErrorIgnoreLevel = kError;
37  cout << "\n Beginning plot macro..." << endl;
38 
39  // file parameters
40  const TString sOutput("debugCorrelators.alexVsF4A_truthVsRecoAfterPhiFix_ptJet30.pp200run6jet40.d18m4y2023.root");
41  const TString sInDenom[NHist] = {"reference/alex_output/pp200run6jet40.SphenixTruthUnscaled_output_trksAndChrgPars.root",
42  "reference/alex_code/SphenixRecoUnscaled_output_derekTestRun.root"};
43  const TString sInNumer[NHist] = {"./output/pp200run6jet40.forDebug_true_withPhiFix.d18m4y2023.root",
44  "./output/pp200run6jet40.forDebug_reco_withPhiFix.d18m4y2023.root"};
45 
46  // denominator parameters
47  const TString sHeadDenom("#bf{Reference Output}");
48  const TString sHistDenom[NHist] = {"EEC2", "EEC2"};
49  const TString sNameDenom[NHist] = {"hReferencTruth", "hReferenceReco"};
50  const TString sLabelDenom[NHist] = {"Truth",
51  "Reco."};
52 
53  // numerator parameters
54  const TString sHeadNumer("#bf{F4A Implementation}");
55  const TString sHistNumer[NHist] = {"hCorrelatorDrAxis_ptJet30", "hCorrelatorDrAxis_ptJet30"};
56  const TString sNameNumer[NHist] = {"hFun4AllTruth", "hFun4AllReco"};
57  const TString sLabelNumer[NHist] = {"Truth",
58  "Reco."};
59 
60  // plot parameters
61  const TString sTitle("");
62  const TString sTitleX("#DeltaR");
63  const TString sTitleY("EEC");
64  const TString sTitleR("[f4a] / [reference]");
65  const TString sNameRatio[NHist] = {"hRatioTruth", "hRatioReco"};
66  const TString sOptDenom[NHist] = {"", "sames"};
67  const TString sOptNumer[NHist] = {"sames", "sames"};
68  const TString sOptRatio[NHist] = {"", "sames"};
69  const TString sText[NTxt] = {"#bf{#it{sPHENIX}} Simulation [Run6]", "JS 40 GeV Jet Sample", "p_{T}^{jet} #in (30, 50) GeV"};
70  const Float_t xPlotRange[NPlot] = {0.001, 1.};
71  const UInt_t fColDen[NHist] = {923, 634};
72  const UInt_t fColNum[NHist] = {921, 895};
73  const UInt_t fMarDen[NHist] = {20, 29};
74  const UInt_t fMarNum[NHist] = {24, 30};
75 
76  // norm/rebin parameters
77  const Bool_t doIntNorm(false);
78  const Bool_t doRebinDenom[NHist] = {false, false};
79  const Bool_t doRebinNumer[NHist] = {false, false};
80  const UInt_t nRebinDenom[NHist] = {2, 2};
81  const UInt_t nRebinNumer[NHist] = {2, 2};
82 
83  // open files
84  TFile *fOutput = new TFile(sOutput.Data(), "recreate");
85  if (!fOutput) {
86  cerr << "PANIC: couldn't open output file!\n" << endl;
87  return;
88  }
89 
90  TFile *fDenom[NHist];
91  TFile *fNumer[NHist];
92  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
93  fDenom[iHist] = new TFile(sInDenom[iHist].Data(), "read");
94  fNumer[iHist] = new TFile(sInNumer[iHist].Data(), "read");
95  if (!fDenom[iHist] || !fNumer[iHist]) {
96  cerr << "PANIC: couldn't open denominator or numerator file #" << iHist << "!\n"
97  << " fDenom = " << fDenom[iHist] << ", fNumer = " << fNumer[iHist] << "\n"
98  << endl;
99  return;
100  }
101  }
102  cout << " Opened files." << endl;
103 
104  // grab histograms
105  TH1D *hDenom[NHist];
106  TH1D *hNumer[NHist];
107  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
108  hDenom[iHist] = (TH1D*) fDenom[iHist] -> Get(sHistDenom[iHist].Data());
109  hNumer[iHist] = (TH1D*) fNumer[iHist] -> Get(sHistNumer[iHist].Data());
110  if (!hDenom[iHist] || !hNumer[iHist]) {
111  cerr << "PANIC: couldn't grab numerator or denominator histogram # " << iHist << "!\n"
112  << " hDenom = " << hDenom[iHist] << ", hNumer = " << hNumer[iHist] << "\n"
113  << endl;
114  return;
115  }
116  hDenom[iHist] -> SetName(sNameDenom[iHist].Data());
117  hNumer[iHist] -> SetName(sNameNumer[iHist].Data());
118  }
119  cout << " Grabbed histograms." << endl;
120 
121  // rebin histograms
122  Bool_t doRebin(false);
123  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
124  doRebin = (doRebinDenom[iHist] || doRebinNumer[iHist]);
125  if (doRebin) break;
126  }
127 
128  if (doRebin) {
129  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
130  if (doRebinDenom[iHist]) hDenom[iHist] -> Rebin(nRebinDenom[iHist]);
131  if (doRebinNumer[iHist]) hNumer[iHist] -> Rebin(nRebinNumer[iHist]);
132  }
133  cout << " Rebinned histograms." << endl;
134  }
135 
136  // normalize by integrals )if needed)
137  if (doIntNorm) {
138  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
139  const Double_t intDenom = hDenom[iHist] -> Integral();
140  const Double_t intNumer = hNumer[iHist] -> Integral();
141  hDenom[iHist] -> Scale(1. / intDenom);
142  hNumer[iHist] -> Scale(1. / intNumer);
143  }
144  cout << " Normalized histograms by integral." << endl;
145  }
146 
147  // calculate ratios
148  TH1D *hRatio[NHist];
149  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
150  hRatio[iHist] = (TH1D*) hDenom[iHist] -> Clone();
151  hRatio[iHist] -> Reset("ICE");
152  hRatio[iHist] -> Divide(hNumer[iHist], hDenom[iHist], 1., 1.);
153  hRatio[iHist] -> SetName(sNameRatio[iHist]);
154  }
155  cout << " Calculated ratios." << endl;
156 
157  // set styles
158  const UInt_t fFil(0);
159  const UInt_t fLin(1);
160  const UInt_t fWid(1);
161  const UInt_t fTxt(42);
162  const UInt_t fAln(12);
163  const UInt_t fCnt(1);
164  const Float_t fLab[NPad] = {0.074, 0.04};
165  const Float_t fTit[NPad] = {0.074, 0.04};
166  const Float_t fOffX[NPad] = {1.1, 1.};
167  const Float_t fOffY[NPad] = {0.7, 1.3};
168  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
169  hDenom[iHist] -> SetMarkerColor(fColDen[iHist]);
170  hDenom[iHist] -> SetMarkerStyle(fMarDen[iHist]);
171  hDenom[iHist] -> SetFillColor(fColDen[iHist]);
172  hDenom[iHist] -> SetFillStyle(fFil);
173  hDenom[iHist] -> SetLineColor(fColDen[iHist]);
174  hDenom[iHist] -> SetLineStyle(fLin);
175  hDenom[iHist] -> SetLineWidth(fWid);
176  hDenom[iHist] -> SetTitle(sTitle.Data());
177  hDenom[iHist] -> SetTitleFont(fTxt);
178  hDenom[iHist] -> GetXaxis() -> SetRangeUser(xPlotRange[0], xPlotRange[1]);
179  hDenom[iHist] -> GetXaxis() -> SetTitle(sTitleX.Data());
180  hDenom[iHist] -> GetXaxis() -> SetTitleFont(fTxt);
181  hDenom[iHist] -> GetXaxis() -> SetTitleSize(fTit[1]);
182  hDenom[iHist] -> GetXaxis() -> SetTitleOffset(fOffX[1]);
183  hDenom[iHist] -> GetXaxis() -> SetLabelFont(fTxt);
184  hDenom[iHist] -> GetXaxis() -> SetLabelSize(fLab[1]);
185  hDenom[iHist] -> GetXaxis() -> CenterTitle(fCnt);
186  hDenom[iHist] -> GetYaxis() -> SetTitle(sTitleY.Data());
187  hDenom[iHist] -> GetYaxis() -> SetTitleFont(fTxt);
188  hDenom[iHist] -> GetYaxis() -> SetTitleSize(fTit[1]);
189  hDenom[iHist] -> GetYaxis() -> SetTitleOffset(fOffY[1]);
190  hDenom[iHist] -> GetYaxis() -> SetLabelFont(fTxt);
191  hDenom[iHist] -> GetYaxis() -> SetLabelSize(fLab[1]);
192  hDenom[iHist] -> GetYaxis() -> CenterTitle(fCnt);
193  hNumer[iHist] -> SetMarkerColor(fColNum[iHist]);
194  hNumer[iHist] -> SetMarkerStyle(fMarNum[iHist]);
195  hNumer[iHist] -> SetFillColor(fColNum[iHist]);
196  hNumer[iHist] -> SetFillStyle(fFil);
197  hNumer[iHist] -> SetLineColor(fColNum[iHist]);
198  hNumer[iHist] -> SetLineStyle(fLin);
199  hNumer[iHist] -> SetLineWidth(fWid);
200  hNumer[iHist] -> SetTitle(sTitle.Data());
201  hNumer[iHist] -> SetTitleFont(fTxt);
202  hNumer[iHist] -> GetXaxis() -> SetRangeUser(xPlotRange[0], xPlotRange[1]);
203  hNumer[iHist] -> GetXaxis() -> SetTitle(sTitleX.Data());
204  hNumer[iHist] -> GetXaxis() -> SetTitleFont(fTxt);
205  hNumer[iHist] -> GetXaxis() -> SetTitleSize(fTit[1]);
206  hNumer[iHist] -> GetXaxis() -> SetTitleOffset(fOffX[1]);
207  hNumer[iHist] -> GetXaxis() -> SetLabelFont(fTxt);
208  hNumer[iHist] -> GetXaxis() -> SetLabelSize(fLab[1]);
209  hNumer[iHist] -> GetXaxis() -> CenterTitle(fCnt);
210  hNumer[iHist] -> GetYaxis() -> SetTitle(sTitleY.Data());
211  hNumer[iHist] -> GetYaxis() -> SetTitleFont(fTxt);
212  hNumer[iHist] -> GetYaxis() -> SetTitleSize(fTit[1]);
213  hNumer[iHist] -> GetYaxis() -> SetTitleOffset(fOffY[1]);
214  hNumer[iHist] -> GetYaxis() -> SetLabelFont(fTxt);
215  hNumer[iHist] -> GetYaxis() -> SetLabelSize(fLab[1]);
216  hNumer[iHist] -> GetYaxis() -> CenterTitle(fCnt);
217  hRatio[iHist] -> SetMarkerColor(fColNum[iHist]);
218  hRatio[iHist] -> SetMarkerStyle(fMarNum[iHist]);
219  hRatio[iHist] -> SetFillColor(fColNum[iHist]);
220  hRatio[iHist] -> SetFillStyle(fFil);
221  hRatio[iHist] -> SetLineColor(fColNum[iHist]);
222  hRatio[iHist] -> SetLineStyle(fLin);
223  hRatio[iHist] -> SetLineWidth(fWid);
224  hRatio[iHist] -> SetTitle(sTitle.Data());
225  hRatio[iHist] -> SetTitleFont(fTxt);
226  hRatio[iHist] -> GetXaxis() -> SetRangeUser(xPlotRange[0], xPlotRange[1]);
227  hRatio[iHist] -> GetXaxis() -> SetTitle(sTitleX.Data());
228  hRatio[iHist] -> GetXaxis() -> SetTitleFont(fTxt);
229  hRatio[iHist] -> GetXaxis() -> SetTitleSize(fTit[0]);
230  hRatio[iHist] -> GetXaxis() -> SetTitleOffset(fOffX[0]);
231  hRatio[iHist] -> GetXaxis() -> SetLabelFont(fTxt);
232  hRatio[iHist] -> GetXaxis() -> SetLabelSize(fLab[0]);
233  hRatio[iHist] -> GetXaxis() -> CenterTitle(fCnt);
234  hRatio[iHist] -> GetYaxis() -> SetTitle(sTitleR.Data());
235  hRatio[iHist] -> GetYaxis() -> SetTitleFont(fTxt);
236  hRatio[iHist] -> GetYaxis() -> SetTitleSize(fTit[0]);
237  hRatio[iHist] -> GetYaxis() -> SetTitleOffset(fOffY[0]);
238  hRatio[iHist] -> GetYaxis() -> SetLabelFont(fTxt);
239  hRatio[iHist] -> GetYaxis() -> SetLabelSize(fLab[0]);
240  hRatio[iHist] -> GetYaxis() -> CenterTitle(fCnt);
241  }
242  cout << " Set styles." << endl;
243 
244  // make legend
245  const UInt_t fColLe(0);
246  const UInt_t fFilLe(0);
247  const UInt_t fLinLe(0);
248  const UInt_t nObjLe(2 * (NHist + 1));
249  const Float_t hObjLe(nObjLe * 0.05);
250  const Float_t yObjLe(0.1 + hObjLe);
251  const Float_t fLegXY[NVtx] = {0.1, 0.1, 0.3, yObjLe};
252 
253  TLegend *leg = new TLegend(fLegXY[0], fLegXY[1], fLegXY[2], fLegXY[3]);
254  leg -> SetFillColor(fColLe);
255  leg -> SetFillStyle(fFilLe);
256  leg -> SetLineColor(fColLe);
257  leg -> SetLineStyle(fLinLe);
258  leg -> SetTextFont(fTxt);
259  leg -> SetTextAlign(fAln);
260  leg -> AddEntry((TObject*)0, sHeadDenom.Data(), "");
261  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
262  leg -> AddEntry(hDenom[iHist], sLabelDenom[iHist], "pf");
263  }
264  leg -> AddEntry((TObject*)0, sHeadNumer.Data(), "");
265  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
266  leg -> AddEntry(hNumer[iHist], sLabelNumer[iHist], "pf");
267  }
268  cout << " Made legend." << endl;
269 
270  // make text
271  const UInt_t fColTx(0);
272  const UInt_t fFilTx(0);
273  const UInt_t fLinTx(0);
274  const UInt_t nObjTx(NTxt);
275  const Float_t hObjTx(nObjTx * 0.05);
276  const Float_t yObjTx(0.1 + hObjTx);
277  const Float_t fTxtXY[NVtx] = {0.3, 0.1, 0.5, yObjTx};
278 
279  TPaveText *txt = new TPaveText(fTxtXY[0], fTxtXY[1], fTxtXY[2], fTxtXY[3], "NDC NB");
280  txt -> SetFillColor(fColTx);
281  txt -> SetFillStyle(fFilTx);
282  txt -> SetLineColor(fColTx);
283  txt -> SetLineStyle(fLinTx);
284  txt -> SetTextFont(fTxt);
285  txt -> SetTextAlign(fAln);
286  for (UInt_t iTxt = 0; iTxt < NTxt; iTxt++) {
287  txt -> AddText(sText[iTxt].Data());
288  }
289  cout << " Made text." << endl;
290 
291  // make line
292  const UInt_t fColLi(923);
293  const UInt_t fLinLi(9);
294  const UInt_t fWidLi(1);
295  const Float_t fLinXY[NVtx] = {xPlotRange[0], 1., xPlotRange[1], 1.};
296 
297  TLine *line = new TLine(fLinXY[0], fLinXY[1], fLinXY[2], fLinXY[3]);
298  line -> SetLineColor(fColLi);
299  line -> SetLineStyle(fLinLi);
300  line -> SetLineWidth(fWidLi);
301  cout << " Made line." << endl;
302 
303  // make plot
304  const UInt_t width(750);
305  const UInt_t height(950);
306  const UInt_t heightNR(750);
307  const UInt_t fMode(0);
308  const UInt_t fBord(2);
309  const UInt_t fGrid(0);
310  const UInt_t fTick(1);
311  const UInt_t fLogX(1);
312  const UInt_t fLogY1(1);
313  const UInt_t fLogY2(1);
314  const UInt_t fFrame(0);
315  const Float_t fMarginL(0.15);
316  const Float_t fMarginR(0.02);
317  const Float_t fMarginT1(0.005);
318  const Float_t fMarginT2(0.02);
319  const Float_t fMarginB1(0.25);
320  const Float_t fMarginB2(0.005);
321  const Float_t fMarginBNR(0.15);
322  const Float_t fPadXY1[NVtx] = {0., 0., 1., 0.35};
323  const Float_t fPadXY2[NVtx] = {0., 0.35, 1., 1.};
324 
325  TCanvas *cPlot = new TCanvas("cPlot", "", width, height);
326  TPad *pPad1 = new TPad("pPad1", "", fPadXY1[0], fPadXY1[1], fPadXY1[2], fPadXY1[3]);
327  TPad *pPad2 = new TPad("pPad2", "", fPadXY2[0], fPadXY2[1], fPadXY2[2], fPadXY2[3]);
328  cPlot -> SetGrid(fGrid, fGrid);
329  cPlot -> SetTicks(fTick, fTick);
330  cPlot -> SetBorderMode(fMode);
331  cPlot -> SetBorderSize(fBord);
332  pPad1 -> SetGrid(fGrid, fGrid);
333  pPad1 -> SetTicks(fTick, fTick);
334  pPad1 -> SetLogx(fLogX);
335  pPad1 -> SetLogy(fLogY1);
336  pPad1 -> SetBorderMode(fMode);
337  pPad1 -> SetBorderSize(fBord);
338  pPad1 -> SetFrameBorderMode(fFrame);
339  pPad1 -> SetLeftMargin(fMarginL);
340  pPad1 -> SetRightMargin(fMarginR);
341  pPad1 -> SetTopMargin(fMarginT1);
342  pPad1 -> SetBottomMargin(fMarginB1);
343  pPad2 -> SetGrid(fGrid, fGrid);
344  pPad2 -> SetTicks(fTick, fTick);
345  pPad2 -> SetLogx(fLogX);
346  pPad2 -> SetLogy(fLogY2);
347  pPad2 -> SetBorderMode(fMode);
348  pPad2 -> SetBorderSize(fBord);
349  pPad2 -> SetFrameBorderMode(fFrame);
350  pPad2 -> SetLeftMargin(fMarginL);
351  pPad2 -> SetRightMargin(fMarginR);
352  pPad2 -> SetTopMargin(fMarginT2);
353  pPad2 -> SetBottomMargin(fMarginB2);
354  cPlot -> cd();
355  pPad1 -> Draw();
356  pPad2 -> Draw();
357  pPad1 -> cd();
358  hRatio[0] -> Draw(sOptRatio[0].Data());
359  for (UInt_t iHist = 1; iHist < NHist; iHist++) {
360  hRatio[iHist] -> Draw(sOptRatio[iHist].Data());
361  }
362  line -> Draw();
363  pPad2 -> cd();
364  hDenom[0] -> Draw(sOptDenom[0].Data());
365  hNumer[0] -> Draw(sOptNumer[0].Data());
366  for(UInt_t iHist = 1; iHist < NHist; iHist++) {
367  hDenom[iHist] -> Draw(sOptDenom[iHist].Data());
368  hNumer[iHist] -> Draw(sOptNumer[iHist].Data());
369  }
370  leg -> Draw();
371  txt -> Draw();
372  fOutput -> cd();
373  cPlot -> Write();
374  cPlot -> Close();
375 
376  TCanvas *cPlotNR = new TCanvas("cPlot_NoRatio", "", width, heightNR);
377  cPlotNR -> SetGrid(fGrid, fGrid);
378  cPlotNR -> SetTicks(fTick, fTick);
379  cPlotNR -> SetLogx(fLogX);
380  cPlotNR -> SetLogy(fLogY2);
381  cPlotNR -> SetBorderMode(fMode);
382  cPlotNR -> SetBorderSize(fBord);
383  cPlotNR -> SetFrameBorderMode(fFrame);
384  cPlotNR -> SetLeftMargin(fMarginL);
385  cPlotNR -> SetRightMargin(fMarginR);
386  cPlotNR -> SetTopMargin(fMarginT1);
387  cPlotNR -> SetBottomMargin(fMarginBNR);
388  hDenom[0] -> Draw(sOptDenom[0].Data());
389  hNumer[0] -> Draw(sOptNumer[0].Data());
390  for(UInt_t iHist = 1; iHist < NHist; iHist++) {
391  hDenom[iHist] -> Draw(sOptDenom[iHist].Data());
392  hNumer[iHist] -> Draw(sOptNumer[iHist].Data());
393  }
394  leg -> Draw();
395  txt -> Draw();
396  fOutput -> cd();
397  cPlotNR -> Write();
398  cPlotNR -> Close();
399  cout << " Made plot." << endl;
400 
401  // save histograms
402  fOutput -> cd();
403  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
404  hDenom[iHist] -> Write();
405  hNumer[iHist] -> Write();
406  hRatio[iHist] -> Write();
407  }
408  cout << " Saved histograms." << endl;
409 
410  // close files
411  fOutput -> cd();
412  fOutput -> Close();
413  for (UInt_t iHist = 0; iHist < NHist; iHist++) {
414  fDenom[iHist] -> cd();
415  fDenom[iHist] -> Close();
416  fNumer[iHist] -> cd();
417  fNumer[iHist] -> Close();
418  }
419  cout << " Finished plot!\n" << endl;
420 
421 }
422 
423 // end ------------------------------------------------------------------------