Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
draw_hists_3D.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file draw_hists_3D.C
1 #include "sPhenixStyle.C"
2 #include "sPhenixStyle.h"
3 
4 TPad* getPad(int j, int k, int l, float cw, float ch, const int nEBins, const int nCutBins);
5 
6 //void drawCanvas_invMass(TCanvas *c, TH1F *hInvMass, int pad_x, int pad_y, TPad *pad, int isEta);
7 void drawCanvas_invMass(TCanvas *c, TH1F *hCorr, int pad_x, int pad_y, TPad *pad, int isEta, float peakPos, float peakW);
8 
9 void SetHistoStyle(TH1F *histo, int cutBin, int eBin, float low, float hi, int isEta);
10 
11 void DrawCutLabel(int x, int y, int isEta);
12 
13 void GetSubtractedDist(TH1F *histOrig, TH1F *histSub, TF1 *invMassFit, TF1 *invMassBG);
14 
15 //float eBins[] = {1,3,5,7,9,11,13,15,17,20};
16 float eBins[] = {1,2,3,4,5,6,7,8,9,10,11,12,13};
17 float eBinsEta[] = {12,14,16,18,20};
18 float eCuts[] = {0.5,0.6,0.7,0.8,0.9,1.,1.1/*,1.2,1.3,1.4,1.5*/};
19 float eCutsEta[] = {3,3.5,4,4.5,5,5.5,6};
20 const int nEtaBins = 5;
21 void draw_hists_3D(const char *input)
22 {
23 
24  SetsPhenixStyle();
25  //TH1::SetDefaultSumw2();
26  //TH2::SetDefaultSumw2();
27  //TH3::SetDefaultSumw2();
28 
29  TFile *fin = new TFile(input);
30 
31  string th1List[] = {"photonE_Reco","isoPhotonE_Reco_Tru","truthPi0E","truthDphoE","pi0ETruthReco"};
32 
33  string th1xAxis[] = {"E_{#gamma}","E_{#gamma}^{Iso}","E_{#pi^{0}}^{Truth}","E_{#gamma_{Dir}}^{Truth}","E_{#pi^{0}}^{Reco}/E_{#pi^{0}}^{Truth}"};
34 
35  string th2List[] = {"clusterChi2","clusterProbPhoton","isoPhotonChi2","isoPhotonProb2","tspE","tspEiso","asymEtruthpi0","deltaREtruthpi0"};
36 
37  string th2xAxis[] = {"Cluster #xhi^{2}","Cluster Photon Prob","Iso Cluster #chi^{2}","Iso Cluster Prob","TSP","TSP","#frac{|E_{1} - E_{2}|}{E_{#pi^0}}","#Delta R [rad]"};
38 
39  string th2yAxis[] = {"E_{Cluster}","E_{Cluster}","E_{Cluster}","E_{Cluster}","E_{Cluster}","E_{Cluster}","E_{#pi^{0}}","E_{#pi^{0}}"};
40 
41  string th3List[] = {"deltaREinvMass","asymEinvMass","invMassEtaE","dphoChi2","dphoProb","pi0Chi2","pi0Prob","etaChi2","etaProb","eChi2","eProb","hChi2","hProb","ePi0InvMassEcut0"};
42 
43  string th3xAxis[] = {"#DeltaR","#frac{|E_{1} - E_{2}|}{E_{#pi^0}}","M_{#gamma#gamma}","#chi^{2}","Cluster Photon Prob","#chi^{2}","Cluster Photon Prob","#chi^{2}","Cluster Photon Prob","#chi^{2}","Cluster Photon Prob","#chi^{2}","Cluster Photon Prob","E_{#pi^{0}}^{Reco}"};
44 
45  string th3yAxis[] = {"E_{#pi^{0}}^{Reco}","E_{#pi^{0}}^{Reco}","#eta"," E_{Cluster}^{#gamma} [GeV]","E_{Cluster}^{#gamma}[GeV]","E_{Cluster}^{#pi^{0}}[GeV]","E_{Cluster}^{#pi^{0}}[GeV]","E_{Cluster}^{#eta}[GeV]","E_{Cluster}^{#eta}[GeV]","E_{Cluster}^{e^{#pm}}[GeV]","E_{Cluster}^{e^{#pm}}[GeV]","E_{Cluster}^{h^{#pm}}[GeV]","E_{Cluster}^{h^{#pm}}[GeV]","M_{#gamma#gamma}"};
46 
47  string th3zAxis[] = {"M_{#gamma#gamma}[GeV/c^{2}]","M_{#gamma#gamma}[GeV/c^{2}]","E^{#gamma}_{Reco}[GeV]","E^{#gamma}_{True}[GeV]","E^{#gamma}_{True}[GeV]","E^{#pi^{0}}_{True}[GeV]","E^{#pi^{0}}_{True}[GeV]","E^{#eta}_{True}[GeV]","E^{#eta}_{True}[GeV]","E^{e^{#pm}}_{True}[GeV]","E_{e^{#pm}}_{True}[GeV]","E_{#pi^{#pm}}_{True}[GeV]","E^{#pi^{#pm}}_{True}[GeV]","E_{Cut}[GeV]"};
48 
49 
50  //loop over and draw the TH1F's
51  /*TCanvas *c1 = new TCanvas();
52  const int nTh1 = (int)sizeof(th1List)/sizeof(th1List[0]);
53  for(int i = 0; i < nTh1; i++)
54  {
55  TH1F *hist = (TH1F*)fin -> Get(th1List[i].c_str());
56  hist -> GetXaxis() -> SetTitle(th1xAxis[i].c_str());
57  gPad -> SetLogy();
58  hist -> Draw();
59  c1 -> SaveAs(Form("plots/%s.pdf",th1List[i].c_str()));
60  }
61  gPad -> SetLogy(0);
62  const int nTh2 = (int)sizeof(th2List)/sizeof(th2List[0]);
63  for(int i = 0; i < nTh2; i++)
64  {
65  gPad -> SetLogz();
66  TH2F *hist = (TH2F*)fin -> Get(th2List[i].c_str());
67  hist -> GetYaxis() -> SetTitle(th2yAxis[i].c_str());
68  hist -> GetXaxis() -> SetTitle(th2xAxis[i].c_str());
69  hist -> Draw("colz");
70  c1 -> SaveAs(Form("plots/%s.pdf",th2List[i].c_str()));
71  }*/
72 
73 
74  const int nTh3 = (int)sizeof(th3List)/sizeof(th3List[0]);
75  const int nEBins = (int)sizeof(eBins)/sizeof(eBins[0]);
76  const int nCutBins = (int)sizeof(eCuts)/sizeof(eCuts[0]);
77  const int nCutBinsEta = (int)sizeof(eCutsEta)/sizeof(eCutsEta[0]);
78  const int nEBinsEta = (int)sizeof(eBinsEta)/sizeof(eBinsEta[0]);
79 
80  TGraphErrors *chi2E[5][2];
81  //[5]: dir pho, pi0, eta, electron, pi+/-
82  //[2]: vs. truth energy, vs. cluster energy
83  TGraphErrors *probE[5][2];
84 
85  /*
86  for(int i = 0; i < nTh3 - 1; i++)
87  {
88  gPad -> SetLogy(0);
89  TH3F *hist = (TH3F*)fin -> Get(th3List[i].c_str());
90  hist -> GetXaxis() -> SetTitle(th3xAxis[i].c_str());
91  hist -> GetYaxis() -> SetTitle(th3yAxis[i].c_str());
92  hist -> GetZaxis() -> SetTitle(th3zAxis[i].c_str());
93 
94  TH2F *hist2D = (TH2F*)hist -> Project3D("xy");
95  c1 -> cd();
96  hist2D -> Draw("colz");
97  c1 -> SaveAs(Form("plots/%s_projxy.pdf",th3List[i].c_str()));
98 
99  TCanvas *cSec = new TCanvas();
100 
101  if(i > 2 && i < nTh3 - 1)
102  {
103  for(int j = 0; j < nEBins - 1; j++)
104  {
105  TH1F *hist1D = (TH1F*)hist2D->ProjectionY(Form("%s_truthE_%g_%g",th3List[i].c_str(), eBins[j], eBins[j+1]), hist2D -> GetXaxis() -> FindBin(eBins[j]), hist2D -> GetXaxis() -> FindBin(eBins[j+1]));
106  //hist1D -> Fit("gaus",
107  //gPad -> WaitPrimitive();
108  hist1D -> Draw();
109  TLatex energy;
110  gPad -> SetLogy();
111  energy.DrawLatexNDC(0.5,0.7,Form("%g < %s < %g",eBins[j], th3yAxis[i].c_str(), eBins[j+1]));
112  cSec -> SaveAs(Form("plots/%s_ClusterE_%g_%g.pdf",th3List[i].c_str(), eBins[j], eBins[j+1]));
113 
114  }
115  }
116 
117  hist2D = (TH2F*)hist -> Project3D("xz");
118  c1 -> cd();
119  hist2D -> Draw("colz");
120  c1 -> SaveAs(Form("plots/%s_projxz.pdf",th3List[i].c_str()));
121 
122  if(i > 2 && i < nTh3 - 1)
123  {
124  for(int j = 0; j < nEBins - 1; j++)
125  {
126  TH1F *hist1D = (TH1F*)hist2D->ProjectionY(Form("%s_truthE_%g_%g",th3List[i].c_str(), eBins[j], eBins[j+1]), hist2D -> GetXaxis() -> FindBin(eBins[j]), hist2D -> GetXaxis() -> FindBin(eBins[j+1]));
127  cSec -> cd();
128  hist1D -> Draw();
129  TLatex energy;
130  gPad -> SetLogy();
131  energy.DrawLatexNDC(0.5,0.7,Form("%g < %s < %g",eBins[j], th3zAxis[i].c_str(), eBins[j+1]));
132  //hist1D -> Fit("gaus",
133  //gPad -> WaitPrimitive();
134  cSec -> SaveAs(Form("plots/%s_TruthE_%g_%g.pdf",th3List[i].c_str(), eBins[j], eBins[j+1]));
135 
136  }
137  }
138 
139  }*/
140 
141  //going to handle the last TH3 separately since
142  //it's going on a stupid big canvas.
143 
144  TH3F *hist[nEtaBins];
145  TH3F *histMatched[nEtaBins];
146  float cw = 2*966, ch = 1.3*637;
147 
148  TCanvas *cPi0Cut[nEtaBins];
149  TPad *pad_invMass[nEBins-1][nCutBins][nEtaBins];
150 
151  TCanvas *cSubMass[nEtaBins];
152  TPad *pad_invMassSub[nEBins-1][nCutBins][nEtaBins];
153 
154  TCanvas *cMatchMass[nEtaBins];
155  TPad *pad_matchMass[nEBins-1][nCutBins][nEtaBins];
156 
157 
158  TGraphErrors *mass[nCutBins][nEtaBins];
159  TGraphErrors *width[nCutBins][nEtaBins];
160  TGraphErrors *yields[nCutBins][nEtaBins];
161  TGraphErrors *yieldsMatched[nCutBins][nEtaBins];
162  TGraphErrors *matchZeroFrac[nCutBins][nEtaBins];
163  TF1 *invMassFit, *invMassFitBG;
164  int colors[] = {1,2,4,kGreen+2, kViolet,kCyan,kOrange+2,kMagenta+2,kAzure-2};
165  TLegend *leggy = new TLegend(.161, .729, .722, .939 );
166  leggy -> SetNColumns(3);
167  leggy -> SetFillStyle(0);
168  leggy -> SetBorderSize(0);
169 
170  float xEndFit = 0.25;
171  float xStartFit = 0.07;
172  float xEndHist = 0.25;
173  float xStartHist = 0.07;
174  float invMassWindow[2] = {0.112, 0.162};
175  float peakPos[nEBins - 1][nCutBins][nEtaBins];
176  float peakW[nEBins - 1][nCutBins][nEtaBins];
177  float sigma = 3.;
178  TLatex marks;
179 
180  for(int k = 0; k < nEtaBins; k++)
181  {
182  cMatchMass[k] = new TCanvas(Form("cMatchMass_Eta%d",k),"Mass from Matched Cluster",cw,ch);
183  cMatchMass[k] -> Range(0,0,1,1);
184 
185  cSubMass[k]= new TCanvas(Form("cMassSubPi0_Eta%d",k),"Subtracted Mass Dist",cw,ch);
186  cSubMass[k] -> Range(0,0,1,1);
187 
188  cPi0Cut[k] = new TCanvas(Form("cMassPi0_Eta%d",k),"Invariant Mass",cw,ch);
189  cPi0Cut[k] -> Range(0,0,1,1);
190 
191  hist[k] = (TH3F*)fin -> Get(Form("ePi0InvMassEcut_0Match_Eta%d",k));
192  histMatched[k] = (TH3F*)fin -> Get(Form("ePi0InvMassEcut_1Match_Eta%d",k));
193 
194  for(int i = 0; i < nEBins - 1; i++)
195  {
196  for(int j = 0; j < nCutBins; j++)
197  {
198  if(i == 0)
199  {
200  mass[j][k] = new TGraphErrors();
201  mass[j][k] -> SetMarkerColor(colors[j]);
202 
203  width[j][k] = new TGraphErrors();
204  width[j][k] -> SetMarkerColor(colors[j]);
205 
206  yields[j][k] = new TGraphErrors();
207  yields[j][k] -> SetMarkerColor(colors[j]);
208 
209  yieldsMatched[j][k] = new TGraphErrors();
210  yieldsMatched[j][k] -> SetMarkerColor(colors[j]);
211 
212  matchZeroFrac[j][k] = new TGraphErrors();
213  matchZeroFrac[j][k] -> SetMarkerColor(colors[j]);
214  }
215 
216  pad_invMass[i][j][k] = getPad(i,j,k,cw,ch,nEBins , nCutBins);//for unsubtracted invariant mass distributions
217  pad_invMass[i][j][k] -> cd();
218 
219  pad_invMassSub[i][j][k] = getPad(i,j,k,cw,ch,nEBins , nCutBins);//for subtracted mass distributions
220 
221  pad_matchMass[i][j][k] = getPad(i,j,k,cw,ch,nEBins, nCutBins);//for distributions composed of clusters matched to truth decay photons
222 
223  TH1F *hist1D = (TH1F*)hist[k] -> ProjectionY(Form("invMass_E%d_Cut%d_Eta%d",i,j,k), hist[k] -> GetXaxis() -> FindBin(eBins[i] + 0.001), hist[k] -> GetXaxis() -> FindBin(eBins[i+1] - 0.0001), hist[k] -> GetZaxis() -> FindBin(eCuts[j] + 0.0001), hist[k] -> GetNbinsZ());
224 
225  TH1F *hist1DMatched = (TH1F*)histMatched[k] -> ProjectionY(Form("invMass_E%d_Cut%d_eta%d_Matched",i,j,k), histMatched[k] -> GetXaxis() -> FindBin(eBins[i] + 0.001), histMatched[k] -> GetXaxis() -> FindBin(eBins[i+1] - 0.0001), histMatched[k] -> GetZaxis() -> FindBin(eCuts[j] + 0.0001), histMatched[k] -> GetNbinsZ());
226 
227 
228  TH1F *invMassSub = (TH1F*)hist1D -> Clone();
229  if(eBins[i+1] <= 13)
230  {
231  invMassFit = new TF1(Form("invMassFit_E%d_Cut%d_Eta%d",i,j,k),"[0] + [1]*x + [2]*pow(x,2)+ [3]*pow(x,3)+ [4]*pow(x,4)+ [8]*pow(x,6) +[5]*exp(-(([7]-x)^2)/(2*[6]^2))");
232 
233 
234  invMassFitBG = new TF1(Form("invMassFit_E%d_Cut%d_Eta%dBG",i,j,k),"[0] + [1]*x + [2]*pow(x,2)+ [3]*pow(x,3)+ [4]*pow(x,4)+ [8]*pow(x,6) +[5]*exp(-(([7]-x)^2)/(2*[6]^2)) - [5]*exp(-(([7]-x)^2)/(2*[6]^2))"); //hehe
235 
236  invMassFit -> FixParameter(8,0);
237 
238  invMassFit -> SetLineWidth(1);
239  invMassFit -> SetParameter(6,0.015);
240 
241  invMassFit -> SetLineColor(2);
242  invMassFit -> SetParameter(7,0.135);
243  invMassFit -> SetParLimits(7, 0.135 - 0.07, 0.135 + 0.07);
244  if((eCuts[j] >= 0.9) && (eBins[i] == 3 || eBins[i] == 4))
245  {
246  invMassFit -> SetParameter(5, hist1D -> GetMaximum()*2);
247 
248  }
249  else invMassFit -> SetParameter(5,hist1D -> GetMaximum());
250  hist1D -> Fit(invMassFit,"Q","",xStartFit, xEndFit);
251  GetSubtractedDist(hist1D, invMassSub, invMassFit, invMassFitBG);
252  invMassFitBG -> SetLineColor(4);
253 
254  float peakPosition = invMassFit -> GetParameter(7);
255  peakPos[i][j][k] = peakPosition;
256 
257  float peakWidth = invMassFit -> GetParameter(6);
258  peakW[i][j][k] = peakWidth;
259  Double_t yieldErr, yieldErrMatched;
260 
261  float yield;
262  if(hist1D -> Integral() == 0) yield = 0;
263  else yield = invMassSub -> IntegralAndError(invMassSub -> GetXaxis() -> FindBin((peakPosition - sigma*peakWidth)+0.0001), invMassSub -> GetXaxis() -> FindBin((peakPosition + sigma*peakWidth)+0.0001), yieldErr, "");
264 
265  float yieldMatched = hist1DMatched -> IntegralAndError(hist1DMatched -> FindFirstBinAbove(0), hist1DMatched -> GetNbinsX(), yieldErrMatched, "");
266  if(yield != 0)
267  {
268  mass[j][k] -> SetPoint(i, (eBins[i]+eBins[i+1])/2, invMassFit -> GetParameter(7));
269  mass[j][k] -> SetPointError(i, 0, invMassFit -> GetParError(7));
270 
271  width[j][k] -> SetPoint(i, (eBins[i]+eBins[i+1])/2, invMassFit -> GetParameter(6));
272  width[j][k] -> SetPointError(i, 0, invMassFit -> GetParError(6));
273 
274  yields[j][k] -> SetPoint(i, (eBins[i]+eBins[i+1])/2, yield);
275  yields[j][k] -> SetPointError(i, 0, yieldErr);
276 
277 
278  }
279  yieldsMatched[j][k] -> SetPoint(i, (eBins[i]+eBins[i+1])/2, yieldMatched);
280 
281  yieldsMatched[j][k] -> SetPointError(i, 0, yieldErrMatched);
282 
283  float zeroMatchYield = hist1DMatched -> GetBinContent(hist1DMatched -> FindFirstBinAbove(0));
284  float zeroMatchYieldErr = hist1DMatched -> GetBinError(46);
285  float matchYieldRatio = zeroMatchYield/yieldMatched;
286  float zeroMatchYieldRatErr = matchYieldRatio*sqrt(pow(yieldErrMatched/yieldMatched,2) + pow(zeroMatchYieldErr/zeroMatchYieldErr,2));
287 
288  matchZeroFrac[j][k] -> SetPoint(i, (eBins[i]+eBins[i+1])/2, matchYieldRatio);
289  matchZeroFrac[j][k] -> SetPointError(i, 0, zeroMatchYieldRatErr);
290  }
291  hist1D -> GetXaxis() -> SetRangeUser(xStartHist,xEndHist);
292 
293  TCanvas *cMassFitCheck = new TCanvas();
294  hist1D -> SetTitle(";M_{#gamma#gamma} [GeV/c^{2}];");
295  hist1D -> Draw("epx0");
296  invMassFitBG -> SetLineWidth(1);
297 
298  invMassFitBG -> Draw("same");
299 
300  marks.DrawLatexNDC(0.5,0.3,Form("#splitline{E_{Min} > %g GeV}{%g > E_{#pi^{0}} > %g}",eCuts[j], eBins[i], eBins[i+1]));
301  cMassFitCheck -> SaveAs(Form("plots/invMassFitCheck_%dE_%dCut_Sigma%g_Eta%d.pdf",i,j,sigma,k));
302 
303 
304  TCanvas *cMassMatchCheck = new TCanvas();
305  hist1DMatched -> SetTitle(";M_{#gamma#gamma} [GeV/c^{2}];");
306  hist1DMatched -> Draw("epx0");
307  gPad -> SetLogy();
308  marks.DrawLatexNDC(0.5,0.85,Form("#splitline{E_{Min} > %g GeV}{%g > E_{#pi^{0}} > %g}",eCuts[j], eBins[i], eBins[i+1]));
309  cMassMatchCheck -> SaveAs(Form("plots/invMassMatchCheck_%dE_%dCut_Sigma%g_Eta%dMatch.pdf",i,j,sigma,k));
310 
311  drawCanvas_invMass(cPi0Cut[k], hist1D, j, i, pad_invMass[i][j][k], 0, peakPos[i][j][k], sigma*peakW[i][j][k]);
312  invMassFitBG -> Draw("same");
313  pad_invMassSub[i][j][k] -> cd();
314 
315  invMassSub -> GetXaxis() -> SetRangeUser(xStartHist, xEndHist);
316  drawCanvas_invMass(cSubMass[k], invMassSub, j, i, pad_invMassSub[i][j][k], 0, peakPos[i][j][k], sigma*peakW[i][j][k]);
317 
318  pad_matchMass[i][j][k] -> cd();
319  drawCanvas_invMass(cMatchMass[k], hist1DMatched, j, i, pad_matchMass[i][j][k], 1, peakPos[i][j][k], sigma*peakW[i][j][k]);
320 
321  }
322  }
323 
324  cPi0Cut[k] -> SaveAs(Form("plots/invMassCuts_Sigma%g_Eta%d.pdf",sigma, k));
325  cSubMass[k] -> SaveAs(Form("plots/invMassCutsSubtracted_Sigma%g_Eta%d.pdf",sigma, k));
326  cMatchMass[k] -> SaveAs(Form("plots/invMassMatched_Sigma%g_Eta%d.pdf",sigma,k));
327  }
328 
329 
330  TCanvas *cMass[nEtaBins];
331  for(int k = 0; k < nEtaBins; k++)
332  {
333  cMass[k] = new TCanvas();
334  for(int i = 0; i < nCutBins; i++)
335  {
336  if(i == 0)
337  {
338  mass[i][k] -> Draw("ap");
339  mass[i][k] -> GetYaxis() -> SetRangeUser(0.130,0.180);
340  mass[i][k] -> GetXaxis() -> SetLimits(0,13);
341  mass[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
342  mass[i][k] -> GetYaxis() -> SetTitle("Mass Peak Location [GeV/c^{2}]");
343  }
344  else mass[i][k] -> Draw("samep");
345  if(k==0)leggy -> AddEntry(mass[i][k], Form("E_{min} > %g",eCuts[i]),"p");
346  }
347 
348  leggy -> Draw("same");
349  cMass[k] -> SaveAs(Form("plots/invMassGraphPeakPos_Sigma%g_Eta%d.pdf",sigma,k));
350  }
351 
352  TCanvas *cWidth[nEtaBins];
353  for(int k = 0; k < nEtaBins; k++)
354  {
355  cWidth[k] = new TCanvas();
356 
357  for(int i = 0; i < nCutBins; i++)
358  {
359  if(i == 0)
360  {
361  width[i][k] -> Draw("ap");
362  width[i][k] -> GetYaxis() -> SetRangeUser(0.013, 0.028);
363  width[i][k] -> GetXaxis() -> SetLimits(0,13);
364  width[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
365  width[i][k] -> GetYaxis() -> SetTitle("Mass Peak Width [GeV/c^{2}]");
366  }
367  else width[i][k] -> Draw("samep");
368  }
369 
370  leggy -> Draw("same");
371  cWidth[k] -> SaveAs(Form("plots/invMassGraphPeakWidth_Sigma%g_Eta%d.pdf",sigma,k));
372  }
373 
374  TCanvas *cYields[nEtaBins];
375  for(int k = 0; k < nEtaBins; k++)
376  {
377  cYields[k] = new TCanvas();
378  for(int i = 0; i < nCutBins; i++)
379  {
380  if(i == 0)
381  {
382  yields[i][k] -> Draw("ap");
383  yields[i][k] -> GetYaxis() -> SetRangeUser(1e3,1e7);
384  yields[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
385  yields[i][k] -> GetYaxis() -> SetTitle("N_{#pi^{0}}");
386  }
387  else yields[i][k] -> Draw("samep");
388  gPad -> SetLogy();
389  }
390 
391  leggy -> Draw("same");
392  cYields[k] -> SaveAs(Form("plots/invMassYield_Sigma%g_Eta%d.pdf",sigma,k));
393  }
394 
395  TCanvas *cYieldsMatch[nEtaBins];
396  for(int k = 0; k < nEtaBins; k++)
397  {
398  cYieldsMatch[k] = new TCanvas();
399  for(int i = 0; i < nCutBins; i++)
400  {
401  if(i == 0)
402  {
403  yieldsMatched[i][k] -> Draw("ap");
404  yieldsMatched[i][k] -> GetYaxis() -> SetRangeUser(1e3,1e7);
405  yieldsMatched[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
406  yieldsMatched[i][k] -> GetYaxis() -> SetTitle("N_{#pi^{0}}");
407  }
408  else yields[i][k] -> Draw("samep");
409  gPad -> SetLogy();
410  }
411 
412  leggy -> Draw("same");
413  cYields[k] -> SaveAs(Form("plots/invMassYield_Sigma%g_Eta%d_Matched.pdf",sigma,k));
414  }
415  TH1F *pi0TruSpec[nEtaBins];
416  //for(int i = 0; i < nEtaBins; i++) pi0TruSpec[i] = (TH1F*)fin -> Get(th1List[2].c_str());
417  for(int i = 0; i < nEtaBins; i++) pi0TruSpec[i] = (TH1F*)fin -> Get(Form("truthPi0E_Eta%d",i));
418  TGraphErrors *pi0Eff[nCutBins][nEtaBins];
419  TCanvas *cEff[nEtaBins];
420  float scaleMax;
421  float scaleMin;
422 
423  for(int k = 0; k < nEtaBins; k++)
424  {
425  cEff[k] = new TCanvas();
426  scaleMax = -9999;
427  scaleMin = 9999;
428  for(int i = 0; i < nCutBins ; i++)
429  {
430  pi0Eff[i][k] = new TGraphErrors();
431  for(int j = 0; j < nEBins; j++)
432  {
433  Double_t pi0TruthYieldErr;
434  float pi0TruthYield = pi0TruSpec[k] -> IntegralAndError(pi0TruSpec[k] -> GetXaxis() -> FindBin(eBins[j]+0.001), pi0TruSpec[k] -> GetXaxis() -> FindBin(eBins[j+1]-0.001), pi0TruthYieldErr, "");
435 
436 
437  Double_t x, pi0YieldMeas, pi0YieldErr;
438  yields[i][k] -> GetPoint(j, x, pi0YieldMeas);
439  pi0YieldErr = yields[i][k] -> GetErrorY(j);
440  float eff = pi0YieldMeas/pi0TruthYield;
441  float err = eff*sqrt(pow(pi0YieldErr/pi0YieldMeas,2) + pow(pi0TruthYieldErr/pi0TruthYield,2));
442  if(eff < scaleMin) scaleMin = eff;
443  if(eff > scaleMax) scaleMax = eff;
444  pi0Eff[i][k] -> SetPoint(j, (eBins[j]+eBins[j+1])/2, eff);
445  pi0Eff[i][k] -> SetPointError(j, 0, err);
446  }
447  pi0Eff[i][k] -> SetMarkerColor(colors[i]);
448  if(i == 0)
449  {
450  pi0Eff[i][k] -> Draw("ap");
451  pi0Eff[i][k] -> GetYaxis() -> SetRangeUser(-0.01,scaleMax*1.4);
452  pi0Eff[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
453  pi0Eff[i][k] -> GetYaxis() -> SetTitle("#epsilon");
454  }
455  else pi0Eff[i][k] -> Draw("samep");
456 
457  }
458  leggy -> Draw("same");
459  cEff[k] -> SaveAs(Form("plots/pi0Eff_Sigma%g_Eta%d.pdf",sigma,k));
460  }
461 
462  TGraphErrors *pi0EffMatch[nCutBins][nEtaBins];
463  TCanvas *cEffMatch[nEtaBins];
464 
465  for(int k = 0; k < nEtaBins; k++)
466  {
467  scaleMax = -9999;
468  scaleMin = 9999;
469  cEffMatch[k] = new TCanvas();
470  for(int i = 0; i < nCutBins ; i++)
471  {
472  pi0EffMatch[i][k] = new TGraphErrors();
473  for(int j = 0; j < nEBins - 1; j++)
474  {
475  Double_t pi0TruthYieldErr;
476  float pi0TruthYield = pi0TruSpec[k] -> IntegralAndError(pi0TruSpec[k] -> GetXaxis() -> FindBin(eBins[j]+0.001), pi0TruSpec [k]-> GetXaxis() -> FindBin(eBins[j+1]-0.001), pi0TruthYieldErr, "");
477 
478  Double_t x, pi0YieldMeas, pi0YieldErr;
479  yieldsMatched[i][k] -> GetPoint(j, x, pi0YieldMeas);
480  pi0YieldErr = yieldsMatched[i][k] -> GetErrorY(j);
481  float eff = pi0YieldMeas/pi0TruthYield;
482  float err = eff*sqrt(pow(pi0YieldErr/pi0YieldMeas,2) + pow(pi0TruthYieldErr/pi0TruthYield,2));
483  if(eff < scaleMin) scaleMin = eff;
484  if(eff > scaleMax) scaleMax = eff;
485  pi0EffMatch[i][k] -> SetPoint(j, (eBins[j]+eBins[j+1])/2, eff);
486  pi0EffMatch[i][k] -> SetPointError(j, 0, err);
487  }
488  pi0EffMatch[i][k] -> SetMarkerColor(colors[i]);
489  if(i == 0)
490  {
491  pi0EffMatch[i][k] -> Draw("ap");
492  pi0EffMatch[i][k] -> GetYaxis() -> SetRangeUser(-0.01,scaleMax*1.25);
493  pi0EffMatch[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
494  pi0EffMatch[i][k] -> GetYaxis() -> SetTitle("#epsilon");
495  }
496  else pi0EffMatch[i][k] -> Draw("samep");
497 
498  }
499  leggy -> Draw("same");
500  cEffMatch[k] -> SaveAs(Form("plots/pi0Eff_Sigma%g_Eta%d_Matched.pdf",sigma,k));
501  }
502 
503  TCanvas *cZeroFrac[nEtaBins];
504  for(int k = 0; k < nEtaBins; k++)
505  {
506  cZeroFrac[k] = new TCanvas();
507  for(int j = 0; j < nCutBins; j++)
508  {
509  if(j == 0)
510  {
511  matchZeroFrac[j][k] -> Draw("ap");
512  //matchZeroFrac[j][k] -> GetYaxis() -> SetRangeUser(-0.1,1.1);
513  matchZeroFrac[j][k] -> SetTitle(";E_{#pi^{0}}^{Truth} [GeV]; N^{Unmatched}/N^{Matched}");
514  }
515  else matchZeroFrac[j][k] -> Draw("samep");
516  }
517  leggy -> Draw("same");
518  cZeroFrac[k] -> SaveAs(Form("plots/matchedZeroYieldFrac_Sigma%g_Eta%d.pdf",sigma,k));
519  }
520 
521 
522 
523  TCanvas *cMiss = new TCanvas();
524  cMiss -> Divide(4);
525  cMiss -> cd();
526 
527  TH3F *unmatchedLocale = (TH3F*)fin -> Get("unmatchedLocale");
528 
529  float eBinsMerge[] = {1,4,5,9,13};
530 
531  for(int i = 0; i < 4; i++)
532  {
533  unmatchedLocale -> GetZaxis() -> SetRangeUser(eBinsMerge[i], eBinsMerge[i+1]);
534  TH2F *locationProj = (TH2F*)unmatchedLocale -> Project3D("yx");
535  locationProj -> SetTitle(Form("%g < E^{#pi^{0}}_{Truth} < %g;#eta;#phi [rad]", eBinsMerge[i], eBinsMerge[i+1]));
536  locationProj -> Draw("colz");
537  marks.DrawLatexNDC(0.2,0.5,Form("%g < E^{#pi^{0}}_{Truth} < %g", eBinsMerge[i], eBinsMerge[i+1]));
538 
539  cMiss -> SaveAs(Form("plots/missLocationBin_%d.pdf",i));
540  }
541 
542  TCanvas *cMassScale = new TCanvas();
543  TH3F *invMassPhoPi0 = (TH3F*)fin -> Get("invMassPhoPi0");
544 
545  float massZones[] = {0,0.100,0.200,1};
546 
547  for(int i = 0; i < 3; i++)
548  {
549  TH1F *pi0EScale = (TH1F*)invMassPhoPi0 -> ProjectionZ(Form("pi0EScale_%d",i),invMassPhoPi0 -> GetXaxis() -> FindBin(massZones[i] + 0.0001), invMassPhoPi0 -> GetXaxis() -> FindBin(massZones[i + 1] - 0.0001), 1, invMassPhoPi0 -> GetNbinsY());
550  pi0EScale -> SetTitle(Form("%g < M_{#gamma#gamma} < %g; E_{Reco}^{#pi^{0}}/E_{True}^{#pi^{0}};",massZones[i], massZones[i+1]));
551  pi0EScale -> Draw();
552  marks.DrawLatexNDC(0.7,0.85,Form("%g < M_{#gamma#gamma} < %g",massZones[i], massZones[i+1]));
553  cMassScale -> SaveAs(Form("plots/pi0EScale_MassBin_%d.pdf",i));
554 
555 
556  TH1F *phoEScale = (TH1F*)invMassPhoPi0 -> ProjectionY(Form("phoEScale_%d",i),invMassPhoPi0 -> GetXaxis() -> FindBin(massZones[i] + 0.0001), invMassPhoPi0 -> GetXaxis() -> FindBin(massZones[i + 1] - 0.0001), 1, invMassPhoPi0 -> GetNbinsZ());
557  phoEScale -> SetTitle(Form("%g < M_{#gamma#gamma} < %g; E_{Reco}^{#gamma}/E_{True}^{#gamma};",massZones[i], massZones[i+1]));
558  phoEScale -> Draw();
559  marks.DrawLatexNDC(0.7,0.85,Form("%g < M_{#gamma#gamma} < %g",massZones[i], massZones[i+1]));
560 
561  cMassScale -> SaveAs(Form("plots/phoEScale_MassBin_%d.pdf",i));
562  }
563 
564  TCanvas *cMassPhi = new TCanvas();
565  TH3F *invMassEPhi = (TH3F*)fin -> Get("invMassEPhi");
566 
567  for(int i = 0; i < 4; i++)
568  {
569  invMassEPhi -> GetYaxis() -> SetRangeUser(eBinsMerge[i],eBinsMerge[i+1]);
570  TH2F *invMassPhi = (TH2F*)invMassEPhi -> Project3D("xz");
571 
572  invMassPhi -> SetTitle(Form("%g < E^{#pi^{0}}_{Truth} < %g; #phi [rad]; M_{#gamma#gamma};",eBinsMerge[i], eBinsMerge[i+1] ));
573  invMassPhi -> Draw("colz");
574  marks.DrawLatexNDC(0.2,0.18,Form("%g < E^{#pi^{0}}_{Truth} < %g",eBinsMerge[i], eBinsMerge[i+1] ));
575 
576  cMassPhi -> SaveAs(Form("plots/invMassPhi_EBin_%d.pdf",i));
577  }
578 
579 
580 
581 
582  /*
583  TCanvas *cEtaCut = new TCanvas("cMassEta","Invariant Mass Eta Cuts",cw,ch);
584  cEtaCut -> Range(0,0,1,1);
585  TPad *pad_invMassEta[nEBinsEta - 1][nCutBinsEta];
586 
587  for(int i = 0; i < nEBinsEta - 1; i++)
588  {
589  for(int j = 0; j < nCutBinsEta; j++)
590  {
591  pad_invMassEta[i][j] = getPad(i,j,cw,ch,nEBinsEta , nCutBinsEta);
592  pad_invMassEta[i][j] -> cd();
593  TH1F *hist1D = (TH1F*)hist -> ProjectionY(Form("invMass_EtaE%d_Cut%d",i,j), hist -> GetXaxis() -> FindBin(eBinsEta[i] + 0.001), hist -> GetXaxis() -> FindBin(eBinsEta[i+1] - 0.0001), hist -> GetZaxis() -> FindBin(eCutsEta[j] + 0.0001), hist -> GetNbinsZ());
594  drawCanvas_invMass(cEtaCut, hist1D, j, i, pad_invMassEta[i][j], 2);
595  }
596  }
597 
598  cEtaCut -> SaveAs("plots/invMassCutsEta.pdf");*/
599 
600 }
601 
602 TPad* getPad(int j, int k, int l, float cw, float ch, const int nEBins, const int nCutBins)
603 {
604  float xpos[nCutBins+1], ypos[nEBins];
605  int nCutOffset = 1;
606  float marginLeft = 0.07;
607  float marginRight = 0.035;
608 
609  float marginTop = 0.03;
610  float marginBottom = 0.06;
611 
612  float marginColumn = 0.04;
613 
614  float w = (1 - marginLeft - marginRight - marginColumn*(nCutBins-2 + nCutOffset))/(nCutBins - 1+ nCutOffset);
615  float h = (1 - marginTop - marginBottom)/(nEBins - 1);
616 
617  for(int m = 0; m < nCutBins + nCutOffset; m++)
618  {
619  if(m == 0) xpos[m] = 0;
620  else if(m == nCutBins - 1+ nCutOffset) xpos[m] = 1.;
621  else xpos[m] = marginLeft + m*w + (m-1)*marginColumn;
622  }
623 
624  for(int m = 0; m < nEBins; m++)
625  {
626  if(m==0) ypos[nEBins -1 -m] = 0;
627  else if (m == nEBins - 1) ypos[nEBins-1-m] = 1.;
628  else ypos[nEBins - 1 -m] = marginBottom +m*h;
629  }
630 
631  TPad *myPad = new TPad(Form("pad%d%d%d",j,k,l),Form("pad%d%d%d",j,k,l),xpos[k],ypos[j+1],xpos[k+1],ypos[j]);
632 
633  if(k==0) myPad->SetLeftMargin(marginLeft/(xpos[k+1] - xpos[k]));
634  else myPad -> SetLeftMargin(marginColumn/(xpos[k+1] - xpos[k]));
635 
636  if(k == nCutBins - 2+ nCutOffset) myPad -> SetRightMargin(marginRight/(xpos[k+1]-xpos[k]));
637  else myPad -> SetRightMargin(0.012);
638 
639  if(j == 0) myPad -> SetTopMargin(marginTop/(ypos[j] - ypos[j+1]));
640  else myPad -> SetTopMargin(0);
641 
642  if(j == nEBins - 2) myPad -> SetBottomMargin(marginBottom/(ypos[j] - ypos[j+1]));
643  else myPad -> SetBottomMargin(0);
644 
645  return myPad;
646 }
647 
648 void drawCanvas_invMass(TCanvas *c, TH1F *hCorr, int pad_x, int pad_y, TPad *pad, int mode, float peakPos, float peakW)
649 {
650 
651  float range = hCorr -> GetMaximum() - hCorr -> GetMinimum();
652  SetHistoStyle(hCorr, pad_x, pad_y, hCorr -> GetMinimum(), hCorr->GetMaximum() + 0.2*range, mode);
653  c -> cd(0);
654  /*TLine *pi0Mass = new TLine(0.135,0,0.135,hCorr->GetMaximum());
655  pi0Mass -> SetLineStyle(7);
656  pi0Mass -> SetLineColor(2);
657  TLine *etaMass = new TLine(0.548,0,0.548,hCorr->GetMaximum());
658  etaMass -> SetLineStyle(7);
659  etaMass -> SetLineColor(2);*/
660  TLine *massWindow[2];
661  //massWindow[0] = new TLine(0.112,0,0.112,hCorr->GetMaximum());
662  pad -> Draw();
663  pad -> cd();
664 
665  hCorr -> Draw("epX0");
666  DrawCutLabel(pad_x,pad_y,mode);
667  if(mode == 1) gPad -> SetLogy();
668  if(peakPos)
669  {
670  massWindow[0] = new TLine(peakPos - peakW,0,peakPos - peakW,hCorr->GetMaximum());
671  massWindow[0] -> SetLineStyle(7);
672  massWindow[0] -> SetLineColor(2);
673 
674  //massWindow[1] = new TLine(0.162,0,0.162,hCorr->GetMaximum());
675  massWindow[1] = new TLine(peakPos + peakW,0,peakPos + peakW,hCorr->GetMaximum());
676  massWindow[1] -> SetLineStyle(7);
677  massWindow[1] -> SetLineColor(2);
678  massWindow[0] -> Draw("same");
679  massWindow[1] -> Draw("same");
680  }
681 
682 
683  //etaMass -> Draw("same");
684  //pi0Mass -> Draw("same");
685 
686 }
687 
688 void SetHistoStyle(TH1F *histo, int cutBin, int eBin, float low, float hi, int mode)
689 {
690  if(hi > low) histo -> GetYaxis() -> SetRangeUser(low,hi);
691  if(mode == 1)
692  {
693  histo -> GetYaxis() -> SetRangeUser(0.1,10*hi);
694  }
695 
696  histo -> SetTitle(";;");
697  int xTitleBin[3][2] = {{11,3},{11,3},{3,3}};
698  int yTitleBin[3][2] = {{5,0},{5,0},{2,0}};
699 
700  float yTitleOffset[3] = {0.3, 0.4, 1};
701  float xTitleOffset[3] = {0.7, 0.7, 0.7};
702 
703  float yTitleSize[3] = {0.5, 0.5, 0.15};
704  float xTitleSize[3] = {0.28, 0.28, 0.13};
705 
706  if(eBin == xTitleBin[mode][0] && cutBin == xTitleBin[mode][1]) histo -> GetXaxis() -> SetTitle("M_{#gamma#gamma}");
707  if(eBin == yTitleBin[mode][0] && cutBin == yTitleBin[mode][1]) histo -> GetYaxis() -> SetTitle("N_{#gamma#gamma}");
708 
709  histo -> GetXaxis() -> SetTitleSize(xTitleSize[mode]);
710  histo -> GetXaxis() -> SetTitleOffset(xTitleOffset[mode]);
711  histo -> GetXaxis() -> CenterTitle();
712 
713  histo -> GetXaxis() -> SetLabelFont(43);
714  histo -> GetXaxis() -> SetLabelSize(18);
715  histo -> GetXaxis() -> SetLabelOffset(0.0);
716  histo -> SetTickLength(0.05,"X");
717  histo -> SetNdivisions(505,"x");
718 
719 
720  histo -> GetYaxis() -> SetNdivisions(505);
721  histo -> GetYaxis() -> SetTitleSize(yTitleSize[mode]);
722  histo -> GetYaxis() -> SetTitleOffset(yTitleOffset[mode]);
723  histo -> GetYaxis() -> CenterTitle();
724 
725  histo -> GetYaxis() -> SetLabelFont(43);
726  histo -> GetYaxis() -> SetLabelSize(19);
727  histo -> GetYaxis() -> SetLabelOffset(0.003);
728  histo -> GetYaxis() -> SetNoExponent();
729  histo -> SetTickLength(0.04,"Y");
730 
731  histo -> SetStats(0);
732  //histo -> SetLineWidth(1.4);
733  histo -> SetLineColor(1);
734  histo -> SetMarkerColor(1);
735  histo -> SetMarkerStyle(20);
736  histo -> SetMarkerSize(0.4);
737 }
738 
739 void DrawCutLabel(int x, int y, int mode)
740 {
741  gPad -> cd();
742  gPad -> Update();
743  float Y2 = gPad -> GetUymax();
744  float Y1 = gPad -> GetUymin();
745  float rangeY = Y2 - Y1;
746  float xtop[3] = {0.8, -0.1, 0};
747  float xbottom[3] = {0.273152, 1.01, 1.1};
748  double ytop[3] = {Y2+0.01*rangeY, Y2+1*rangeY, Y2+0.01*rangeY};
749  double ybottom[3] = { 0.98*Y2, 0.98*Y2, 0.9*Y2};
750 
751  TText *textCut;
752  if(y == 0)
753  {
754  // if(mode < 2)
755  // {
756  // if(x <= 4)textCut = new TText(0.08, Y2+0.01*rangeY, Form("Min E > %g GeV", eCuts[x]));
757  // else textCut = new TText(0.08, Y2+0.01*rangeY, Form("Min E > %g GeV", eCuts[x]));
758  // }
759  //else textCut = new TText(0, Y2+0.01*rangeY, Form("Min E > %g GeV", eCutsEta[x]));
760  textCut = new TText(xtop[mode], ytop[mode], Form("Min E > %g GeV", eCuts[x]));
761  textCut -> SetTextFont(43);
762  textCut -> SetTextSize(18);
763  textCut -> Draw();
764  }
765  if(x == sizeof(eCuts)/sizeof(eCuts[0]) - 1)
766  {
767  //if(mode < 2)textCut = new TText(0.273152, 0.98*Y2, Form("%g-%gGeV",eBins[y], eBins[y+1]));
768  //else textCut = new TText(1.1, 0.9*Y2, Form("%g-%gGeV",eBinsEta[y], eBinsEta[y+1]));
769  if(mode < 2)textCut = new TText(xbottom[mode], ybottom[mode], Form("%g-%gGeV",eBins[y], eBins[y+1]));
770  else textCut = new TText(xbottom[mode], ybottom[mode], Form("%g-%gGeV",eBinsEta[y], eBinsEta[y+1]));
771  textCut -> SetTextAngle(270);
772  textCut -> SetTextFont(43);
773  textCut -> SetTextSize(18);
774  textCut -> Draw();
775  }
776 
777 }
778 
779 void GetSubtractedDist(TH1F *histOrig, TH1F *histSub, TF1 *invMassFit, TF1 *invMassBG)
780 {
781 
782 
783  for(int i = 0; i < invMassFit -> GetNpar(); i++)
784  {
785  invMassBG -> SetParameter(i, invMassFit -> GetParameter(i));
786  }
787 
788  for(int i = 1; i < histSub -> GetNbinsX() + 1; i++)
789  {
790  histSub -> SetBinContent(i, histOrig -> GetBinContent(i) - invMassBG -> Eval(histOrig -> GetBinCenter(i)));
791  }
792 
793 
794 }