simple-tof-analysis
 All Classes Namespaces Functions Variables Groups Pages
Plotting.h
1 #include <TBox.h>
2 #include <TPaveStats.h>
3 #include <TList.h>
4 #include <TLatex.h>
5 #include <TPaveText.h>
6 #include <TPad.h>
7 #include <TAttPad.h>
8 
9 #include "FitSpectra.h"
10 
11 
12 using namespace std;
13 
14 
15 
16 const std::vector<Color_t> vParticleColor = {kOrange, kMagenta, kRed, kBlue, kBlack, kYellow, kGreen+3, kAzure+10};
17 
18 
19 
20 
21 //__________________________________________________________________________________________________________________________
22 void RedrawBorder(bool isYLinear = false)
23 {
24  gPad->Update();
25  gPad->RedrawAxis();
26  TLine l;
27  l.SetLineWidth(3);
28 
29  double ymin = pow(10, gPad->GetUymin());
30  double ymax = pow(10, gPad->GetUymax());
31 
32  if (isYLinear == true)
33  {
34  ymin = gPad->GetUymin();
35  ymax = gPad->GetUymax();
36  }
37 
38  l.DrawLine(gPad->GetUxmin(), ymin, gPad->GetUxmax(), ymin); // x1, y1, x2, y2 //bottom line
39  l.DrawLine(gPad->GetUxmin(), ymin, gPad->GetUxmin(), ymax); // x1, y1, x2, y2 //left line
40  l.DrawLine(gPad->GetUxmin(), ymax, gPad->GetUxmax(), ymax); // x1, y1, x2, y2 //top line
41  l.DrawLine(gPad->GetUxmax(), ymin, gPad->GetUxmax(), ymax); // x1, y1, x2, y2 //right line
42 }
43 
44 
45 
46 //__________________________________________________________________________________________________________________________
47 void PrintDxDy (TH2D * h, const string name, const string year)
48 {
49  TCanvas * c1 = new TCanvas("c1", "c1", 1200, 1200); //1800, 1200
50 
51  gStyle->SetOptStat(0); // set stat box to show number of entries
52  gPad->SetGridx(1);
53  gPad->SetGridy(1);
54 
55  c1->SetLogz(1);
56  c1->SetTopMargin(0.04);
57  c1->SetBottomMargin(0.09);
58  c1->SetRightMargin(0.15);
59  c1->SetLeftMargin(0.09);
60 
61  h->GetZaxis()->SetTitle("Counts");
62  h->GetXaxis()->SetTitleOffset(1.2);
63  h->GetYaxis()->SetTitleOffset(1.1);
64  h->GetZaxis()->SetTitleOffset(1.3);
65 
66  h->SetAxisRange(-6, 6, "X");
67  h->SetAxisRange(-6, 6, "Y");
68 
69  //h->SetMinimum(1);
70  //h->SetMaximum(zmax);
71 
72  h->Draw("colz");
73  c1->Update();
74 
75  string filename = name + h->GetName() + "_" + year + ".png";
76  c1->SaveAs(filename.c_str());
77  delete c1;
78 }
79 
80 
81 
82 
83 //__________________________________________________________________________________________________________________________
84 void Print1DMassSqFits (const TH1D * h1, vector<TH1D*> vFits, const string name, const string year, bool titleOn = false)
85 {
86  if (h1 == NULL)
87  {
89  cerr<<"[WARNING]: Specified histogram not found!"<<endl;
90  return;
91  }
92 
93  TH1D * h = (TH1D*)h1->Clone();
94 
95  int rebin = 1;
96  double xMin = -2.5;
97  gStyle->SetOptStat(0); // set stat box to show number of entries
98 
99  TCanvas * c1 = new TCanvas("c1", "c1", 1800, 1200); //1800, 1200
100  TPad * pad0 = new TPad("pad0", "The top pad with 70% of the height", 0, 0.3, 1, 1); //xlow, ylow, xup, yup,
101  TPad * pad1 = new TPad("pad1", "The bottom pad with 30% of the height", 0, 0, 1, 0.3);
102  pad0->Draw();
103  pad1->Draw();
104 
105  pad0->cd();
106  pad0->SetLogy(1);
107  pad0->SetTicky(1);
108  pad0->SetTickx(1);
109  pad0->SetGridx(1);
110  pad0->SetGridy(1);
111  pad0->SetTopMargin(0.03);
112  pad0->SetBottomMargin(0.005);
113  pad0->SetRightMargin(0.02);
114  pad0->SetLeftMargin(0.1);
115 
116  if (titleOn == true)
117  {
118  pad0->SetTopMargin(0.08);
119  gStyle->SetTitleX(0.5);
120  gStyle->SetTitleAlign(23);
121  gStyle->SetTitleBorderSize(0);
122  }
123 
124  h->Rebin(rebin);
125  h->GetXaxis()->SetTitleOffset(1.3);
126  h->GetYaxis()->SetTitle("Counts");
127  h->GetYaxis()->SetTitleOffset(0.65);
128  h->GetYaxis()->SetTitleSize(0.07);
129  h->SetLineWidth(3);
130  h->SetMinimum(0.2);
131  h->SetAxisRange(xMin, 4.9, "X");
132  h->Draw();
133 
134  TLegend * leg = new TLegend(0.15, 0.6, 0.4, 0.9); //x1, y1, x2, y2
135  // leg->SetNColumns(1);
136  // leg->SetColumnSeparation(0.1);
137  // leg->SetEntrySeparation(0.05);
138  leg->SetFillStyle(0);
139  leg->SetBorderSize(0);
140  leg->AddEntry(h, "Data", "LE");
141 
142  vector<string> legLabel = {"Fit total", "Fit kaon", "Fit proton", "Fit deuteron", "Fit pion", "Fit electron"};
143  vector<Color_t> colors = {kRed, kGreen, kBlack, kAzure+10, kBlue, kMagenta};
144  for (int k = int(vFits.size())-1; k >= 0; k--)
145  {
146  if (k != 0)
147  vFits[k]->SetLineStyle(9);
148 
149  vFits[k]->Rebin(rebin);
150  vFits[k]->SetLineColor(colors[k]);
151  vFits[k]->SetLineWidth(3);
152  vFits[k]->Draw("SAME");
153  leg->AddEntry(vFits[k], legLabel[k].c_str(), "l"); //"p" "l";
154  }
155 
156  leg->Draw();
157  RedrawBorder();
158  pad0->Update();
159 
160  pad1->cd();
161  pad1->SetLogy(1);
162  pad1->SetTickx(1);
163  pad1->SetTicky(1);
164  pad1->SetGridx(1);
165  pad1->SetGridy(1);
166  pad1->SetTopMargin(0.005);
167  pad1->SetBottomMargin(0.35);
168  pad1->SetRightMargin(0.02);
169  pad1->SetLeftMargin(0.1);
170 
171  TH1D * hRatio = GetRatioTH1D(h, vFits[0]);
172  hRatio->SetTitle("");
173  hRatio->SetLineColor(kBlack);
174  hRatio->SetLineWidth(3);
175  hRatio->SetMarkerSize(5);
176  hRatio->SetMinimum(0.1);
177  hRatio->SetMaximum(5);
178  hRatio->SetAxisRange(xMin, 4.9, "X");
179  hRatio->GetXaxis()->SetTitleOffset(1.0);
180  hRatio->GetXaxis()->SetTitleSize(0.16);
181  hRatio->GetXaxis()->SetLabelSize(0.12);
182  hRatio->GetYaxis()->SetTitle("Data/Model");
183  hRatio->GetYaxis()->SetTitleOffset(0.3);
184  hRatio->GetYaxis()->SetTitleSize(0.15);
185  hRatio->GetYaxis()->SetLabelSize(0.08);
186  hRatio->Draw("E1 C P X0");
187 
188 
189  TLine * line = new TLine(xMin, 1, 5.0, 1); //x1, y1, x2, y2
190  line->SetLineColor(kRed);
191  line->SetLineWidth(3);
192  // line->SetLineStyle(2);
193  line->Draw("same");
194 
195  RedrawBorder(); //RedrawBorder(true);
196  pad1->Update();
197 
198  c1->Update();
199  string filename = name + h->GetName() + "_" + year + ".png";
200  c1->SaveAs(filename.c_str());
201 
202  delete h;
203  delete c1;
204 }
205 
206 
207 
208 
209 
210 //__________________________________________________________________________________________________________________________
211 void Print1DMassSqFitsNoRatio (const TH1D * h1,
212  vector<TH1D*> vFits,
213  const string name,
214  const string year,
215  bool titleOn = false,
216  double d_norm_PercentErr = 0,
217  double chiSquaredperNDF = 0)
218 {
219  if (h1 == NULL)
220  {
222  cerr<<"[WARNING]: Specified histogram not found!"<<endl;
223  return;
224  }
225 
226  TH1D * h = (TH1D*)h1->Clone();
227 
228  int rebin = 1;
229  double xMin = -2.5;
230  gStyle->SetOptStat(0); // set stat box to show number of entries
231 
232  TCanvas * c1 = new TCanvas("c1", "c1", 1200, 1200); //1800, 1200
233  c1->SetLogy(1);
234  c1->SetTicky(1);
235  c1->SetTickx(1);
236  c1->SetGridx(1);
237  c1->SetGridy(1);
238  c1->SetTopMargin(0.03);
239  c1->SetBottomMargin(0.12);
240  c1->SetRightMargin(0.02);
241  c1->SetLeftMargin(0.11);
242 
243  if (titleOn == true)
244  {
245  c1->SetTopMargin(0.08);
246  gStyle->SetTitleX(0.55);
247  gStyle->SetTitleAlign(23);
248  gStyle->SetTitleBorderSize(0);
249 
251  //gStyle->SetTitleFontSize(0.08);
252  // gStyle->SetTitleSize(0.1, "t");
253  gStyle->SetTitleW(0.9);
254  }
255 
256  h->Rebin(rebin);
257  h->GetXaxis()->SetTitle("TOF m^{2} (GeV/c^{2})^{2}");
258  h->GetXaxis()->SetTitleOffset(0.70);
259  h->GetXaxis()->SetTitleSize(0.07);
260  h->GetYaxis()->SetTitle("Counts");
261  h->GetYaxis()->SetTitleOffset(0.75);
262  h->GetYaxis()->SetTitleSize(0.07);
263  h->SetMinimum(0.2);
264  h->SetAxisRange(xMin, 4.9, "X");
265  h->SetLineWidth(3);
266  h->SetLineColor(kBlack);
267  h->SetMarkerStyle(20);
268  h->SetMarkerSize(2);
269  h->Draw("P E1 X0");
270 
271  TLegend * leg = new TLegend(0.12, 0.60, 0.4, 0.90); //x1, y1, x2, y2
272  // leg->SetNColumns(1);
273  // leg->SetColumnSeparation(0.1);
274  // leg->SetEntrySeparation(0.05);
275  leg->SetFillStyle(0);
276  leg->SetBorderSize(0);
277  leg->AddEntry(h, "Data", "PE");
278 
279  vector<string> legLabel = {"Fit total", "Fit kaon", "Fit proton", "Fit deuteron", "Fit pion", "Fit electron"};
280  vector<Color_t> colors = {kRed, kMagenta, kBlue, kAzure+10, kYellow, kGreen};
281  for (int k = int(vFits.size())-1; k >= 0; k--)
282  {
283  if (k != 0)
284  vFits[k]->SetLineStyle(9);
285 
286  if (k == 0)
287  vFits[k]->SetLineWidth(5);
288  else
289  vFits[k]->SetLineWidth(3);
290 
291  vFits[k]->Rebin(rebin);
292  vFits[k]->SetLineColor(colors[k]);
293  vFits[k]->Draw("SAME");
294  leg->AddEntry(vFits[k], legLabel[k].c_str(), "l"); //"p" "l";
295  }
296 
297  leg->Draw();
298  RedrawBorder();
299 
301  TPaveText * pt1 = NULL;
302  if (vFits[3]->GetMaximum() > 5)
303  {
304  pt1 = new TPaveText(0.71, 0.40, 0.95, 0.60, "NB NDC"); //x1, y1, x2, y2
305  pt1->SetTextSize(0.06);
306  pt1->SetTextColor(kAzure+10);
307  pt1->SetFillColor(0);
308  pt1->SetFillStyle(0);
309  pt1->SetBorderSize(0);
310  // pt1->SetTextAlign(12);
311  pt1->AddText("Deuteron");
312  pt1->Draw();
313  }
314 
316  TPaveText * pt2 = new TPaveText(0.60, 0.77, 0.7, 0.90, "NB NDC"); //x1, y1, x2, y2
317  pt2->SetTextSize(0.06);
318  pt2->SetTextColor(kBlue);
319  pt2->SetFillColor(0);
320  pt2->SetFillStyle(0);
321  pt2->SetBorderSize(0);
322  // pt2->SetTextAlign(12);
323  pt2->AddText("Proton");
324  pt2->Draw();
325 
327  TPaveText * pt3 = NULL;
328  if (d_norm_PercentErr > 0 && vFits[3]->GetMaximum() > 0.5)
329  {
330  pt3 = new TPaveText(0.70, 0.55, 0.94, 0.60, "NB NDC"); //x1, y1, x2, y2
331  pt3->SetTextSize(0.025);
332  pt3->SetFillColor(0);
333  pt3->SetFillStyle(0);
334  pt3->SetBorderSize(0);
335  pt3->AddText(Form("Amplitude error = %5.2f%%", d_norm_PercentErr));
336  pt3->Draw();
337 
338  pt3 = new TPaveText(0.70, 0.60, 0.94, 0.65, "NB NDC"); //x1, y1, x2, y2
339  pt3->SetTextSize(0.025);
340  pt3->SetFillColor(0);
341  pt3->SetFillStyle(0);
342  pt3->SetBorderSize(0);
343  pt3->AddText(Form("Chi-sq/NDF = %5.2f", chiSquaredperNDF));
344  pt3->Draw();
345  }
346 
347 
348 
349  h->Draw("P E1 X0 SAME");
350  c1->Update();
351  string filename = name + h->GetName() + "_" + year + ".png";
352  c1->SaveAs(filename.c_str());
353 
354  delete pt1;
355  delete pt2;
356  delete pt3;
357  delete leg;
358  delete h;
359  delete c1;
360 }
361 
362 
363 
364 
365 
366 
367 
368 
369 //__________________________________________________________________________________________________________________________
370 void Print1DMassSqFits (TH1D * h1, vector<TF1*> vFits, const string name, const string year, bool titleOn = false)
371 {
372  if (h1 == NULL)
373  {
375  cerr<<"[WARNING]: Specified histogram not found!"<<endl;
376  return;
377  }
378 
379  TH1D * h = (TH1D*)h1->Clone();
380 
381  TCanvas * c1 = new TCanvas("c1", "c1", 1800, 1200); //1800, 1200
382  gStyle->SetOptStat(0); // set stat box to show number of entries
383  gPad->SetGridx(1);
384  gPad->SetGridy(1);
385 
386  c1->SetLogy(1);
387  c1->SetTopMargin(0.03);
388  c1->SetBottomMargin(0.1);
389  c1->SetRightMargin(0.05);
390  c1->SetLeftMargin(0.1);
391 
392  if (titleOn == true)
393  {
394  c1->SetTopMargin(0.08);
395  gStyle->SetTitleX(0.5);
396  gStyle->SetTitleAlign(23);
397  gStyle->SetTitleBorderSize(0);
398  }
399 
400  h->GetYaxis()->SetTitle("Counts");
401  h->GetXaxis()->SetTitleOffset(1.3);
402  h->GetYaxis()->SetTitleOffset(1.3);
403 
404  h->SetLineWidth(3);
405  h->SetMinimum(0.1);
406  h->Draw();
407 
408  std::vector<Color_t> colors = {kRed, kGreen, kBlack, kAzure+10, kBlue, kMagenta};
409  for (int k = int(vFits.size())-1; k >= 0; k--)
410  {
411  if (k != 0)
412  vFits[k]->SetLineStyle(9);
413 
414  vFits[k]->SetLineColor(colors[k]);
415  vFits[k]->SetLineWidth(3);
416  vFits[k]->Draw("HIST SAME L X0"); //"E1 C P X0 SAME", "CSAME"
417  }
418 
419  RedrawBorder();
420  c1->Update();
421 
422  string filename = name + h->GetName() + "_" + year + ".png";
423  c1->SaveAs(filename.c_str());
424 
425  delete h;
426  delete c1;
427 }
428 
429 
430 
431 
432 
433 //__________________________________________________________________________________________________________________________
434 void Print1DHist (const TH1D * h1,
435  const string &name,
436  const string &year,
437  int RebinX = 1,
438  bool titleOn = false,
439  bool logY = 1,
440  double zMin = 1,
441  double zMax = 0,
442  double xMin = 0,
443  double xMax = 0,
444  const string &YAxisTitle = "Counts")
445 {
446  if (h1 == NULL)
447  {
449  cerr<<"[WARNING]: Specified histogram "<<name<<" not found!"<<endl;
450  return;
451  }
452 
453  TH1D * h = (TH1D*)h1->Clone();
454 
455  TCanvas * c1 = new TCanvas("c1", "c1", 1200, 1200); //1800, 1200
456  gStyle->SetOptStat(0); // set stat box to show number of entries
457  gPad->SetGridx(1);
458  gPad->SetGridy(1);
459 
460  c1->SetLogy(logY);
461  c1->SetTopMargin(0.04);
462  c1->SetBottomMargin(0.1);
463  c1->SetRightMargin(0.05);
464  c1->SetLeftMargin(0.1);
465 
466  if (titleOn == true)
467  {
468  c1->SetTopMargin(0.08);
469  gStyle->SetTitleX(0.5);
470  gStyle->SetTitleAlign(23);
471  gStyle->SetTitleBorderSize(0);
472  }
473 
474  h->GetYaxis()->SetTitle(YAxisTitle.c_str());
475  h->GetXaxis()->SetTitleOffset(1.3);
476  h->GetYaxis()->SetTitleOffset(1.3);
477 
478  h->Rebin(RebinX);
479 
480  h->SetLineWidth(3);
481  h->SetMinimum(zMin); //1
482  if (zMax > 0) h->SetMaximum(zMax);
483  if (xMax != xMin) h->SetAxisRange(xMin, xMax, "X");
484 
485 
486  h->Draw(); //"HIST CP"
487  RedrawBorder(true);
488  c1->Update();
489 
490  string filename = name + h->GetName() + "_" + year + ".png";
491  c1->SaveAs(filename.c_str());
492 
493  delete h;
494  delete c1;
495 }
496 
497 
498 
499 
500 
501 //__________________________________________________________________________________________________________________________
502 void Print2DHist (const TH2D * h1,
503  const string &name,
504  const string &year,
505  double yMin = 0,
506  double yMax = 0,
507  const string &YAxisTitle = "Counts")
508 {
509  if (h1 == NULL)
510  {
512  cerr<<"[WARNING]: Specified histogram not found!"<<endl;
513  return;
514  }
515 
516  TH1D * h = (TH1D*)h1->Clone();
517 
518  TCanvas * c1 = new TCanvas("c1", "c1", 1800, 1200); //1800, 1200
519  gStyle->SetOptStat(0); // set stat box to show number of entries
520  gPad->SetGridx(1);
521  gPad->SetGridy(1);
522 
523  c1->SetLogz(1);
524  c1->SetTopMargin(0.03);
525  c1->SetBottomMargin(0.10);
526  c1->SetRightMargin(0.15);
527  c1->SetLeftMargin(0.1);
528 
529  h->GetZaxis()->SetTitle(YAxisTitle.c_str());
530  h->GetXaxis()->SetTitleOffset(1.3);
531  h->GetYaxis()->SetTitleOffset(1.3);
532 
533  //h->SetMinimum(1);
534  //h->SetMaximum(zmax);
535  if (yMax != yMin) h->SetAxisRange(yMin, yMax, "Y");
536 
537  h->Draw("colz");
538  c1->Update();
539 
540  string filename = name + h->GetName() + "_" + year + ".png";
541  c1->SaveAs(filename.c_str());
542  delete c1;
543 }
544 
545 
546 
547 
548 
549 //__________________________________________________________________________________________________________________________
550 void Print1DHistVector (const vector<TH1D*> vHist1D,
551  const string &name,
552  const string &year,
553  bool titleOn = false,
554  bool logY = 1,
555  double zMin = 1,
556  double zMax = 0,
557  double xMin = 0,
558  double xMax = 0,
559  const string &YAxisTitle = "Counts")
560 {
561  TCanvas * c1 = new TCanvas("c1", "c1", 1200, 1200); //1800, 1200
562  gStyle->SetOptStat(0); // set stat box to show number of entries
563  gPad->SetGridx(1);
564  gPad->SetGridy(1);
565 
566  c1->SetLogy(logY);
567  c1->SetTopMargin(0.04);
568  c1->SetBottomMargin(0.1);
569  c1->SetRightMargin(0.05);
570  c1->SetLeftMargin(0.1);
571 
572  if (titleOn == true)
573  {
574  c1->SetTopMargin(0.08);
575  gStyle->SetTitleX(0.5);
576  gStyle->SetTitleAlign(23);
577  gStyle->SetTitleBorderSize(0);
578  }
579 
580 
581  TH1D * h = new TH1D();
582  TLegend * leg = new TLegend(0.75, 0.60, 0.94, 0.95); //x1, y1, x2, y2
583 
584  for (int k = 0; k < int(vHist1D.size()); k++)
585  {
586  h = (TH1D*)vHist1D[k]->Clone();
587 
588  h->GetYaxis()->SetTitle(YAxisTitle.c_str());
589  h->GetXaxis()->SetTitleOffset(1.3);
590  h->GetYaxis()->SetTitleOffset(1.3);
591 
592  h->SetLineWidth(3);
593  h->SetLineColor(vParticleColor[k]);
594 
595  h->SetMinimum(zMin); //1
596  if (zMax > 0) h->SetMaximum(zMax);
597  if (xMax != xMin) h->SetAxisRange(xMin, xMax, "X");
598 
599  if (k == 0) h->Draw("HIST CP"); //"HIST CP"
600  else h->Draw("HIST CP same");
601 
602  int p0_percent_change = -10*(k-3);
603  string p0_label;
604  if (p0_percent_change > 0) p0_label = " + " + to_string(p0_percent_change) + "%";
605  else if (p0_percent_change == 0) p0_label = " ";
606  else p0_label = " - " + to_string(-p0_percent_change)+ "%";
607  leg->AddEntry(h, Form("p_{0}%s", p0_label.c_str()), "l");
608  }
609 
610  leg->SetBorderSize(0);
611  leg->Draw();
612 
613  RedrawBorder(true);
614  c1->Update();
615  string filename = name + vHist1D[0]->GetName() + "_" + year + "_overlaid.png";
616  c1->SaveAs(filename.c_str());
617 
618  delete h;
619  delete c1;
620 }
621 
622 
623 
624 
625 //__________________________________________________________________________________________________________________________
626 void make_y_pT_plots(const TH2D * h1,
627  const string& name,
628  const string& year,
629  double zmax=350e3,
630  double zmin=10,
631  bool logZ=1,
632  string ZAxisTitle="Counts",
633  bool printText = false,
634  bool zoomDeuteron = false)
635 {
636  TH2D * h = (TH2D*)h1->Clone();
637 
638  // gROOT->Reset();
639  // TStyle * plain = new TStyle("plain","plain");
640  // plain->SetCanvasBorderMode(0);
641  // plain->SetPadBorderMode(0);
642  // plain->SetPadColor(0);
643  // plain->SetCanvasColor(0);
644  // plain->SetTitleColor(1);
645  // plain->SetStatColor(0);
646  // plain->SetTitleFillColor(0);
647 
648  // gROOT->SetStyle("plain");
649 
650  // int NRGBs = 4;
651  // int NCont = 100;
652  // double stops[] = {0.00, 0.33, 0.67, 1.00};
653 
654  // double red[] = { 1, 1.00, 1.00, 0 };
655  // double green[] = { 1, 1, 0, 0.00 };
656  // double blue[] = { 1, 0, 0.00, 0.00 };
657 
658  // TColor::CreateGradientColorTable(NRGBs, stops, red,
659  // green, blue, NCont);
660  // gStyle->SetNumberContours(NCont);
661 
662 
663  TCanvas *c1 = new TCanvas("c1", "c1", 1200, 1200);
664  gPad->SetGridx(1);
665  gPad->SetGridy(1);
666  gStyle->SetOptStat(0); // set stat box to show number of entries
667 
668  c1->SetLogy(0);
669  c1->SetLogx(0);
670  c1->SetLogz(logZ);
671  c1->SetTopMargin(0.04);
672  c1->SetBottomMargin(0.10);
673  c1->SetRightMargin(0.17);
674  c1->SetLeftMargin(0.13);
675 
676  h->GetZaxis()->SetTitle(ZAxisTitle.c_str());
677  h->GetYaxis()->SetTitle("p_{T} (GeV/c)");
678 
679  h->GetXaxis()->SetTitleOffset(0.95);
680  h->GetYaxis()->SetTitleOffset(1.2);
681 
682 
683  if (zoomDeuteron == false)
684  {
685  h->SetAxisRange(0, 1.9, "Y"); //pT 0, 0.5; 0, 1.9 //special deuteron plot
686  h->SetAxisRange(-2.0, 1.9, "X"); //y -1.5, -0.2; -2.0, 1.9 //special deuteron plot
687  }
688  else
689  {
690  h->SetAxisRange(0, 1.0, "Y"); //pT 0, 0.5; 0, 1.9 //special deuteron plot
691  h->SetAxisRange(-1.7, 0.3, "X"); //y -1.5, -0.2; -2.0, 1.9 //special deuteron plot
692  }
693 
694  h->SetMinimum(zmin);
695  h->SetMaximum(zmax);
696 
697  h->GetXaxis()->SetTitleSize(0.05);
698  h->GetYaxis()->SetTitleSize(0.05);
699  h->GetZaxis()->SetTitleSize(0.05);
700 
701  if (printText == true)
702  {
703  gStyle->SetPaintTextFormat("3.1f"); //5.3f
704  h->SetMarkerSize(1.2);
705  h->Draw("colz TEXT01"); //TEXT89
706  }
707  else
708  {
709  h->Draw("colz"); //colz TEXT
710  }
711 
712 
713  c1->Update();
714 
715  string filename = name + h->GetName() + "_" + year + ".png";
716  c1->SaveAs(filename.c_str());
717  delete c1;
718  // gROOT->Reset();
719 }
720 
721 
722 
723 
724 
725 //__________________________________________________________________________________________________________________________
726 void makeNormalizationPlots(TH1D * hist_all, const TH1D * hist_d, const string &name, const string &year, bool IsPrintRMS = false)
727 {
728  TCanvas * c1 = new TCanvas("c1", "c1", 1800, 1200);
729 
730  double integral_all = hist_all->Integral();
731  double integral_d = hist_d->Integral();
732  double normalization = integral_all/integral_d;
733 
734  TH1D * hist_d_scaled = (TH1D*) hist_d->Clone();
735  hist_d_scaled->Scale(normalization); //scale h_empty_scaled by normalization
736 
737 
738  gStyle->SetOptStat(0); // set stat box to show number of entries
739  gPad->SetGridx(1);
740  gPad->SetGridy(1);
741  c1->SetLogy(1);
742 
743  c1->SetTopMargin(0.03);
744  c1->SetBottomMargin(0.1);
745  c1->SetRightMargin(0.05);
746  c1->SetLeftMargin(0.1);
747 
748  hist_all->SetLineColor(kBlack);
749  hist_all->SetLineWidth(3);
750  hist_all->GetXaxis()->SetTitleOffset(1.3);
751  hist_all->GetYaxis()->SetTitleOffset(1.3);
752  hist_all->Draw();
753 
754  hist_d_scaled->SetLineColor(kRed);
755  hist_d_scaled->SetLineWidth(3);
756  hist_d_scaled->Draw("same");
757 
758  c1->Update();
759  string filename = name + hist_all->GetName() + "_OverlayScaled_" + year + ".png";
760  c1->SaveAs(filename.c_str());
761 
762  if (IsPrintRMS == true)
763  cout<<"RMS hist_all = "<<hist_all->GetStdDev()<<" +/- "<<hist_all->GetStdDevError()<<". RMS hist_d_scaled = "<<hist_d->GetStdDev()<<" +/- "<<hist_d->GetStdDevError()<<endl;
764 
765  delete hist_d_scaled;
766  delete c1;
767 }
768 
769 
770 
771 
772 
773 
774 
775 //__________________________________________________________________________________________________________________________
776 void makeMomentumTOFmassSqPlots (TH2D * h, const string &name, const string &year)
777 {
778  if (h == NULL)
779  {
781  cerr<<"[WARNING]: Specified histogram not found!"<<endl;
782  return;
783  }
784 
785  TCanvas * c1 = new TCanvas("c1", "c1", 1800, 1200);
786 
787  gStyle->SetOptStat(0); // set stat box to show number of entries
788  gPad->SetGridx(1);
789  gPad->SetGridy(1);
790 
791  c1->SetLogz(1);
792  c1->SetTopMargin(0.04);
793  c1->SetBottomMargin(0.09);
794  c1->SetRightMargin(0.12);
795  c1->SetLeftMargin(0.07);
796 
797  h->GetZaxis()->SetTitle("Counts");
798  h->GetXaxis()->SetTitleOffset(1.0);
799  h->GetYaxis()->SetTitleOffset(1.0);
800  h->GetZaxis()->SetTitleOffset(0.9);
801 
802  h->RebinY(2); //2, 4
803 
804  h->SetAxisRange(4, 15, "X");
805  h->SetMinimum(1);
806  //h->SetMaximum(zmax);
807 
808  h->Draw("colz");
809 
810  c1->Update();
811 
812  string filename = name + year + ".png";
813  c1->SaveAs(filename.c_str());
814  delete c1;
815 }
816 
817 
818 
819 //__________________________________________________________________________________________________________________________
820 void PrintPhaseSpace (TH2D * h1, const string &name, const string &year, bool trimPhaseSpace = true, int RebinX = 1, int RebinY = 1, bool printText = true)
821 {
822  TH2D * h = (TH2D*)h1->Clone();
823 
824  TCanvas * c1 = new TCanvas("c1", "c1", 1800, 1200);
825  gStyle->SetOptStat(0); // set stat box to show number of entries
826  gPad->SetGridx(1);
827  gPad->SetGridy(1);
828 
829  c1->SetLogz(1);
830  c1->SetTopMargin(0.04);
831  c1->SetBottomMargin(0.09);
832  c1->SetRightMargin(0.12);
833  c1->SetLeftMargin(0.07);
834 
835  h->GetZaxis()->SetTitle("Counts");
836  h->GetXaxis()->SetTitleOffset(1.0);
837  h->GetYaxis()->SetTitleOffset(1.0);
838  h->GetZaxis()->SetTitleOffset(0.9);
839 
840  h->RebinX(RebinX);
841  h->RebinY(RebinY); // 2 //4 //h->RebinY(8);
842 
843  if (trimPhaseSpace == true)
844  {
845  h->SetAxisRange(4, 15, "X"); //p
846  h->SetAxisRange(0, 1.9, "Y"); //pT
847  }
848 
849  h->SetMinimum(1);
850  //h->SetMaximum(zmax);
851 
852 
853  if (printText == true)
854  {
855  gStyle->SetPaintTextFormat("4.0f");
856  h->SetMarkerSize(1.2);
857  h->Draw("colz TEXT"); //"colz TEXT89"
858  }
859  else
860  {
861  h->Draw("colz"); //colz TEXT
862  }
863 
864  string filename = name + h->GetName() + "_" + year + ".png";
865  c1->Update();
866  c1->SaveAs(filename.c_str());
867 
868  delete h;
869  delete c1;
870 }
871 
872 
873 
874 
875 //__________________________________________________________________________________________________________________________
876 void SimpleCleanUpDEdxTOFmassSqPlot (TH2D * h1)
877 {
878  double entries = 0; //to keep track of histogram entries after bin subtraction
879 
880  for (int i = 0; i <= h1->GetNbinsX(); i++)
881  {
882  double dEdx_value = h1->GetXaxis()->GetBinCenter(i);
883 
884  for (int j = 0; j <= h1->GetNbinsY(); j++)
885  {
886  double m2_value = h1->GetYaxis()->GetBinCenter(j);
887 
888  if (dEdx_value > 1.16 && m2_value > 2.5)
889  {
890  h1->SetBinContent(i, j, 0);
891  h1->SetBinError(i, j, 0);
892  }
893 
894  entries = entries + h1->GetBinContent(i, j);
895  }
896  }
897 
898  h1->SetEntries(entries);
899 }
900 
901 
902 
903 
904 //__________________________________________________________________________________________________________________________
905 void BetheBlochCleanUpDEdxTOFmassSqPlot (TH2D * h1, double pBinCenter)
906 {
907  double pion_dEdx_cut = BetheBlochAntoniMod(pBinCenter/0.335); //see output plot of MakeBetheBlochPlot
908  double entries = 0; //to keep track of histogram entries after bin subtraction
909 
910  for (int i = 0; i <= h1->GetNbinsX(); i++)
911  {
912  double dEdx_value = h1->GetXaxis()->GetBinCenter(i);
913 
914  for (int j = 0; j <= h1->GetNbinsY(); j++)
915  {
916  double m2_value = h1->GetYaxis()->GetBinCenter(j);
917 
918  if (dEdx_value > pion_dEdx_cut && m2_value > 2.5)
919  {
920  h1->SetBinContent(i, j, 0);
921  h1->SetBinError(i, j, 0);
922  }
923 
924  entries = entries + h1->GetBinContent(i, j);
925  }
926  }
927 
928  h1->SetEntries(entries);
929 }
930 
931 
932 
933 
934 
935 //__________________________________________________________________________________________________________________________
936 void makedEdxTOFmassSqPlots (TH2D * h1, double pBinCenter, int RebinY, const string &name, string year, bool titleOn = false, bool cleanUpPions = false, string drawBox = "", int nD = 0)
937 {
938  if (h1 == NULL)
939  {
941  cerr<<"[WARNING]: Specified histogram not found!"<<endl;
942  return;
943  }
944 
945  TH2D * h = (TH2D*)h1->Clone();
946 
947  TCanvas * c1 = new TCanvas("c1", "c1", 1800, 1200);
948  gStyle->SetOptStat(0);
949  gPad->SetGridx(1);
950  gPad->SetGridy(1);
951 
952  c1->SetLogz(1);
953  c1->SetBottomMargin(0.09);
954  c1->SetRightMargin(0.12);
955  c1->SetLeftMargin(0.07);
956 
957  if (titleOn == false)
958  {
959  c1->SetTopMargin(0.04);
960  h->SetTitle("");
961  }
962  else c1->SetTopMargin(0.1);
963 
964  h->GetZaxis()->SetTitle("Counts");
965  h->GetXaxis()->SetTitleOffset(1.0);
966  h->GetYaxis()->SetTitleOffset(1.0);
967  h->GetZaxis()->SetTitleOffset(0.9);
968 
969 
970  // h->RebinX(2);
971  h->RebinY(RebinY); // 2 //4 //h->RebinY(8);
972 
973 
974  if (pBinCenter < 0) cleanUpPions = false;
975  if (cleanUpPions == true)
976  {
977  //SimpleCleanUpDEdxTOFmassSqPlot(h);
978  BetheBlochCleanUpDEdxTOFmassSqPlot(h, pBinCenter);
979  year = year + "_PionRemoved";
980  }
981 
982  h->SetAxisRange(0.8, 1.78, "X");
983  h->SetMinimum(1);
984  //h->SetMaximum(zmax);
985 
986  h->Draw("colz");
987 
988 
989  if (drawBox == "deuteron")
990  {
991  TBox * DeuteronBox = new TBox(0.9, 3.16, 1.2, 4.2); //x1, y1, x2, y2
992  DeuteronBox->SetFillColorAlpha(0, 0.);
993  DeuteronBox->SetFillStyle(0);
994  DeuteronBox->SetLineColor(kRed);
995  DeuteronBox->SetLineWidth(4);
996  DeuteronBox->Draw();
997 
998  if (nD > 0)
999  {
1002  gStyle->SetOptStat("e");
1003  TPaveStats * st = (TPaveStats*)c1->GetPrimitive("stats");
1004  st->SetX1NDC(0.60); //new x start position
1005  st->SetX2NDC(0.85); //new x end position
1006  st->SetY1NDC(0.85); //new y start position
1007  st->SetY2NDC(0.90); //new y end position
1008 
1009 
1011  TPaveText * pt = new TPaveText(0.60, 0.80, 0.85, 0.84, "NB NDC"); //x1, y1, x2, y2
1012  pt->SetTextSize(0.04);
1013  pt->SetFillColor(0);
1014  // pt->SetTextAlign(12);
1015  pt->AddText(Form("nDeuterons = %d", nD));
1016  pt->Draw();
1017  }
1018  }
1019 
1020  if (drawBox == "piMinus")
1021  {
1022  TBox * PiMinusBox = new TBox(1.22, -0.5, 1.40, 0.5); //x1, y1, x2, y2
1023  PiMinusBox->SetFillColorAlpha(0, 0.);
1024  PiMinusBox->SetFillStyle(0);
1025  PiMinusBox->SetLineColor(kBlack);
1026  PiMinusBox->SetLineWidth(4);
1027  PiMinusBox->Draw();
1028  }
1029 
1030 
1031  string filename = name + year + ".png";
1032  c1->Update();
1033  c1->SaveAs(filename.c_str());
1034  delete c1;
1035 }
1036 
1037 
1038 
1039 
1040 //__________________________________________________________________________________________________________________________
1041 void Make_YPt_ParticleSpectra (const TH2D * h,
1042  const string &particle_name,
1043  const string &MC_type,
1044  double year,
1045  double yAxis_max,
1046  double yAxis_min,
1047  int rapidity_bin_max_data,
1048  int rapidity_bin_min_data,
1049  bool makeFit = false,
1050  bool isLogYAxis = true,
1051  bool doScaling = true,
1052  bool showEPOS = false)
1053 {
1054  double pT_max = 1.6;
1055  if (particle_name == "deuteron" && makeFit == 0) pT_max = 1.0;
1056  else if (particle_name == "deuteron" && makeFit == 1) pT_max = 1.0;
1057  else if (particle_name == "proton" && makeFit == 0) pT_max = 1.6;
1058  else if (particle_name == "proton" && makeFit == 1) pT_max = 1.8;
1059 
1060  TMultiGraph * MG1 = new TMultiGraph();
1061  MG1->SetTitle("");
1062 
1063 
1064  /*
1065  vector<TGraph*> SzymonGraphs = GetDataPointsSzymon(particle_name);
1066  for (int k = 0; k < int(SzymonGraphs.size()); k++)
1067  {
1068  //MG1->Add(SzymonGraphs[k], "cp"); //comment this out for no Szymon comparison
1069  }
1070  */
1071 
1072 
1073  //_________________________________________________________
1074  vector<TGraphErrors*> vDataGraphs = GetGraphsfromTH2D(h, year, particle_name, rapidity_bin_max_data, rapidity_bin_min_data, doScaling);
1075  for (int k = 0; k < int(vDataGraphs.size()); k++)
1076  {
1077  vDataGraphs[k] = RemoveZeroFromGraph(vDataGraphs[k]);
1078  vDataGraphs[k]->SetLineWidth(5);
1079  vDataGraphs[k]->SetLineColor(vParticleColor[k]);
1080  vDataGraphs[k]->SetMarkerSize(2);
1081  vDataGraphs[k]->SetMarkerStyle(20); //24, 20
1082  vDataGraphs[k]->SetMarkerColor(vParticleColor[k]);
1083 
1084  if (makeFit == true)
1085  MG1->Add(vDataGraphs[k], "p");
1086  else
1087  MG1->Add(vDataGraphs[k], "cp");
1088  }
1089 
1090 
1091  //_________________________________________________________
1092  TCanvas * c2 = new TCanvas("c2", "c2", 1600, 1600);
1093  c2->SetLogy(isLogYAxis);
1094  c2->SetTicky(1);
1095  c2->SetTickx(1);
1096  c2->SetLeftMargin(0.19);
1097  c2->SetRightMargin(0.025);
1098  c2->SetBottomMargin(0.12);
1099  c2->SetTopMargin(0.02); if (!isLogYAxis) c2->SetTopMargin(0.04);
1100  MG1->Draw("ap");
1101 
1102 
1103  //_________________________________________________________
1104  vector<TF1*> vDataFits;
1105  if (makeFit == true)
1106  {
1107  vector<TGraphErrors*> vGFitConfidenceBands; //to store graphs with fits and their confidence interval bands.
1108  vector<int> vMissingPointIndex; // to store the index of any missing points that are added, for later use.
1109  vDataFits = GetPtSpectraFitAndEval(vDataGraphs, vGFitConfidenceBands, vMissingPointIndex, particle_name);
1110  // cerr<<"Missing point index: "; printVector(vMissingPointIndex); cerr<<endl; /// FIXME, but not high priority.
1111 
1112  for (int k = 0; k < int(vDataFits.size()); k++)
1113  {
1114  vDataFits[k]->SetLineColor(vParticleColor[k]);
1115  vDataFits[k]->SetLineWidth(3);
1116  vDataFits[k]->Draw("same");
1117 
1118  if (particle_name == "deuteron")
1119  {
1120  vGFitConfidenceBands[k]->SetFillColorAlpha(vParticleColor[k]-2, 0.30);
1121  vGFitConfidenceBands[k]->SetFillStyle(1001);
1122  vGFitConfidenceBands[k]->Draw("e3 same");
1123  }
1124  }
1125  }
1126 
1127 
1128  //_________________________________________________________
1129  if (showEPOS == true)
1130  {
1131  cerr<<"\n[INFO] Adding EPOS comparisons..."<<endl;
1132  const vector<double> vP0 = {109.45, 101.03, 92.61, 84.19, 75.77, 67.35, 58.94};
1133 
1135  double pT_max_EPOS = 1.0;
1136  int j_min = 3; int j_max = 4;
1137  int y_bin_min_EPOS = 48; int y_bin_max_EPOS = 51;
1138  if (particle_name == "deuteron")
1139  {
1140  pT_max_EPOS = 0.6;
1141  j_min = 0; j_max = int(vP0.size());
1142  y_bin_min_EPOS = 45; y_bin_max_EPOS = 48;
1143  }
1144 
1146  for (int j = j_min; j < j_max; j++)
1147  {
1148  vector<TGraphErrors*> vG_EPOS = GetEPOS_Y_PT_Spectra_Graphs(particle_name, vP0[j], 2, doScaling, pT_max_EPOS, y_bin_max_EPOS, y_bin_min_EPOS);
1149  if (vDataGraphs.size() != vG_EPOS.size())
1150  cerr<<"\n\n[WARNING] DataGraphs size != vG_EPOS size"<<endl;
1151 
1152 
1154  if (particle_name == "proton")
1155  {
1156  vector<TGraphErrors*> vGFitConfidenceBands; //to store graphs with fits and their confidence interval bands.
1157  vector<TF1*> vEPOSTestFits = GetPtSpectraFit(vG_EPOS, vGFitConfidenceBands, particle_name, true, 0.139, 0.1, 1.5, false);
1158  for (int k = 0; k < int(vG_EPOS.size()); k++)
1159  {
1160  vEPOSTestFits[k]->SetLineColor(vParticleColor[k]);
1161  vEPOSTestFits[k]->SetLineWidth(1);
1162  vEPOSTestFits[k]->Draw("c same");
1163  }
1164  }
1165 
1166 
1167  for (int k = 0; k < int(vG_EPOS.size()); k++)
1168  {
1169  vG_EPOS[k]->SetLineColor(vParticleColor[k]);
1170  vG_EPOS[k]->SetLineWidth(1);
1171  vG_EPOS[k]->SetLineStyle(2);
1172  vG_EPOS[k]->Draw("c same");
1173  }
1174  }
1175 
1177  if (particle_name == "deuteron")
1178  {
1179  vector<TGraphAsymmErrors*> vG_EPOSBands = GetEPOS_Y_PT_SpectraBand_Graphs(particle_name);
1180  if (vG_EPOSBands.size() != vDataGraphs.size()) cerr<<"\n\n[WARNING] vG_EPOSBands size != vDataGraphs size"<<endl;
1181  for (int k = 0; k < int(vG_EPOSBands.size()); k++)
1182  {
1183  if (!isLogYAxis)
1184  vG_EPOSBands[k]->SetFillColorAlpha(vParticleColor[k]-2, 0.3); //(vParticleColor[k]-2, 0.3)
1185  else
1186  vG_EPOSBands[k]->SetFillColorAlpha(vParticleColor[k]-2, 0.3); //(vParticleColor[k]-2, 0.3)
1187 
1188  vG_EPOSBands[k]->SetFillStyle(3008);
1189  vG_EPOSBands[k]->Draw("e3 same");
1190  }
1191  }
1192 
1193  cerr<<"[INFO] Completed EPOS comparisons.\n\n"<<endl;
1194  }
1195 
1196 
1197  //_________________________________________________________
1198  TLegend * leg = GetTLegend(particle_name, year, "ExtraPoints");
1199  for (int k = 0; k < int(vDataGraphs.size()); k++)
1200  {
1201  if (makeFit == true)
1202  {
1203  if (doScaling)
1204  leg->AddEntry(vDataGraphs[k], Form("#times10^{%d} y = %0.1f, T=%4.1f MeV", -1*k, -1.9+2*double(k+rapidity_bin_min_data-1)/10, 1000*vDataFits[k]->GetParameter(1)), "l");
1205  else
1206  leg->AddEntry(vDataGraphs[k], Form("y = %0.1f, T=%4.1f MeV", -1.9+2*double(k+rapidity_bin_min_data-1)/10, 1000*vDataFits[k]->GetParameter(1)), "l");
1207  }
1208 
1209  else
1210  {
1211  if (doScaling)
1212  leg->AddEntry(vDataGraphs[k], Form("#times10^{%d} y = %0.1f", -1*k, -1.9+2*double(k+rapidity_bin_min_data-1)/10), "pl");
1213  else
1214  leg->AddEntry(vDataGraphs[k], Form("y = %0.1f", -1.9+2*double(k+rapidity_bin_min_data-1)/10), "pl");
1215  }
1216  }
1217  leg->SetBorderSize(0);
1218  leg->Draw();
1219 
1220 
1221  //_________________________________________________________
1222  MG1->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1223  MG1->GetXaxis()->SetTitleOffset(1.0);
1224  MG1->GetXaxis()->SetTitleSize(0.05);
1225  MG1->GetYaxis()->SetTitle("#frac{d^{2}n}{dydp_{T}} (GeV/c)^{-1}");
1226  MG1->GetYaxis()->SetTitleOffset(1.65);
1227  MG1->GetYaxis()->SetTitleSize(0.05);
1228  MG1->GetHistogram()->SetMaximum(yAxis_max); // along
1229  MG1->GetHistogram()->SetMinimum(yAxis_min); // Y
1230  MG1->GetXaxis()->SetLimits(0.0, pT_max);
1231  c2->Update();
1232 
1233  string name = "plots_spectra/" + particle_name + "/" + "spectra_" + to_string(int(year)) + "_" + particle_name + "_" + MC_type + ".png";
1234  c2->SaveAs(name.c_str());
1235 
1236  bool makePDFtoPNG = false;
1237  if (makePDFtoPNG == true && makeFit == true && showEPOS == true)
1238  {
1239  name = "plots_spectra/" + particle_name + "/" + "spectra_" + to_string(int(year)) + "_" + particle_name + "_" + MC_type + ".pdf";
1240  c2->SaveAs(name.c_str());
1241 
1242  std::string command = std::string("pdftoppm -png -rx 300 -ry 300 -cropbox ") + name + " " + name;
1243  int conversion_exit_code = system(command.c_str());
1244  cerr<<"'"<<command<<"' returned with code "<<conversion_exit_code<<endl;
1245  }
1246 
1247  delete c2;
1248 }
1249 
1250 
1251 
1252 
1253 //__________________________________________________________________________________________________________________________
1254 void Make_YPt_1DParticleSpectraStacked (const TH2D * h,
1255  const string &particle_name,
1256  const string &MC_type,
1257  double year,
1258  double yAxis_max,
1259  double yAxis_min,
1260  int rapidity_bin_max_data,
1261  int rapidity_bin_min_data,
1262  bool makeFit = false,
1263  bool isLogYAxis = true,
1264  bool doScaling = true)
1265 {
1266  double pT_max = 1.6;
1267  if (particle_name == "deuteron") pT_max = 1.0;
1268 
1269  TMultiGraph * MG1 = new TMultiGraph();
1270  MG1->SetTitle("");
1271 
1272 
1273  //_________________________________________________________
1274  vector<TGraphErrors*> vDataGraphs = GetGraphsfromTH2D(h, year, particle_name, rapidity_bin_max_data, rapidity_bin_min_data, doScaling);
1275  for (int k = 0; k < int(vDataGraphs.size()); k++)
1276  {
1277  vDataGraphs[k] = RemoveZeroFromGraph(vDataGraphs[k]);
1278  vDataGraphs[k]->SetLineWidth(5);
1279  vDataGraphs[k]->SetLineColor(vParticleColor[k]);
1280  vDataGraphs[k]->SetMarkerSize(2);
1281  vDataGraphs[k]->SetMarkerStyle(20); //24, 20
1282  vDataGraphs[k]->SetMarkerColor(vParticleColor[k]);
1283  }
1284 
1285 
1286  //_________________________________________________________
1287  vector<TGraphErrors*> vGFitConfidenceBands; //to store graphs with fits and their confidence interval bands.
1288  vector<int> vMissingPointIndex; // to store the index of any missing points that are added, for later use.
1289  vector<TF1*> vDataFits;
1290  if (makeFit == true)
1291  {
1292  vDataFits = GetPtSpectraFitAndEval(vDataGraphs, vGFitConfidenceBands, vMissingPointIndex, particle_name);
1293 
1294  if (vDataFits.size() != vGFitConfidenceBands.size())
1295  cerr<<"\n\n[ERROR] Size of vGFitConfidenceBands = "<<vGFitConfidenceBands.size()<<" is not equal to size of vDataFits = "<<vDataFits.size()<<endl;
1296  }
1297 
1298 
1299  //_________________________________________________________
1300  TCanvas * c2 = new TCanvas("c2", "c2", 3000, 3000);
1301  int nPlots = int(vDataGraphs.size());
1302 
1303  for (int k = 0; k < nPlots; k++)
1304  {
1305  double fraction = 0.22;
1306  double y_up = 1-(k+0)*fraction;
1307  double y_low = 1-(k+1)*fraction; if (k == nPlots-1) y_low = 0;
1308  TPad * pad = new TPad(Form("pad_%d", k), "", 0.0, y_low, 1.0, y_up); //xlow, ylow, xup, yup,
1309  pad->Draw();
1310  pad->cd();
1311 
1312  pad->SetLogy(isLogYAxis);
1313  pad->SetTickx(1);
1314  pad->SetTicky(1);
1315  pad->SetGridx(1);
1316  pad->SetGridy(1);
1317  pad->SetLeftMargin(0.15);
1318  pad->SetRightMargin(0.02);
1319  pad->SetBottomMargin(0.04); if (k == nPlots-1) pad->SetBottomMargin(0.35);
1320  pad->SetTopMargin(0.10);
1321 
1322  vDataGraphs[k]->Draw("ap");
1323  vDataGraphs[k]->GetXaxis()->SetTitle(""); if (k == nPlots-1) vDataGraphs[k]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1324  vDataGraphs[k]->GetXaxis()->SetTitleOffset(0); if (k == nPlots-1) vDataGraphs[k]->GetXaxis()->SetTitleOffset(0.95);
1325  vDataGraphs[k]->GetXaxis()->SetTitleSize(0); if (k == nPlots-1) vDataGraphs[k]->GetXaxis()->SetTitleSize(0.16);
1326  vDataGraphs[k]->GetXaxis()->SetLabelSize(0); if (k == nPlots-1) vDataGraphs[k]->GetXaxis()->SetLabelSize(0.10);
1327 
1328  vDataGraphs[k]->GetYaxis()->SetTitle("#frac{d^{2}n}{dydp_{T}} (GeV/c)^{-1}");
1329  vDataGraphs[k]->GetYaxis()->SetLabelFont(63);
1330  vDataGraphs[k]->GetYaxis()->SetTitleFont(63);
1331  vDataGraphs[k]->GetYaxis()->SetTitleOffset(2.2); //if (k == 1) vDataGraphs[k]->GetYaxis()->SetTitleOffset(0.40);
1332  vDataGraphs[k]->GetYaxis()->SetTitleSize(0); if (k == 1) vDataGraphs[k]->GetYaxis()->SetTitleSize(90);
1333  vDataGraphs[k]->GetYaxis()->SetLabelSize(60); //if (k == 1) vDataGraphs[k]->GetYaxis()->SetLabelSize(1);
1334 
1335  vDataGraphs[k]->GetHistogram()->SetMaximum(yAxis_max); //if (k == nPlots-1) vDataGraphs[k]->GetHistogram()->SetMaximum(0.5*yAxis_max);
1336  vDataGraphs[k]->GetHistogram()->SetMinimum(yAxis_min); //if (k == nPlots-1) vDataGraphs[k]->GetHistogram()->SetMaximum(yAxis_min); // Y
1337  vDataGraphs[k]->GetXaxis()->SetLimits(0.0, pT_max);
1338 
1339 
1340  //______________________
1341  if (makeFit == true)
1342  {
1343  vDataFits[k]->SetLineColor(vParticleColor[k]);
1344  vDataFits[k]->SetLineWidth(3);
1345  // vDataFits[k]->GetXaxis()->SetLimits(0.0, pT_max);
1346  vDataFits[k]->Draw("same");
1347 
1348  vGFitConfidenceBands[k]->SetFillColorAlpha(vParticleColor[k]-2, 0.30);
1349  vGFitConfidenceBands[k]->SetFillStyle(1001); //3004
1350  // vGFitConfidenceBands[k]->GetXaxis()->SetLimits(0.0, pT_max);
1351  vGFitConfidenceBands[k]->Draw("e3 psame");
1352  }
1353 
1354 
1355  //______________________
1356  TLegend * leg = GetTLegend(particle_name, year, "stackedspectra");
1357  leg->AddEntry(vDataGraphs[k], Form("y = %0.1f", -1.9+2*double(k+rapidity_bin_min_data-1)/10), "l");
1358  leg->SetBorderSize(0);
1359  leg->SetFillStyle(0);
1360  leg->Draw();
1361 
1362  c2->cd();
1363  }
1364 
1365 
1366  //_________________________________________________________
1367  c2->Update();
1368  string name = "plots_spectra/" + particle_name + "/" + "spectra_" + to_string(int(year)) + "_" + particle_name + "_" + MC_type + "_stacked.png";
1369  c2->SaveAs(name.c_str());
1370  delete c2;
1371 }
1372 
1373 
1374 
1375 
1376 //__________________________________________________________________________________________________________________________
1377 void Make_PtIntegrated_rapidity_spectra (const TH2D * h,
1378  const string &particle_name,
1379  const string &MC_type,
1380  double year,
1381  double yAxis_max,
1382  double yAxis_min,
1383  int rapidity_bin_max_data,
1384  int rapidity_bin_min_data,
1385  bool isLogYAxis = true)
1386 {
1391  cerr<<"\n\n[INFO] Making pT-integrated dn/dy spectra plots for "<<particle_name<<" "<<MC_type<<endl;
1392 
1393  //_________________________________________________________
1394  TCanvas * c2 = new TCanvas("c2", "c2", 2000, 2000);
1395  c2->SetLogy(isLogYAxis);
1396  c2->SetTicky(1);
1397  c2->SetTickx(1);
1398  c2->SetLeftMargin(0.18);
1399  c2->SetRightMargin(0.035);
1400  c2->SetBottomMargin(0.12);
1401  c2->SetTopMargin(0.02); if (!isLogYAxis) c2->SetTopMargin(0.04);
1402  TMultiGraph * MG1 = new TMultiGraph();
1403  MG1->SetTitle("");
1404 
1405 
1406  //_________________________________________________________
1407  vector<TGraphErrors*> vDataGraphs = GetGraphsfromTH2D(h, year, particle_name, rapidity_bin_max_data, rapidity_bin_min_data, false);
1408 
1409  //clean vDataGraphs, and save vDataGraphs to file
1410  TFile * fOut = new TFile (Form("plots/Spectra_TGraphErrors_%s.root", MC_type.c_str()),"RECREATE");
1411  for (int i = 0; i < int(vDataGraphs.size()); i++)
1412  {
1413  vDataGraphs[i] = RemoveZeroFromGraph(vDataGraphs[i]);
1414  vDataGraphs[i]->Write();
1415  }
1416  fOut->Close();
1417  delete fOut;
1418 
1419 
1420  vector<TGraphErrors*> vGFitConfidenceBands; //to store confidence interval bands of the fits.
1421  vector<TGraphErrors*> vIntegratedYSpectraGraphs = GetFitEvaluatedIntegratedYSpectraGraphs(vDataGraphs, vGFitConfidenceBands);
1422  for (int k = 0; k < int(vIntegratedYSpectraGraphs.size()); k++)
1423  {
1424  vIntegratedYSpectraGraphs[k]->SetLineColor(vParticleColor[k]);
1425  vIntegratedYSpectraGraphs[k]->SetLineWidth(3);
1426 
1427  MG1->Add(vIntegratedYSpectraGraphs[k], "cp");
1428  // vIntegratedYSpectraGraphs[k]->Draw("same");
1429  }
1430 
1431  MG1->Draw("ap");
1432 
1433  auto GFitConfidenceBand = MergeSinglePointTGraphErrors(vGFitConfidenceBands);
1434  GFitConfidenceBand->SetFillColorAlpha(14, 0.30);
1435  GFitConfidenceBand->SetFillStyle(1001); //3004
1436  GFitConfidenceBand->Draw("e3 same");
1437 
1438 
1439  //_________________________________________________________
1440  const vector<double> vP0 = {58.94, 67.35, 75.77, 84.19, 92.61, 101.03, 109.45};
1441  // const vector<double> vP0 = {84.19, 92.61, 101.03, 109.45};
1442  vector<TH1D*> vHistEPOS = GetEPOSYGraphs(particle_name, vP0);
1443  for (int k = 0; k < int(vHistEPOS.size()); k++)
1444  {
1445  vHistEPOS[k]->SetLineColor(kGreen-10+2*k);
1446  vHistEPOS[k]->SetLineWidth(3);
1447  vHistEPOS[k]->SetLineStyle(9);
1448  vHistEPOS[k]->SetAxisRange(-1.4, -0.2, "X");
1449  vHistEPOS[k]->Draw("HIST c same");
1450 
1451  TLegend * leg_p0 = GetTLegend(particle_name, k, "p0");
1452  leg_p0->SetFillStyle(0);
1453  leg_p0->SetBorderSize(0);
1454  leg_p0->AddEntry("", Form("%4.1f", vP0[k]), ""); //"p" "l";
1455  leg_p0->Draw();
1456  }
1457 
1458 
1459  //_________________________________________________________
1460  TLegend * leg = GetTLegend(particle_name, year, "p0");
1461  leg->SetFillStyle(0);
1462  leg->SetBorderSize(0);
1463  leg->AddEntry("","p_{0} (MeV/c)","");
1464  leg->Draw();
1465 
1466  /*
1467  leg = GetTLegend(particle_name, year, "header");
1468  leg->SetNColumns(2);
1469  leg->SetColumnSeparation(.01);
1470  // leg->SetBorderSize(2);
1471  leg->AddEntry("","High-stat.","");
1472  leg->Draw();
1473  */
1474 
1475 
1476  //_________________________________________________________
1477  MG1->GetXaxis()->SetLimits(-1.4, 0.2);
1478  MG1->GetHistogram()->SetMaximum(yAxis_max); // along
1479  MG1->GetHistogram()->SetMinimum(yAxis_min); // Y
1480  MG1->GetXaxis()->SetTitle("y_{c.m.}");
1481  MG1->GetXaxis()->SetTitleOffset(0.7);
1482  MG1->GetXaxis()->SetTitleSize(0.07);
1483  MG1->GetYaxis()->SetTitle("#frac{dn}{dy}");
1484  MG1->GetYaxis()->SetTitleOffset(1.3);
1485  MG1->GetYaxis()->SetTitleSize(0.06);
1486 
1487 
1488  c2->Update();
1489  RedrawBorder(true);
1490  string name;
1491  name = "plots/hist_spectra_yCM_integrated/spectra_y_" + to_string(int(year)) + "_" + particle_name + "_" + MC_type + ".png";
1492  c2->SaveAs(name.c_str());
1493  name = "plots/hist_spectra_yCM_integrated/spectra_y_" + to_string(int(year)) + "_" + particle_name + "_" + MC_type + ".pdf";
1494  c2->SaveAs(name.c_str());
1495  delete c2;
1496 }
1497 
1498 
1499 
1500 
1501 //__________________________________________________________________________________________________________________________
1502 void Make_YPt_EPOSSpectra (vector<TGraphErrors*> vGraphs,
1503  vector<TF1*> vFits,
1504  const string &name,
1505  const string &particle_name,
1506  const string &MC_type,
1507  double year,
1508  double yAxis_max,
1509  double yAxis_min,
1510  int rapidity_bin_min_data = 45, //41
1511  bool isLogYAxis = true,
1512  bool doScaling = true
1513  )
1514 {
1515  if (vGraphs.size() != vFits.size()) cerr<<"\n\n[WARNING] DataGraphs size != vFits size"<<endl;
1516 
1517  double pT_max = 3.5;
1518 
1519  TMultiGraph * MG1 = new TMultiGraph();
1520  MG1->SetTitle("");
1521 
1522 
1523  //_________________________________________________________
1524  for (int k = 0; k < int(vGraphs.size()); k++)
1525  {
1526  vGraphs[k] = RemoveZeroFromGraph(vGraphs[k]);
1527  vGraphs[k]->SetLineWidth(5);
1528  vGraphs[k]->SetLineColor(vParticleColor[k]);
1529  vGraphs[k]->SetMarkerSize(2.5);
1530  vGraphs[k]->SetMarkerStyle(20); //24, 20
1531  vGraphs[k]->SetMarkerColor(vParticleColor[k]);
1532 
1533  MG1->Add(vGraphs[k], "p");
1534  }
1535 
1536 
1537  //_________________________________________________________
1538  TCanvas * c2 = new TCanvas("c2", "c2", 1200, 1200);
1539  c2->SetLogy(isLogYAxis);
1540  c2->SetTicky(1);
1541  c2->SetTickx(1);
1542  c2->SetLeftMargin(0.19);
1543  c2->SetRightMargin(0.025);
1544  c2->SetBottomMargin(0.12);
1545  c2->SetTopMargin(0.02); if (!isLogYAxis) c2->SetTopMargin(0.04);
1546  MG1->Draw("ap");
1547 
1548 
1549  //_________________________________________________________
1550  for (int k = 0; k < int(vFits.size()); k++)
1551  {
1552  vFits[k]->SetLineColor(vParticleColor[k]);
1553  vFits[k]->SetLineWidth(3);
1554  vFits[k]->Draw("c same");
1555  }
1556 
1557 
1558  //_________________________________________________________
1559  TLegend * leg = GetTLegend(particle_name, year, "ExtraPoints");
1560  for (int k = 0; k < int(vGraphs.size()); k++)
1561  {
1562  if (doScaling)
1563  leg->AddEntry(vGraphs[k], Form("#times10^{%d} y = %0.1f, T=%4.1f MeV", -1*k, -9.9+2*double(k+rapidity_bin_min_data-1)/10, 1000*vFits[k]->GetParameter(1)), "l");
1564  else
1565  leg->AddEntry(vGraphs[k], Form("y = %0.1f, T=%4.1f MeV", -9.9+2*double(k+rapidity_bin_min_data-1)/10, 1000*vFits[k]->GetParameter(1)), "l");
1566  }
1567  leg->SetBorderSize(0);
1568  leg->Draw();
1569 
1570 
1571  //_________________________________________________________
1572  MG1->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1573  MG1->GetXaxis()->SetTitleOffset(1.0);
1574  MG1->GetXaxis()->SetTitleSize(0.05);
1575  MG1->GetYaxis()->SetTitle("#frac{d^{2}n}{dydp_{T}} (GeV/c)^{-1}");
1576  MG1->GetYaxis()->SetTitleOffset(1.65);
1577  MG1->GetYaxis()->SetTitleSize(0.05);
1578  MG1->GetHistogram()->SetMaximum(yAxis_max); // along
1579  MG1->GetHistogram()->SetMinimum(yAxis_min); // Y
1580  MG1->GetXaxis()->SetLimits(0.0, pT_max);
1581  c2->Update();
1582  string filename = name + particle_name + "_" + MC_type + ".png";
1583  c2->SaveAs(filename.c_str());
1584 
1585  delete c2;
1586 }
Generated by   doxygen 1.8.5