simple-tof-analysis
 All Classes Namespaces Functions Variables Groups Pages
Utilities.h
1 using namespace std;
2 
3 #include "TVirtualFitter.h"
4 
5 
6 //__________________________________________________________________________________________________________________________
7 vector<double> GetErrorN(double n, bool debug = false)
8 {
14  if (n == 0)
15  {
16  vector<double> result;
17  result.push_back(0); //error_low
18  result.push_back(0); //error_high
19 
20  if (debug == true)
21  cerr<<"| Setting errors to zero."<<endl;
22  return result;
23  }
24  else if (n < 10000)
25  {
26  vector<double> result;
27  result.push_back(-0.5 + sqrt(n+0.25)); //error_low
28  result.push_back(0.5 + sqrt(n+0.25)); //error_high
29 
30  return result;
31  }
32  else
33  {
34  vector<double> result;
35  result.push_back(sqrt(n)); //error_low
36  result.push_back(sqrt(n)); //error_high
37 
38  return result;
39  }
40 }
41 
42 
43 
44 
45 //__________________________________________________________________________________________________________________________
46 void SetEmptyBinsToOne (TH1D * h)
47 {
48  for (int j = 0; j <= h->GetNbinsX(); j++)
49  {
50  double bin_content = h->GetBinContent(j);
51  if (bin_content < 1)
52  bin_content = 1;
53 
54  h->SetBinContent(j, bin_content);
55  }
56 }
57 
58 
59 
60 
61 //__________________________________________________________________________________________________________________________
62 void NearestNeighborSmoothing (TH1D * h, double xmin, double xmax, int nPassMax = 10)
63 {
64  int iPass = 0;
65  while(iPass < nPassMax)
66  {
67  for (int j = 0; j <= h->GetNbinsX(); j++)
68  {
69  if (h->GetBinCenter(j) >= xmin && h->GetBinCenter(j) < xmax)
70  {
71  double bin_content_xlow = h->GetBinContent(j-1);
72  double bin_content_xup = h->GetBinContent(j+1);
73  double bin_content = (bin_content_xlow + bin_content_xup)/2;
74  h->SetBinContent(j, bin_content);
75  }
76  }
77  ++iPass;
78  }
79 }
80 
81 
82 
83 
84 
85 //__________________________________________________________________________________________________________________________
86 TH1D * GetScaledTH1D(double scale_factor, const TH1D * h1)
87 {
88  TGraph * g = new TGraph(h1);
89  // g->SetBit(TGraph::kIsSortedX);
90 
91  TH1D * h = (TH1D*)h1->Clone();
92  h->Reset("M");
93 
94  for (int j = 1; j <= h->GetNbinsX(); j++)
95  {
96  double bin_content = scale_factor*g->Eval(h->GetBinCenter(j), 0, "S");
97  // double bin_error = GetErrorN(bin_content)[0];
98  h->SetBinContent(j, bin_content);
99  h->SetBinError(j, 0); //h->SetBinError(j, bin_error);
100  }
101 
102  return h;
103 }
104 
105 
106 
107 
108 //__________________________________________________________________________________________________________________________
109 TH1D * Normalize_1D_Histogram (TH1D* h1, double nEvts, bool diagnostics = false)
110 {
115  // h1->Sumw2();
116 
117  for (int i = 0; i <= h1->GetNbinsX(); i++)
118  {
119  if (h1->GetBinContent(i) > 0)
120  {
121  double scaleFactor = 1/(nEvts*h1->GetXaxis()->GetBinWidth(i));
122  double bin_content = h1->GetBinContent(i);
123  double bin_error = h1->GetBinError(i);
124  h1->SetBinContent( i, scaleFactor*bin_content);
125  h1->SetBinError( i, scaleFactor*bin_error);
126 
127  if (diagnostics == true)
128  {
129  cout<<"Bin "<<i<<": "<<h1->GetBinContent(i)<<", "<<h1->GetBinError(i)<<endl;
130  }
131  }
132  }
133 
134  if (diagnostics == true)
135  {
136  cout<<"\n";
137  }
138 
139 
140  return h1;
141 }
142 
143 
144 
145 
146 //__________________________________________________________________________________________________________________________
147 TH2D * Normalize_2D_Histogram (TH2D * h1, double nEvts, const TString& option)
148 {
153  double scaleFactor = 0;
154 
155  for (int i = 0; i <= h1->GetNbinsX(); i++)
156  {
157  for (int j = 0; j <= h1->GetNbinsY(); j++)
158  {
159  if (h1->GetBinContent(i, j) > 0)
160  {
161  double bin_content = h1->GetBinContent(i, j);
162  double bin_error = h1->GetBinError(i);
163 
164  if (option=="probability")
165  scaleFactor = 1/(nEvts);
166  else
167  scaleFactor = 1/(nEvts * h1->GetXaxis()->GetBinWidth(i) * h1->GetYaxis()->GetBinWidth(j));
168 
169  h1->SetBinContent(i, j, scaleFactor*bin_content);
170  h1->SetBinError (i, j, scaleFactor*bin_error);
171  }
172  }
173  }
174 
175  return h1;
176 }
177 
178 
179 
180 
181 //__________________________________________________________________________________________________________________________
182 void CorrectTH1DErrors (TH1D * h)
183 {
184  for (int j = 0; j <= h->GetNbinsX(); j++)
185  {
186  double bin_content = h->GetBinContent(j);
187  double bin_error = GetErrorN(bin_content)[0];
188  h->SetBinError(j, bin_error);
189  }
190 }
191 
192 
193 
194 
195 //__________________________________________________________________________________________________________________________
196 void CorrectTH2DErrors (TH2D * h)
197 {
198  for (int i = 0; i <= h->GetNbinsX(); i++)
199  {
200  for (int j = 0; j <= h->GetNbinsY(); j++)
201  {
202  double bin_content = h->GetBinContent(i, j);
203  double bin_error = GetErrorN(bin_content)[0];
204  h->SetBinError(i, j, bin_error);
205  }
206  }
207 }
208 
209 
210 
211 
212 //__________________________________________________________________________________________________________________________
213 TH2D * GetCleanedUpHistogram (TH2D * h1, double MinBinCount = 1)
214 {
215  h1->SetName(Form("%s_cleaned", h1->GetName()));
216 
217  for (int i = 0; i <= h1->GetNbinsX(); i++)
218  {
219  for (int j = 0; j <= h1->GetNbinsY(); j++)
220  {
221  if (h1->GetBinContent(i, j) < MinBinCount)
222  {
223  h1->SetBinContent(i, j, 0);
224  h1->SetBinError(i, j, 0);
225  }
226  }
227  }
228 
229  return h1;
230 }
231 
232 
233 
234 
235 //__________________________________________________________________________________________________________________________
236 TH2D * GetCleanedUpHistogram (const TH2D * h, double MinBinCount = 1)
237 {
238  TH2D * h1 = (TH2D*)h->Clone();
239  h1->SetName(Form("%s_cleaned", h->GetName()));
240 
241  for (int i = 0; i <= h1->GetNbinsX(); i++)
242  {
243  for (int j = 0; j <= h1->GetNbinsY(); j++)
244  {
245  if (h1->GetBinContent(i, j) < MinBinCount)
246  {
247  h1->SetBinContent(i, j, 0);
248  h1->SetBinError(i, j, 0);
249  }
250  }
251  }
252 
253  return h1;
254 }
255 
256 
257 
258 
259 //__________________________________________________________________________________________________________________________
260 TH2D * GetIntegratedSpectraHistogram (TH2D * h1_original, int yMin, int yMax)
261 {
262  //____________________________________________
263  //make backup copy
264  TH2D * h1 = (TH2D*)h1_original->Clone();
265 
266  for (int i = 0; i <= h1->GetNbinsX(); i++)
267  {
268  for (int j = 0; j <= h1->GetNbinsY(); j++)
269  {
270  if (j < yMin || j > yMax)
271  {
272  h1->SetBinContent(i, j, 0);
273  h1->SetBinError(i, j, 0);
274  }
275  }
276  }
277 
278  return h1;
279 }
280 
281 
282 
283 
284 //__________________________________________________________________________________________________________________________
285 double GetIntegralTH2D (const TH2D* h1)
286 {
287  double integral = 0;
288  for (int i = 1; i <= h1->GetNbinsX(); i++)
289  {
290  for (int j = 1; j <= h1->GetNbinsY(); j++)
291  {
292  integral = integral + h1->GetBinContent(i, j);
293  }
294  }
295 
296  return integral;
297 }
298 
299 
300 
301 
302 //________________________________________________________________________________________________________
303 void printVector(const vector<int> &vec)
304 {
305  for (double x : vec)
306  cout<<x<<' ';
307  cout<<endl;
308 }
309 
310 
311 
312 
313 //__________________________________________________________________________________________________________________________
314 std::vector<int> GetPPtFromVolumeId(int volumeID, int ixx_MAX, bool debug = false)
315 {
316  int ipt = (int)volumeID/(int)ixx_MAX + 1;
317  int ipp = volumeID % ixx_MAX;
318 
319  if (ipp == 0)
320  {
321  ipp = ixx_MAX;
322  ipt = ipt - 1;
323  }
324 
325  if (debug == true) cerr<<"\nVerify (ipp, ipt) = "<<ipp<<", "<<ipt<<endl;
326 
327  std::vector<int> ipp_ipt = {ipp, ipt};
328  return ipp_ipt;
329 }
330 
331 
332 
333 
334 //__________________________________________________________________________________________________________________________
335 TH2D * MultiplyWithCorrectionFactor(const TH2D * h1_original, const TH2D * hist_corrections, double MAX_CORRECTION = 50.0)
336 {
337  //____________________________________________
338  //make backup copy
339  TH2D * h1 = (TH2D*)h1_original->Clone();
340 
341  for (int i = 0; i <= h1->GetNbinsX(); i++)
342  {
343  for (int j = 0; j <= h1->GetNbinsY(); j++)
344  {
345  if (h1->GetBinContent(i, j) > 0 && hist_corrections->GetBinContent(i, j) < MAX_CORRECTION)
346  {
347  if(h1->GetXaxis()->GetBinCenter(i) != hist_corrections->GetXaxis()->GetBinCenter(i))
348  cerr<<"\n[ERROR]: x bin centers do not match!";
349  if(h1->GetYaxis()->GetBinCenter(j) != hist_corrections->GetYaxis()->GetBinCenter(j))
350  cerr<<"\n[ERROR]: y bin centers do not match!";
351  if(h1->GetXaxis()->GetBinWidth(i) != hist_corrections->GetXaxis()->GetBinWidth(i))
352  cerr<<"\n[ERROR]: x bin widths do not match!";
353  if(h1->GetYaxis()->GetBinWidth(j) != hist_corrections->GetYaxis()->GetBinWidth(j))
354  cerr<<"\n[ERROR]: y bin widths do not match!";
355 
356 
357  double new_bin_content = h1->GetBinContent(i, j) * hist_corrections->GetBinContent(i, j);
358 
359  // double new_bin_error = h1->GetBinError(i, j) * hist_corrections->GetBinContent(i, j);
360  double new_bin_error = GetErrorN(h1->GetBinContent(i, j))[0] * hist_corrections->GetBinContent(i, j);
361 
362  h1->SetBinContent(i, j, new_bin_content);
363  h1->SetBinError(i, j, new_bin_error);
364  }
365 
366  else
367  {
368  h1->SetBinContent(i, j, 0);
369  h1->SetBinError(i, j, 0);
370  }
371  }
372  }
373 
374  return h1;
375 }
376 
377 
378 
379 
380 
381 //__________________________________________________________________________________________________________________________
382 // TGraph * RemoveZeroFromGraph(const TGraph * G, bool diagnostics = false)
383 
384 template <class TGraphAll>
385 TGraphAll * RemoveZeroFromGraph (const TGraphAll * G, bool diagnostics = false)
386 {
387  // template <class SomeType>
388  // SomeType sum (SomeType a, SomeType b)
389  // {
390  // return a+b;
391  // }
392 
393  TGraphAll * G_with_zeros = (TGraphAll*)G->Clone();
394 
395  // Remove points whose y value is zero in G_with_zeros
396  int n = G_with_zeros->GetN();
397  Double_t ax[n],ay[n];
398  vector<int> point_indices_to_delete;
399  for (int i_bin = 0; i_bin < n; i_bin++)
400  {
401  G_with_zeros->GetPoint(i_bin,ax[i_bin],ay[i_bin]);
402  if (diagnostics == true)
403  cout<<i_bin<<", "<<ax[i_bin]<<", "<<ay[i_bin]<<endl;
404 
405  if (ay[i_bin] == 0)
406  {
407  if (diagnostics == true)
408  cout<<"Bin to delete: "<<i_bin<<endl;
409 
410  point_indices_to_delete.push_back(i_bin);
411  }
412  }
413 
414  for (int k = 0; k < int(point_indices_to_delete.size()); k++)
415  {
416  if (diagnostics == true)
417  cerr<<"Deleting bin "<<point_indices_to_delete[k];
418 
419  G_with_zeros->RemovePoint(point_indices_to_delete[k]); if (diagnostics == true) cerr<<" ...done."<<endl;
420  std::transform(std::begin(point_indices_to_delete),std::end(point_indices_to_delete),std::begin(point_indices_to_delete),[](int x){return x-1;});
421  }
422 
423  if (diagnostics == true)
424  {
425  cout<<endl;
426  n = G_with_zeros->GetN();
427  Double_t cx[n],cy[n];
428  for (int i_bin = 0; i_bin < n; i_bin++)
429  {
430  G_with_zeros->GetPoint(i_bin,cx[i_bin],cy[i_bin]);
431  cout<<i_bin<<", "<<cx[i_bin]<<", "<<cy[i_bin]<<endl;
432  }
433  cout<<endl<<endl<<endl;
434  }
435 
436  return G_with_zeros;
437 }
438 
439 
440 
441 
442 //__________________________________________________________________________________________________________________________
443 TGraphErrors * SetTGraphXErrorZero (const TGraphErrors * G, bool ShowDiagnostics = false)
444 {
445 
446  TGraphErrors * G_with_errorX_zero = (TGraphErrors*)G->Clone();
447 
448  // Remove points whose y value is zero in G_with_zeros
449  int n = G_with_errorX_zero->GetN();
450  Double_t ax[n], ay[n];
451 
452  for (int i_bin = 0; i_bin < n; i_bin++)
453  {
454  if (ShowDiagnostics) {};
455 
456  G_with_errorX_zero->GetPoint(i_bin, ax[i_bin], ay[i_bin]);
457  G_with_errorX_zero->SetPointError(i_bin, 0, G_with_errorX_zero->GetErrorY(i_bin));
458  }
459 
460  return G_with_errorX_zero;
461 }
462 
463 
464 
465 
466 //__________________________________________________________________________________________________________________________
467 TGraphErrors * MergeSinglePointTGraphErrors (vector<TGraphErrors*> vGFitConfidenceBands)
468 {
469  //each graph in the vector of graphs has a single point in it.
470 
471  int n = vGFitConfidenceBands.size();
472  Double_t ax[n], ay[n], ex[n], ey[n];
473 
474  for (int i = 0; i < n; i++)
475  {
476  vGFitConfidenceBands[i]->GetPoint(0, ax[i], ay[i]);
477  ex[i] = 0;
478  ey[i] = vGFitConfidenceBands[i]->GetErrorY(0);
479 
480  cout<<"Merged graph input: "<<ax[i]<<", "<<ay[i]<<", "<<ex[i]<<", "<<ey[i]<<endl;
481  }
482 
483  auto gr = new TGraphErrors(n, ax, ay, ex, ey);
484  return gr;
485 }
486 
487 
488 
489 
490 
491 
492 //__________________________________________________________________________________________________________________________
493 void AddMissingPointAndErrorToGraphFromFit (TGraphErrors * G,
494  const TFitResultPtr fit_ptr,
495  const TF1 * FitD,
496  vector<int> &vMissingPointIndex,
497  bool debugMode = false)
498 {
499  if (debugMode)
500  cerr<<"\nChecking graph for missing points..."<<endl;
501 
502  vector<double> vXPoints = {0.05, 0.15, 0.25, 0.35, 0.45, 0.55};
503  int n = G->GetN();
504  Double_t ax[n], ay[n];
505 
507  for (unsigned int k = 0; k < vXPoints.size(); k++)
508  {
511  bool point_missing = true;
512 
513  for (int i_bin = 0; i_bin < n; i_bin++)
514  {
515  G->GetPoint(i_bin, ax[i_bin], ay[i_bin]);
516 
518  if (abs(ax[i_bin] - vXPoints[k]) < 0.01)
519  {
520  if (debugMode)
521  cerr<<"x-point "<<ax[i_bin]<<" found."<<endl;
522 
523  point_missing = false;
524  break;
525  }
526  }
527 
529  if (point_missing == true)
530  {
531  if (debugMode)
532  cerr<<"Adding missing point "<<vXPoints[k];
533 
534  vMissingPointIndex.push_back(k);
535 
536  double x[1] = {vXPoints[k]};
537  double err[1]; // error on the function at point x0
538 
540  fit_ptr->GetConfidenceIntervals(1, 1, 1, x, err, 0.683, false);
541 
542  if (debugMode)
543  cerr<<": Function value at "<<x[0]<<" = "<<FitD->Eval(x[0])<<" +/- "<<err[0]<<endl;
544 
545  int temp_n = G->GetN();
546  G->SetPoint(temp_n, vXPoints[k], FitD->Eval(x[0]));
547  G->SetPointError(temp_n, 0, err[0]);
548  }
549  }
550 
551  G->Sort();
552 
553  if (debugMode)
554  cerr<<"Done."<<endl;
555 }
556 
557 
558 
559 
560 //__________________________________________________________________________________________________________________________
561 TH1D * GetRatioTH1D(const TH1D * h1, const TH1D * h2)
562 {
563  if (h1->GetNbinsX() != h2->GetNbinsX())
564  cerr<<"\n\n\n[ERROR]: The bins for calculating the ratio histogram do no match!\n\n"<<endl;
565 
566  TH1D * hUp = (TH1D*)h1->Clone();
567  TH1D * hDown = (TH1D*)h2->Clone();
568 
569  TH1D * hRatio = (TH1D*)hUp->Clone();
570  hRatio->Reset("M");
571 
572  CorrectTH1DErrors(hUp);
573  CorrectTH1DErrors(hDown);
574 
575  for (int i = 0; i <= hUp->GetNbinsX(); i++)
576  {
577  if (hUp->GetBinCenter(i) - hDown->GetBinCenter(i) < 0.01)
578  {
579  if (hDown->GetBinContent(i) > 0)
580  {
581  double ratio = hUp->GetBinContent(i)/hDown->GetBinContent(i);
582  double error_ratio = ratio * sqrt(pow(hUp->GetBinError(i)/hUp->GetBinContent(i), 2) + pow(hDown->GetBinError(i)/hDown->GetBinContent(i), 2));
583 
584  hRatio->SetBinContent(i, ratio);
585  hRatio->SetBinError(i, error_ratio);
586  }
587  else
588  {
589  hRatio->SetBinContent(i, 0);
590  hRatio->SetBinError(i, 0);
591  }
592  }
593 
594  // cout<<hRatio->GetBinContent(i)<<endl;
595  }
596 
597  delete hUp;
598  delete hDown;
599 
600  return hRatio;
601 }
602 
603 
604 
605 //__________________________________________________________________________________________________________________________
606 TH2D * GetRatioTH2D(const TH2D * hNumerator, const TH2D * hDenominator)
607 {
609 
610  TH2D * hN = (TH2D*)hNumerator->Clone();
611  TH2D * hD = (TH2D*)hDenominator->Clone();
612 
613  TH2D * hRatio = (TH2D*)hN->Clone();
614  hRatio->Reset();
615  hRatio->GetZaxis()->SetTitle("Ratio");
616 
617  // sanity check
618  int n_xbins_hN = hN->GetXaxis()->GetNbins();
619  int n_xbins_hD = hD->GetXaxis()->GetNbins();
620  int n_ybins_hN = hN->GetYaxis()->GetNbins();
621  int n_ybins_hD = hD->GetYaxis()->GetNbins();
622  if (n_xbins_hN != n_xbins_hD || n_ybins_hN != n_ybins_hD)
623  {
624  cerr<<"\n\n\n[ERROR]: nXBins or nYBins do not match!"<<endl;
625  return hRatio;
626  }
627 
628 
629  for (int i = 1; i < n_xbins_hN ; i++) //loop over x-axis
630  {
631  for (int j = 1; j < n_ybins_hN ; j++) //loop over y-axis
632  {
633  //sanity check #1 on x-bin-center
634  if (hN->GetXaxis()->GetBinCenter(i) != hD->GetXaxis()->GetBinCenter(i))
635  {
636  cerr<<"\n\n\n[ERROR]: bin centers for bin"<<i<<" in the x-axis do not match!"<<endl;
637  return hRatio;
638  }
639 
640  if (hN->GetYaxis()->GetBinCenter(j) != hD->GetYaxis()->GetBinCenter(j))
641  {
642  cerr<<"\n\n\n[ERROR]: bin centers for bin"<<j<<" in the y-axis do not match!"<<endl;
643  return hRatio;
644  }
645 
646 
647  // set bin content to ratio if possible
648  if (hD->GetBinContent(i, j) != 0)
649  hRatio->SetBinContent(i, j, double(hN->GetBinContent(i, j))/(double)hD->GetBinContent(i, j));
650  }
651  }
652 
653  return hRatio;
654 }
655 
656 
657 
658 
659 
660 
661 //__________________________________________________________________________________________________________________________
662 TGraphErrors * GetdndyGraphFromFit (const TFitResultPtr fit_ptr,
663  const TF1 * FitD,
664  const TGraphErrors * DataGraph,
665  TGraphErrors * GConfidenceBand,
666  int nXPoints,
667  const int ySeries,
668  bool debugMode = false)
669 {
670  if (debugMode)
671  cerr<<"\n\nChecking graph for missing points..."<<endl;
672 
673  vector<double> vXPoints; //= {0.05, 0.15, 0.25, 0.35, 0.45, 0.55};
674  for (int k = 0; k < nXPoints; k++)
675  {
676  vXPoints.push_back(0.05 + 0.1*k);
677  }
678 
679 
680  double integrated_y = 0;
681  double y_fit_errSq = 0;
682  double y_stat_errSq = 0;
683  double pT_bin_width = 0.1;
684 
686  for (unsigned int k = 0; k < vXPoints.size(); k++)
687  {
688  if (debugMode)
689  cerr<<"Adding point "<<vXPoints[k];
690 
691  double x[1] = {vXPoints[k]};
692  double fit_err[1]; // error on the function at point x0
693 
695  fit_ptr->GetConfidenceIntervals(1, 1, 1, x, fit_err, 0.683, false);
696 
697  if (debugMode)
698  cerr<<": Function value at "<<x[0]<<" = "<<FitD->Eval(x[0])<<" +/- "<<fit_err[0]<<endl;
699 
700  integrated_y = integrated_y + pT_bin_width*FitD->Eval(x[0]);
701  y_fit_errSq = y_fit_errSq + fit_err[0]*fit_err[0];
702  if (DataGraph->GetErrorY(k) > 0)
703  y_stat_errSq = y_stat_errSq + pT_bin_width*DataGraph->GetErrorY(k)*DataGraph->GetErrorY(k);
704 
705  // cerr<<"\n"<<DataGraph->GetPointY(k)<<", "<<DataGraph->GetErrorY(k)<<endl;
706  }
707 
708  cerr<<"y = "<<integrated_y<<endl;
709  cerr<<"y_fit_err = "<<sqrt(y_fit_errSq)<<endl;
710  cerr<<"y_stat_err = "<<sqrt(y_stat_errSq)<<endl;
711 
712 
713  // for graph with fit confidence band
714  GConfidenceBand->SetPoint(0, -1.1+2*double(ySeries)/10, integrated_y);
715  GConfidenceBand->SetPointError(0, 0, sqrt(y_fit_errSq));
716 
717  // to return graph with statistical error
718  Double_t x[1], y[1], ex[1], ey[1];
719  x[0] = -1.1+2*double(ySeries)/10;
720  y[0] = integrated_y;
721  ex[0] = 0;
722  ey[0] = sqrt(y_stat_errSq);
723  auto gr = new TGraphErrors(1, x, y, ex, ey);
724  return gr;
725 }
726 
727 
728 
729 
730 
731 
732 
733 
734 //__________________________________________________________________________________________________________________________
735 vector<TGraphErrors*> GetGraphsfromTH2D(const TH2D * h,
736  double year,
737  const string &particle_name,
738  int rapidity_bin_max_data,
739  int rapidity_bin_min_data = 11,
740  bool doScaling = true)
741 {
742  vector<TGraphErrors*> dataGraphs;
743 
744 
745 
746  for (int yBin = rapidity_bin_min_data; yBin <= rapidity_bin_max_data; yBin++) //along rapidity
747  {
748  auto htemp = h->ProjectionY(Form("py_bin_%d", yBin), yBin, yBin);
749  // htemp->SaveAs(Form("py_bin_%d.root", yBin));
750  for (int pT_bin = 1; pT_bin <= htemp->GetNbinsX(); pT_bin++) //along pT
751  {
752  if (doScaling == true)
753  {
754  htemp->SetBinContent(pT_bin, pow (10, -(yBin - rapidity_bin_min_data))*htemp->GetBinContent(pT_bin));
755  htemp->SetBinError(pT_bin, pow (10, -(yBin - rapidity_bin_min_data))*htemp->GetBinError(pT_bin));
756  }
757 
758 
759  //_____________________________________________________
760  if (particle_name=="kMinus")
761  {
762  if (year == 2010)
763  {
764  if (pT_bin > 15)
765  htemp->SetBinContent(pT_bin, 0);
766 
767  if (yBin > 18 && pT_bin > 13)
768  htemp->SetBinContent(pT_bin, 0);
769 
770  if (yBin > 19 && pT_bin > 9)
771  htemp->SetBinContent(pT_bin, 0);
772  }
773  }
774 
775 
776  //_____________________________________________________
777  else if (particle_name=="kPlus")
778  {
779  if (year == 2010)
780  {
781  if (pT_bin > 15)
782  htemp->SetBinContent(pT_bin, 0);
783 
784  if (yBin > 18 && pT_bin > 13)
785  htemp->SetBinContent(pT_bin, 0);
786 
787  if (yBin > 19 && pT_bin > 9)
788  htemp->SetBinContent(pT_bin, 0);
789  }
790  }
791 
792 
793  //_____________________________________________________
794  else if (particle_name=="piMinus")
795  {
796  if (year == 2010)
797  {
798  if (pT_bin > 15)
799  htemp->SetBinContent(pT_bin, 0);
800 
801  if (yBin == 11 && pT_bin < 5)
802  htemp->SetBinContent(pT_bin, 0);
803 
804  if (yBin == 12 && pT_bin < 4)
805  htemp->SetBinContent(pT_bin, 0);
806 
807  if (yBin == 13 && pT_bin < 3)
808  htemp->SetBinContent(pT_bin, 0);
809 
810  if (yBin == 14 && pT_bin < 3)
811  htemp->SetBinContent(pT_bin, 0);
812 
813  if (yBin == 15 && pT_bin < 2)
814  htemp->SetBinContent(pT_bin, 0);
815 
816  if (yBin == 20 && pT_bin > 12)
817  htemp->SetBinContent(pT_bin, 0);
818 
819  if (yBin == 21 && pT_bin > 11)
820  htemp->SetBinContent(pT_bin, 0);
821 
822  if (yBin == 22 && pT_bin > 10)
823  htemp->SetBinContent(pT_bin, 0);
824  }
825  }
826 
827 
828  //_____________________________________________________
829  else if (particle_name=="piPlus")
830  {
831  if (year == 2010)
832  {
833  if (pT_bin > 15)
834  htemp->SetBinContent(pT_bin, 0);
835 
836  if (yBin == 11 && pT_bin < 5)
837  htemp->SetBinContent(pT_bin, 0);
838 
839  if (yBin == 12 && pT_bin < 4)
840  htemp->SetBinContent(pT_bin, 0);
841 
842  if (yBin == 13 && pT_bin < 3)
843  htemp->SetBinContent(pT_bin, 0);
844 
845  if (yBin == 14 && pT_bin < 3)
846  htemp->SetBinContent(pT_bin, 0);
847 
848  if (yBin == 15 && pT_bin < 2)
849  htemp->SetBinContent(pT_bin, 0);
850 
851  if (yBin == 21 && pT_bin > 12)
852  htemp->SetBinContent(pT_bin, 0);
853 
854  if (yBin == 22 && pT_bin > 11)
855  htemp->SetBinContent(pT_bin, 0);
856  }
857  }
858 
859 
860  //_____________________________________________________
861  else if (particle_name=="proton")
862  {
863  if (year == 2010)
864  {
865  if (pT_bin > 15)
866  htemp->SetBinContent(pT_bin, 0);
867 
868  // if (yBin > 18 && pT_bin > 13)
869  // htemp->SetBinContent(pT_bin, 0);
870 
871  // if (yBin > 19)
872  // htemp->SetBinContent(pT_bin, 0);
873  }
874  }
875 
876 
877  //_____________________________________________________
878  else if (particle_name=="antiproton")
879  {
880  if (year == 2010)
881  {
882  if (pT_bin > 15)
883  htemp->SetBinContent(pT_bin, 0);
884 
885  if (yBin > 18 && pT_bin > 13)
886  htemp->SetBinContent(pT_bin, 0);
887 
888  if (yBin > 19 && pT_bin > 9)
889  htemp->SetBinContent(pT_bin, 0);
890  }
891  }
892 
893 
894  //_____________________________________________________
895  else if (particle_name=="deuteron")
896  {
897  if (year == 2010)
898  {
899  if (pT_bin > 6) //10
900  htemp->SetBinContent(pT_bin, 0);
901  }
902  }
903 
904 
905  //_____________________________________________________
906  else
907  {
908  cerr<<"\n\n\nERROR in GetGraphsfromTH2D(): string "<<particle_name<<" does not match any known name"<<endl;
909  }
910 
911  }
912  TGraphErrors * Gtemp = new TGraphErrors (htemp);
913  TGraphErrors * GtempZero = SetTGraphXErrorZero(Gtemp);
914  delete Gtemp;
915  dataGraphs.push_back(GtempZero);
916  }
917 
918  return dataGraphs;
919 }
920 
921 
922 
923 
924 //__________________________________________________________________________________________________________________________
925 TLegend * GetTLegend(const string &particle_name, double year, const string &legend_type)
926 {
927  TLegend * leg = new TLegend();
928 
929  //_____________________________________________________
930  if (particle_name=="kMinus")
931  {
932  if (legend_type=="points" && year == 2009) leg = new TLegend(0.70,0.20,0.94,0.91); //x1, y1, x2, y2
933  else if (legend_type=="points" && year == 2010) leg = new TLegend(0.70,0.15,0.94,0.84);
934  else if (legend_type=="header" && year == 2009) leg = new TLegend(0.60,0.90,0.95,0.97);
935  else if (legend_type=="header" && year == 2010) leg = new TLegend(0.60,0.90,0.95,0.97);
936  else if (legend_type=="ratiop" && year == 2009) leg = new TLegend(0.70,0.07,0.90,0.86);
937  else if (legend_type=="ratiop" && year == 2010) leg = new TLegend(0.70,0.02,0.90,0.75);
938  // else if (legend_type=="headerWithEPOS" && year == 2009) leg = new TLegend(0.49,0.92,0.99,0.99);
939  // else if (legend_type=="headerWithEPOS" && year == 2010) leg = new TLegend(0.49,0.92,0.99,0.99);
940  }
941 
942 
943  //_____________________________________________________
944  else if (particle_name=="kPlus")
945  {
946  if (legend_type=="points" && year == 2009) leg = new TLegend(0.70,0.15,0.94,0.89); //x1, y1, x2, y2
947  else if (legend_type=="points" && year == 2010) leg = new TLegend(0.70,0.17,0.94,0.86);
948  else if (legend_type=="header" && year == 2009) leg = new TLegend(0.60,0.90,0.95,0.97);
949  else if (legend_type=="header" && year == 2010) leg = new TLegend(0.60,0.90,0.95,0.97);
950  else if (legend_type=="ratiop" && year == 2009) leg = new TLegend(0.70,0.07,0.90,0.86);
951  else if (legend_type=="ratiop" && year == 2010) leg = new TLegend(0.70,0.02,0.90,0.75);
952  // else if (legend_type=="headerWithEPOS" && year == 2009) leg = new TLegend(0.49,0.92,0.99,0.99);
953  // else if (legend_type=="headerWithEPOS" && year == 2010) leg = new TLegend(0.49,0.92,0.99,0.99);
954  }
955 
956 
957  //_____________________________________________________
958  else if (particle_name=="piMinus")
959  {
960  if (legend_type=="points" && year == 2009) leg = new TLegend(0.70,0.15,0.94,0.89); //x1, y1, x2, y2
961  else if (legend_type=="points" && year == 2010) leg = new TLegend(0.70,0.12,0.94,0.82);
962  else if (legend_type=="header" && year == 2009) leg = new TLegend(0.60,0.90,0.95,0.97); //0.65,0.92,0.95,0.99
963  else if (legend_type=="header" && year == 2010) leg = new TLegend(0.60,0.90,0.95,0.97); //0.65,0.92,0.95,0.99
964  else if (legend_type=="ratiop" && year == 2009) leg = new TLegend(0.70,0.07,0.90,0.86);
965  else if (legend_type=="ratiop" && year == 2010) leg = new TLegend(0.70,0.02,0.90,0.75);
966  // else if (legend_type=="headerWithEPOS" && year == 2009) leg = new TLegend(0.49,0.92,0.99,0.99);
967  // else if (legend_type=="headerWithEPOS" && year == 2010) leg = new TLegend(0.49,0.92,0.99,0.99);
968  }
969 
970 
971  //_____________________________________________________
972  else if (particle_name=="piPlus")
973  {
974  if (legend_type=="points" && year == 2009) leg = new TLegend(0.70,0.15,0.94,0.90); //x1, y1, x2, y2
975  else if (legend_type=="points" && year == 2010) leg = new TLegend(0.70,0.12,0.94,0.83);
976  else if (legend_type=="header" && year == 2009) leg = new TLegend(0.60,0.90,0.95,0.97); //0.65,0.92,0.95,0.99
977  else if (legend_type=="header" && year == 2010) leg = new TLegend(0.60,0.90,0.95,0.97); //0.65,0.92,0.95,0.99
978  else if (legend_type=="ratiop" && year == 2009) leg = new TLegend(0.70,0.07,0.90,0.86);
979  else if (legend_type=="ratiop" && year == 2010) leg = new TLegend(0.70,0.02,0.90,0.75);
980  // else if (legend_type=="headerWithEPOS" && year == 2009) leg = new TLegend(0.49,0.92,0.99,0.99);
981  // else if (legend_type=="headerWithEPOS" && year == 2010) leg = new TLegend(0.49,0.92,0.99,0.99);
982  }
983 
984 
985  //_____________________________________________________
986  else if (particle_name=="proton")
987  {
988  if (legend_type=="points" && year == 2009) leg = new TLegend(0.70,0.21,0.94,0.91); //x1, y1, x2, y2
989  else if (legend_type=="points" && year == 2010) leg = new TLegend(0.70,0.22,0.94,0.87); //(0.70,0.14,0.94,0.87);
990  else if (legend_type=="header" && year == 2009) leg = new TLegend(0.60,0.90,0.95,0.97); //0.65,0.92,0.95,0.99
991  else if (legend_type=="header" && year == 2010) leg = new TLegend(0.60,0.90,0.95,0.97); //0.65,0.92,0.95,0.99
992  else if (legend_type=="ratiop" && year == 2009) leg = new TLegend(0.70,0.07,0.90,0.86);
993  else if (legend_type=="ratiop" && year == 2010) leg = new TLegend(0.70,0.02,0.90,0.75);
994  // else if (legend_type=="headerWithEPOS" && year == 2009) leg = new TLegend(0.49,0.92,0.99,0.99);
995  // else if (legend_type=="headerWithEPOS" && year == 2010) leg = new TLegend(0.49,0.92,0.99,0.99);
996  else if (legend_type=="stackedspectra" && year == 2010) leg = new TLegend(0.70,0.40,0.90,0.80);
997  else if (legend_type=="ExtraPoints" && year == 2010) leg = new TLegend(0.65,0.13,0.96,0.90);
998  }
999 
1000 
1001  //_____________________________________________________
1002  else if (particle_name=="antiproton")
1003  {
1004  if (legend_type=="points" && year == 2009) leg = new TLegend(0.70,0.18,0.94,0.88); //x1, y1, x2, y2
1005  else if (legend_type=="points" && year == 2010) leg = new TLegend(0.70,0.13,0.94,0.81);
1006  // else if (legend_type=="points" && year == 2010) leg = new TLegend(0.70,0.12,0.94,0.85); //for no Szymon comparison
1007  else if (legend_type=="header" && year == 2009) leg = new TLegend(0.60,0.90,0.95,0.97); //0.60,0.88,0.95,0.95
1008  else if (legend_type=="header" && year == 2010) leg = new TLegend(0.60,0.90,0.95,0.97); //0.60,0.88,0.95,0.95
1009  else if (legend_type=="ratiop" && year == 2009) leg = new TLegend(0.70,0.07,0.90,0.86);
1010  else if (legend_type=="ratiop" && year == 2010) leg = new TLegend(0.70,0.02,0.90,0.75);
1011  // else if (legend_type=="headerWithEPOS" && year == 2009) leg = new TLegend(0.49,.92,0.99,0.99);
1012  // else if (legend_type=="headerWithEPOS" && year == 2010) leg = new TLegend(0.49,.92,0.99,0.99);
1013  }
1014 
1015 
1016  //_____________________________________________________
1017  else if (particle_name=="deuteron")
1018  {
1019  if (legend_type=="points" && year == 2009) leg = new TLegend(0.70,0.13,0.94,0.90); //x1, y1, x2, y2
1020  else if (legend_type=="points" && year == 2010) leg = new TLegend(0.70,0.13,0.94,0.90);
1021  // else if (legend_type=="points" && year == 2010) leg = new TLegend(0.70,0.12,0.94,0.85); //for no Szymon comparison
1022  else if (legend_type=="header" && year == 2009) leg = new TLegend(0.60,0.90,0.95,0.97); //0.60,0.88,0.95,0.95
1023  else if (legend_type=="header" && year == 2010) leg = new TLegend(0.60,0.90,0.95,0.97); //0.60,0.88,0.95,0.95
1024  else if (legend_type=="ratiop" && year == 2009) leg = new TLegend(0.70,0.07,0.90,0.86);
1025  else if (legend_type=="ratiop" && year == 2010) leg = new TLegend(0.70,0.02,0.90,0.75);
1026  // else if (legend_type=="headerWithEPOS" && year == 2009) leg = new TLegend(0.49,.92,0.99,0.99);
1027  // else if (legend_type=="headerWithEPOS" && year == 2010) leg = new TLegend(0.49,.92,0.99,0.99);
1028  else if (legend_type=="stackedspectra" && year == 2010) leg = new TLegend(0.70,0.40,0.90,0.80);
1029  else if (legend_type=="ExtraPoints" && year == 2010) leg = new TLegend(0.65,0.13,0.96,0.90);
1030 
1031  // else if (legend_type=="p0" && year == 6) leg = new TLegend(0.80,0.52,0.97,0.57); //x1, y1, x2, y2
1032  // else if (legend_type=="p0" && year == 5) leg = new TLegend(0.80,0.43,0.97,0.48);
1033  // else if (legend_type=="p0" && year == 4) leg = new TLegend(0.80,0.36,0.97,0.41);
1034  // else if (legend_type=="p0" && year == 3) leg = new TLegend(0.80,0.29,0.97,0.34);
1035  // else if (legend_type=="p0" && year == 2) leg = new TLegend(0.80,0.24,0.97,0.29);
1036  // else if (legend_type=="p0" && year == 1) leg = new TLegend(0.80,0.20,0.97,0.25);
1037  // else if (legend_type=="p0" && year == 0) leg = new TLegend(0.80,0.16,0.97,0.21);
1038  else if (legend_type=="p0" && year == 6) leg = new TLegend(0.80,0.82,0.97,0.87); //x1, y1, x2, y2
1039  else if (legend_type=="p0" && year == 5) leg = new TLegend(0.80,0.67,0.97,0.72);
1040  else if (legend_type=="p0" && year == 4) leg = new TLegend(0.80,0.54,0.97,0.59);
1041  else if (legend_type=="p0" && year == 3) leg = new TLegend(0.80,0.44,0.97,0.48);
1042  else if (legend_type=="p0" && year == 2) leg = new TLegend(0.80,0.35,0.97,0.40);
1043  else if (legend_type=="p0" && year == 1) leg = new TLegend(0.80,0.27,0.97,0.32);
1044  else if (legend_type=="p0" && year == 0) leg = new TLegend(0.80,0.21,0.97,0.26);
1045 
1046  // else if (legend_type=="p0" && year == 2010) leg = new TLegend(0.75,0.52,0.95,0.70); //x1, y1, x2, y2
1047  else if (legend_type=="p0" && year == 2010) leg = new TLegend(0.75,0.85,0.95,0.97); //x1, y1, x2, y2
1048 
1049  else if (legend_type=="p0percent" && year == 2010) leg = new TLegend(0.75,0.60,0.94,0.95); //x1, y1, x2, y2
1050  }
1051 
1052 
1053  //_____________________________________________________
1054  else
1055  {
1056  cerr<<"\n\n\nERROR in GetTLegend(): string "<<particle_name<<" does not match any known name"<<endl;
1057  }
1058 
1059 
1060  if (legend_type=="headerWithEPOS" && year == 2009) leg = new TLegend(0.22,0.845,0.955,0.995); //x1, y1, x2, y2
1061  if (legend_type=="headerWithEPOS" && year == 2010) leg = new TLegend(0.22,0.845,0.955,0.995); //x1, y1, x2, y2
1062 
1063 
1064  return leg;
1065 }
1066 
1067 
1068 
1069 
1070 //__________________________________________________________________________________________________________________________
1071 vector<TH1D*> GetEPOSYGraphs(const string &particle_name, const vector<double> &vP0)
1072 {
1073  TFile * fEPOSpp158 = new TFile ("inputfiles_EPOSAfterburner/pp_158_combo.root", "READ");
1074  if (fEPOSpp158->IsOpen())
1075  cout<<"\n[INFO] Opened MC flatphasespace file "<<fEPOSpp158->GetName()<<" in READ mode."<<endl;
1076 
1077  auto hEvents = (TH1D*)fEPOSpp158->Get("Hist_header");
1078  double nEvents = hEvents->GetEntries();
1079 
1080  vector<TH1D*> vHistEPOS;
1081 
1082  for (unsigned int k = 0; k < vP0.size(); k++)
1083  {
1084  TH1D * Hist_y_deuteron = (TH1D*)fEPOSpp158->Get(Form("Hist_y_%s_dp0_%.2f", particle_name.c_str(), vP0[k]));
1085  Hist_y_deuteron->Rebin(2);
1086  Hist_y_deuteron = Normalize_1D_Histogram(Hist_y_deuteron, nEvents);
1087 
1088  vHistEPOS.push_back(Hist_y_deuteron);
1089  }
1090 
1091  return vHistEPOS;
1092 }
1093 
1094 
1095 
1096 //__________________________________________________________________________________________________________________________
1097 vector<TGraphErrors*> GetEPOS_Y_PT_Spectra_Graphs(const string &particle = "deuteron",
1098  const double p0 = 84.19,
1099  int rebinY = 2,
1100  bool doScaling = false,
1101  double pT_max_EPOS = 1.0,
1102  int rapidity_bin_max_data = 45,
1103  int rapidity_bin_min_data = 48)
1104 {
1105  //Get 2D y-pT spectra graphs from the saved antihelium EPOS simulation files.
1106 
1113  TFile * fEPOSpp158 = new TFile ("inputfiles_EPOSAfterburner/pp_158_combo.root", "READ");
1114  if (fEPOSpp158->IsOpen())
1115  cerr<<"[INFO] Opened MC flatphasespace file "<<fEPOSpp158->GetName()<<" in READ mode."<<endl;
1116 
1117  TH1D * hEvents = (TH1D*)fEPOSpp158->Get("Hist_header");
1118  double nEvents = hEvents->GetEntries();
1119  cerr<<"[INFO] Number of p+p events simulated with EPOS+Afterburner: "<<nEvents<<endl;
1120 
1121  string histname = "Hist_pT_y_" + particle + "_";
1122  string str_p0 = string(Form("%.2f", p0));
1123  if (particle == "antiproton" || particle == "piMinus" || particle == "kMinus") histname = histname + "dbarp0_" + str_p0;
1124  else if (particle == "proton" || particle == "piPlus" || particle == "kPlus" || particle == "deuteron") histname = histname + "dp0_" + str_p0;
1125  else
1126  cerr<<"\n\n\nERROR in "<<__PRETTY_FUNCTION__<<": string "<<particle<<" does not match any saved name."<<endl;
1127 
1128 
1129  TH2D * h2D_pT_y = (TH2D*)fEPOSpp158->Get(histname.c_str());
1130  h2D_pT_y->RebinY(rebinY); //20, along the rapidity axis
1131  h2D_pT_y = Normalize_2D_Histogram(h2D_pT_y, nEvents, "");
1132 
1133 
1134  string filenamestring = "plots_EPOS/" + particle + "/Y_PT_Spectra_Graphs_" + particle + "_" + str_p0 + "_DataCheck.root";
1135  TFile * fOut = new TFile(filenamestring.c_str(), "RECREATE");
1136  fOut->cd();
1137  h2D_pT_y->Write();
1138 
1139 
1140  cerr<<"[INFO] In "<<histname<<" from "<<fEPOSpp158->GetName()<<", getting y-Bins: ";
1141  vector<TGraphErrors*> EPOSGraphs;
1142 
1143  for (int yBin = rapidity_bin_min_data; yBin <= rapidity_bin_max_data; yBin++) //along rapidity //61
1144  {
1145  auto htemp = h2D_pT_y->ProjectionX(Form("y_bin_%d_y=%.2f", yBin, h2D_pT_y->GetYaxis()->GetBinCenter(yBin)), yBin, yBin);
1146  cerr<<"["<<h2D_pT_y->GetYaxis()->GetBinLowEdge(yBin)<<", "<<h2D_pT_y->GetYaxis()->GetBinLowEdge(yBin)+h2D_pT_y->GetYaxis()->GetBinWidth(yBin)<<"], ";
1147  htemp->Write();
1148 
1149  for (int pT_bin = 1; pT_bin <= htemp->GetNbinsX(); pT_bin++) //along pT
1150  {
1152  if (doScaling == true)
1153  {
1154  if (htemp->GetBinContent(pT_bin) > 0)
1155  htemp->SetBinContent(pT_bin, 1.0 * pow (10, -(yBin - rapidity_bin_min_data))*htemp->GetBinContent(pT_bin));
1156  }
1157 
1159  pT_max_EPOS = 1.5;
1160  if (htemp->GetBinCenter(pT_bin) > pT_max_EPOS)
1161  {
1162  htemp->SetBinContent(pT_bin, 0);
1163  htemp->SetBinError(pT_bin, 0);
1164  }
1165 
1167  // htemp->SetBinError(pT_bin, 0);
1168  }
1169 
1170  htemp->Write();
1171  TGraphErrors * Gtemp = new TGraphErrors(htemp);
1172  Gtemp = RemoveZeroFromGraph(Gtemp);
1173  TGraphErrors * GtempZero = SetTGraphXErrorZero(Gtemp);
1174  EPOSGraphs.push_back(GtempZero);
1175  GtempZero->Write();
1176  }
1177 
1178 
1179  delete h2D_pT_y;
1180  delete hEvents;
1181  delete fEPOSpp158;
1182  delete fOut;
1183  cerr<<"\n[INFO] Graphs returned."<<endl;
1184  return EPOSGraphs;
1185 }
1186 
1187 
1188 
1189 
1190 //__________________________________________________________________________________________________________________________
1191 vector<TGraphErrors*> GetEPOS_Y_PT_Spectra_Graphs(TH2D * h2D_pT_y,
1192  const string &particle = "deuteron",
1193  const double p0 = 84.19,
1194  int rapidity_bin_min_data = 45, //41
1195  int rapidity_bin_max_data = 48, //50
1196  bool doScaling = false,
1197  bool doCleanup = true
1198  )
1199 {
1200  //Get 2D y-pT spectra graphs from the saved antihelium EPOS simulation files.
1201 
1202  string str_p0 = string(Form("%.2f", p0));
1203  string filenamestring = "plots_EPOS/" + particle + "/Y_PT_Spectra_Graphs_" + particle + "_" + str_p0 + ".root";
1204  TFile * fOut = new TFile(filenamestring.c_str(), "RECREATE");
1205  fOut->cd();
1206  h2D_pT_y->Write();
1207 
1208  cerr<<"[INFO] From "<<h2D_pT_y->GetName()<<", getting y-Bins: ";
1209  vector<TGraphErrors*> EPOSGraphs;
1210 
1211  for (int yBin = rapidity_bin_min_data; yBin <= rapidity_bin_max_data; yBin++) //along rapidity //61
1212  {
1213  auto htemp = h2D_pT_y->ProjectionX(Form("y_bin_%d_y=%.2f", yBin, h2D_pT_y->GetYaxis()->GetBinCenter(yBin)), yBin, yBin);
1214  cerr<<"["<<h2D_pT_y->GetYaxis()->GetBinLowEdge(yBin)<<", "<<h2D_pT_y->GetYaxis()->GetBinLowEdge(yBin)+h2D_pT_y->GetYaxis()->GetBinWidth(yBin)<<"], ";
1215  htemp->Write();
1216 
1217  for (int pT_bin = 1; pT_bin <= htemp->GetNbinsX(); pT_bin++) //along pT
1218  {
1220  if (doScaling == true)
1221  {
1222  if (htemp->GetBinContent(pT_bin) > 0)
1223  htemp->SetBinContent(pT_bin, 1.0 * pow (10, -(yBin - rapidity_bin_min_data))*htemp->GetBinContent(pT_bin));
1224  }
1225 
1227  if (doCleanup == true)
1228  {
1229  if (particle == "deuteron" && pT_bin > 15)
1230  {
1231  htemp->SetBinContent(pT_bin, 0);
1232  htemp->SetBinError(pT_bin, 0);
1233  }
1234 
1235  if (particle == "proton" && pT_bin > 15)
1236  {
1237  htemp->SetBinContent(pT_bin, 0);
1238  htemp->SetBinError(pT_bin, 0);
1239  }
1240  }
1241 
1243  // htemp->SetBinError(pT_bin, 0);
1244  }
1245 
1246  htemp->Write();
1247  TGraphErrors * Gtemp = new TGraphErrors(htemp);
1248  Gtemp = RemoveZeroFromGraph(Gtemp);
1249  TGraphErrors * GtempZero = SetTGraphXErrorZero(Gtemp);
1250  EPOSGraphs.push_back(GtempZero);
1251  GtempZero->Write();
1252  }
1253 
1254 
1255  cerr<<"\n[INFO] Graphs returned.\n[INFO] Intermediate histograms saved to "<<fOut->GetName()<<endl;
1256  delete fOut;
1257  return EPOSGraphs;
1258 }
1259 
1260 
1261 
1262 
1263 //__________________________________________________________________________________________________________________________
1264 vector<TGraphAsymmErrors*> GetTGraphAsymm(const TH2D* h,
1265  const TH2D* h_high,
1266  const TH2D* h_low,
1267  const string &particle_name,
1268  int rapidity_bin_min_data = 45,
1269  int rapidity_bin_max_data = 48,
1270  bool doScaling = true,
1271  bool printNumbers = false)
1272 {
1273  vector<TGraphAsymmErrors*> dataGraphs;
1274  TH1D * h_copy = (TH1D*)h->Clone();
1275  TH1D * h_high_copy = (TH1D*)h_high->Clone();
1276  TH1D * h_low_copy = (TH1D*)h_low->Clone();
1277 
1278  int pT_max_bin = 6;
1279 
1280  if (printNumbers == true)
1281  {
1282  fprintf(stderr, "\n\nParticle: %10s\n", particle_name.c_str());
1283  fprintf(stderr," ");
1284  for (int pT_bin = 1; pT_bin <= 15; pT_bin++)
1285  {
1286  fprintf(stderr, "%7s%2.1f", "pT=", 0.1*pT_bin);
1287  }
1288  // fprintf(stderr,"\n");
1289  }
1290 
1291 
1292  for (int yBin = rapidity_bin_min_data; yBin <= rapidity_bin_max_data; yBin++) //along rapidity (y-axis)
1293  {
1294  std::vector<double> x;
1295  std::vector<double> x_low;
1296  std::vector<double> x_high;
1297  std::vector<double> y;
1298  std::vector<double> y_errorLow;
1299  std::vector<double> y_errorHigh;
1300 
1301  if (printNumbers == true)
1302  {
1303  fprintf(stderr, "\nFor y=%2.1f", -1.1+2*double(yBin-rapidity_bin_min_data)/10);
1304  }
1305 
1306  for (int pT_bin = 1; pT_bin <= pT_max_bin; pT_bin++) //along pT_bin (x-axis)
1307  {
1308  x.push_back(h_copy->GetXaxis()->GetBinCenter(pT_bin));
1309  x_low.push_back(0);
1310  x_high.push_back(0);
1311 
1312  std::vector<double> unsorted_bin_Counts = {h_low_copy->GetBinContent(pT_bin, yBin), h_copy->GetBinContent(pT_bin, yBin), h_high_copy->GetBinContent(pT_bin, yBin)};
1313 
1315  if (doScaling == true)
1316  {
1317  for (int k = 0; k < int(unsorted_bin_Counts.size()); k++)
1318  {
1319  if (unsorted_bin_Counts[k] > 0)
1320  unsorted_bin_Counts[k] = 1.0 * pow (10, -(yBin - rapidity_bin_min_data)) * unsorted_bin_Counts[k];
1321  }
1322  }
1323 
1324  if (std::is_sorted(unsorted_bin_Counts.begin(), unsorted_bin_Counts.end()))
1325  {
1326  //cout<<"Bin contents in bin "<<pT_bin<<" are sorted!"<<endl;
1327  }
1328  else
1329  {
1330  //cout<<"Bin contents in bin "<<pT_bin<<" are unsorted! low: "<<unsorted_bin_Counts[0]<<", middle: "<<unsorted_bin_Counts[1]<<", high: "<<unsorted_bin_Counts[2]<<". Sorting!"<<endl;
1331  sort(unsorted_bin_Counts.begin(), unsorted_bin_Counts.end());
1332  }
1333 
1334  y_errorLow.push_back(unsorted_bin_Counts[1]-unsorted_bin_Counts[0]);
1335  y.push_back(unsorted_bin_Counts[1]);
1336  y_errorHigh.push_back(unsorted_bin_Counts[2]-unsorted_bin_Counts[1]);
1337 
1338 
1339  if (printNumbers == true)
1340  {
1341  fprintf(stderr, "%10.2e", unsorted_bin_Counts[1]);
1342  }
1343  }
1344 
1345  TVectorD tv_x(x.size(), &x[0]);
1346  TVectorD tv_x_high(x_high.size(), &x_high[0]);
1347  TVectorD tv_x_low(x_low.size(), &x_low[0]);
1348  TVectorD tv_y(y.size(), &y[0]);
1349  TVectorD tv_y_errorLow(y_errorLow.size(), &y_errorLow[0]);
1350  TVectorD tv_y_errorHigh(y_errorHigh.size(), &y_errorHigh[0]);
1351 
1352  TGraphAsymmErrors * G = new TGraphAsymmErrors(tv_x, tv_y, tv_x_low, tv_x_high, tv_y_errorLow, tv_y_errorHigh);
1353  dataGraphs.push_back(G);
1354 
1355  x.clear();
1356  x_low.clear();
1357  x_high.clear();
1358  y.clear();
1359  y_errorLow.clear();
1360  y_errorHigh.clear();
1361  }
1362 
1363  if (printNumbers == true)
1364  {
1365  fprintf(stderr,"\n\n");
1366  }
1367 
1368  delete h_copy;
1369  delete h_high_copy;
1370  delete h_low_copy;
1371  return dataGraphs;
1372 }
1373 
1374 
1375 
1376 
1377 //__________________________________________________________________________________________________________________________
1378 vector<TGraphAsymmErrors*> GetEPOS_Y_PT_SpectraBand_Graphs(const string &particle = "deuteron",
1379  int rebinY = 2)
1380 {
1381  //Get 2D y-pT spectra band graphs from the saved antihelium EPOS simulation files.
1382 
1383  TFile * fEPOSpp158 = new TFile ("inputfiles_EPOSAfterburner/pp_158_combo.root", "READ");
1384  if (fEPOSpp158->IsOpen())
1385  cout<<"[INFO] Opened MC flatphasespace file "<<fEPOSpp158->GetName()<<" in READ mode."<<endl;
1386 
1387  TH1D * hEvents = (TH1D*)fEPOSpp158->Get("Hist_header");
1388  double nEvents = hEvents->GetEntries();
1389  cout<<"[INFO] Number of p+p events simulated with EPOS+Afterburner: "<<nEvents<<endl;
1390 
1391  string histname = "Hist_pT_y_" + particle + "_";
1392  if (particle == "antiproton" || particle == "piMinus" || particle == "kMinus") histname = histname + "dbarp0_";
1393  else if (particle == "proton" || particle == "piPlus" || particle == "kPlus" || particle == "deuteron") histname = histname + "dp0_";
1394  else
1395  cerr<<"\n\n\nERROR in "<<__PRETTY_FUNCTION__<<": string "<<particle<<" does not match any saved name."<<endl;
1396 
1397  string histname_UPP = histname + "109.45";
1398  cerr<<"[INFO] Getting "<<histname_UPP<<" from "<<fEPOSpp158->GetName()<<endl;
1399  TH2D * h2D_pT_y_UPP = (TH2D*)fEPOSpp158->Get(histname_UPP.c_str());
1400  h2D_pT_y_UPP->RebinY(rebinY); //20, along the rapidity axis
1401  h2D_pT_y_UPP = Normalize_2D_Histogram(h2D_pT_y_UPP, nEvents, "");
1402 
1403  string histname_MID = histname + "84.19";
1404  cerr<<"[INFO] Getting "<<histname_MID<<" from "<<fEPOSpp158->GetName()<<endl;
1405  TH2D * h2D_pT_y_MID = (TH2D*)fEPOSpp158->Get(histname_MID.c_str());
1406  h2D_pT_y_MID->RebinY(rebinY); //20, along the rapidity axis
1407  h2D_pT_y_MID = Normalize_2D_Histogram(h2D_pT_y_MID, nEvents, "");
1408 
1409  string histname_LOW = histname + "58.94";
1410  cerr<<"[INFO] Getting "<<histname_LOW<<" from "<<fEPOSpp158->GetName()<<endl;
1411  TH2D * h2D_pT_y_LOW = (TH2D*)fEPOSpp158->Get(histname_LOW.c_str());
1412  h2D_pT_y_LOW->RebinY(rebinY); //20, along the rapidity axis
1413  h2D_pT_y_LOW = Normalize_2D_Histogram(h2D_pT_y_LOW, nEvents, "");
1414 
1415  vector<TGraphAsymmErrors*> EPOSBandGraphs = GetTGraphAsymm(h2D_pT_y_MID, h2D_pT_y_UPP, h2D_pT_y_LOW, "deuteron");
1416 
1417  delete h2D_pT_y_UPP;
1418  delete h2D_pT_y_MID;
1419  delete h2D_pT_y_LOW;
1420  delete hEvents;
1421  delete fEPOSpp158;
1422  return EPOSBandGraphs;
1423 }
Generated by   doxygen 1.8.5