Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CommonTools.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CommonTools.h
1 void myText(Double_t x,Double_t y,Color_t color,const char *text, Double_t tsize = 0.05, double angle = -1);
2 void DivideCanvas(TVirtualPad *p, int npads);
3 TGraphErrors* FitProfile(const TH2 *h2);
4 
5 void myText(Double_t x,Double_t y,Color_t color,
6  const char *text, Double_t tsize, double angle) {
7 
8  TLatex l; //l.SetTextAlign(12);
9  l.SetTextSize(tsize);
10  l.SetNDC();
11  l.SetTextColor(color);
12  if (angle > 0) l.SetTextAngle(angle);
13  l.DrawLatex(x,y,text);
14 }
15 
17 TLine *HorizontalLine(TVirtualPad *p, Double_t y)
18 {
19  p->Update();
20 
21  Double_t xMin = p->GetUxmin();
22  Double_t xMax = p->GetUxmax();
23 
24  if (p->GetLogx())
25  {
26  xMin = std::pow(10, xMin);
27  xMax = std::pow(10, xMax);
28  }
29 
30  TLine *line = new TLine(xMin, y, xMax, y);
31  line->SetLineStyle(2);
32  line->SetLineWidth(1);
33  line->SetLineColor(1);
34  return line;
35 }
36 
38 TGraphErrors *
39 FitResolution(const TH2F *h2, const bool normalize_mean = true, const int param = 2)
40 {
41  TProfile *p2 = h2->ProfileX();
42 
43  int n = 0;
44  double x[1000];
45  double ex[1000];
46  double y[1000];
47  double ey[1000];
48 
49  for (int i = 1; i <= h2->GetNbinsX(); i++)
50  {
51  TH1D *h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
52 
53  if (h1->GetSum() < 10)
54  continue;
55 
56  TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
57  p2->GetBinError(i) * 4);
58 
59  TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
60  -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
61 
62  fgaus.SetParameter(1, p2->GetBinContent(i));
63  fgaus.SetParameter(2, p2->GetBinError(i));
64 
65  h1->Fit(&fgaus, "MQ0");
66 
67  x[n] = p2->GetBinCenter(i);
68  ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
69 
70  const double norm = normalize_mean ? fgaus.GetParameter(1) : 1;
71 
72  y[n] = fgaus.GetParameter(param) / norm;
73  ey[n] = fgaus.GetParError(param) / norm;
74 
75  n++;
76  delete h1;
77  }
78 
79  TGraphErrors *ge = new TGraphErrors(n, x, y, 0, ey);
80  ge->SetName(TString(h2->GetName()) + "_FitResolution");
81 
82  ge->SetLineColor(kBlue + 3);
83  ge->SetMarkerColor(kBlue + 3);
84  ge->SetLineWidth(2);
85  ge->SetMarkerStyle(kFullCircle);
86  ge->SetMarkerSize(1);
87  return ge;
88 }
89 
91 TGraphErrors *
92 FitProfile(const TH2 *h2)
93 {
94  TProfile *p2 = h2->ProfileX();
95 
96  int n = 0;
97  double x[1000];
98  double ex[1000];
99  double y[1000];
100  double ey[1000];
101 
102  for (int i = 1; i <= h2->GetNbinsX(); i++)
103  {
104  TH1D *h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
105 
106  if (h1->GetSum() < 10)
107  continue;
108 
109  TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
110  p2->GetBinError(i) * 4);
111 
112  TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
113  -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
114 
115  fgaus.SetParameter(1, p2->GetBinContent(i));
116  fgaus.SetParameter(2, p2->GetBinError(i));
117 
118  h1->Fit(&fgaus, "MQ0");
119 
120  // f2.SetParameters(fgaus.GetParameter(0) / 2, fgaus.GetParameter(1),
121  // fgaus.GetParameter(2), fgaus.GetParameter(0) / 2,
122  // fgaus.GetParameter(2) / 4, 0);
123  //
124  // h1->Fit(&f2, "MQ0");
125 
126  // new TCanvas;
127  // h1->Draw();
128  // fgaus.Draw("same");
129  // break;
130 
131  x[n] = p2->GetBinCenter(i);
132  ex[n] = p2->GetBinWidth(i) / 2;
133  y[n] = fgaus.GetParameter(1);
134  ey[n] = fgaus.GetParameter(2);
135 
136  // p2->SetBinContent(i, fgaus.GetParameter(1));
137  // p2->SetBinError(i, fgaus.GetParameter(2));
138 
139  n++;
140  delete h1;
141  }
142 
143  TGraphErrors *ge = new TGraphErrors(n, x, y, ex, ey);
144  ge->SetName(TString(h2->GetName()) + "_FitProfile");
145  ge->SetLineColor(kBlue + 3);
146  ge->SetMarkerColor(kBlue + 3);
147  ge->SetLineWidth(2);
148  ge->SetMarkerStyle(kFullCircle);
149  ge->SetMarkerSize(1);
150  return ge;
151 }
152 
154 TH1 *GetBinominalRatio(TH1 *h_pass, TH1 *h_n_trial, bool process_zero_bins = false)
155 {
156  assert(h_pass);
157  assert(h_n_trial);
158 
159  assert(h_pass->GetNbinsX() == h_n_trial->GetNbinsX());
160  assert(h_pass->GetNbinsY() == h_n_trial->GetNbinsY());
161  assert(h_pass->GetNbinsZ() == h_n_trial->GetNbinsZ());
162 
163  TH1 *h_ratio = (TH1 *) h_pass->Clone(TString(h_pass->GetName()) + "_Ratio");
164  assert(h_ratio);
165  h_ratio->Divide(h_n_trial); // a rough estimation first, also taking care of the overflow bins and zero bins
166 
167  for (int x = 1; x <= h_n_trial->GetNbinsX(); ++x)
168  for (int y = 1; y <= h_n_trial->GetNbinsY(); ++y)
169  for (int z = 1; z <= h_n_trial->GetNbinsZ(); ++z)
170  {
171  const double n_trial = h_n_trial->GetBinContent(x, y, z);
172 
173  if (n_trial > 0)
174  {
175  const double p = h_pass->GetBinContent(x, y, z) / n_trial;
176 
177  // Wilson score interval
178  h_ratio->SetBinContent(x, y, z, //
179  (p + 1 / (2 * n_trial)) / (1 + 1 / n_trial));
180  h_ratio->SetBinError(x, y,
181  z, //
182  std::sqrt(
183  1. / n_trial * p * (1 - p) + 1. / (4 * n_trial * n_trial)) /
184  (1 + 1 / n_trial));
185  }
186  else if (process_zero_bins)
187  {
188  h_ratio->SetBinContent(x, y, z, 0.5);
189  h_ratio->SetBinError(x, y, z, 0.5);
190  }
191  }
192 
193  return h_ratio;
194 }
195 
196 
198 TLine *VerticalLine(TVirtualPad *p, Double_t x)
199 {
200  p->Update();
201 
202  Double_t yMin = p->GetUymin();
203  Double_t yMax = p->GetUymax();
204 
205  if (p->GetLogy())
206  {
207  yMin = std::pow(10, yMin);
208  yMax = std::pow(10, yMax);
209  }
210 
211  TLine *line = new TLine(x, yMin, x, yMax);
212  line->SetLineStyle(2);
213  line->SetLineWidth(1);
214  line->SetLineColor(1);
215  return line;
216 }
217 double DrawReference(TH1 *hnew, TH1 *href, bool draw_href_error = false, bool do_kstest = true)
218 {
219  hnew->SetLineColor(kBlue + 3);
220  hnew->SetMarkerColor(kBlue + 3);
221  // hnew->SetLineWidth(2);
222  hnew->SetMarkerStyle(kFullCircle);
223  // hnew->SetMarkerSize(1);
224 
225  if (href)
226  {
227  if (draw_href_error)
228  {
229  href->SetLineColor(kGreen + 1);
230  href->SetFillColor(kGreen + 1);
231  href->SetLineStyle(kSolid);
232  href->SetMarkerColor(kGreen + 1);
233  // href->SetLineWidth(2);
234  href->SetMarkerStyle(kDot);
235  href->SetMarkerSize(0);
236  }
237  else
238  {
239  href->SetLineColor(kGreen + 1);
240  href->SetFillColor(kGreen + 1);
241  href->SetLineStyle(0);
242  href->SetMarkerColor(kGreen + 1);
243  href->SetLineWidth(0);
244  href->SetMarkerStyle(kDot);
245  href->SetMarkerSize(0);
246  }
247  }
248 
249  hnew->Draw(); // set scale
250 
251  double ks_test = numeric_limits<double>::signaling_NaN();
252 
253  if (href)
254  {
255  if (draw_href_error)
256  {
257  href->DrawClone("E2 same");
258  href->SetFillStyle(0);
259  href->SetLineWidth(8);
260  href->DrawClone("HIST same ][");
261  }
262  else
263  href->Draw("HIST same");
264  hnew->Draw("same"); // over lay data points
265 
266  if (do_kstest)
267  ks_test = hnew->KolmogorovTest(href, "");
268  }
269 
270  // ---------------------------------
271  // now, make summary header
272  // ---------------------------------
273  if (href)
274  {
275  gPad->SetTopMargin(.14);
276  TLegend *legend = new TLegend(0, .93, 0, 1, hnew->GetTitle(), "NB NDC");
277  legend->Draw();
278  legend = new TLegend(0, .86, .3, .93, NULL, "NB NDC");
279  legend->AddEntry(href, Form("Reference"), "f");
280  legend->Draw();
281  legend = new TLegend(0.3, .86, 1, .93, NULL, "NB NDC");
282 
283  if (do_kstest)
284  {
285  TLegendEntry *le = legend->AddEntry(hnew, Form("New: KS-Test P=%.3f", ks_test), "lpe");
286  if (ks_test >= 1)
287  le->SetTextColor(kBlue + 1);
288  else if (ks_test >= .2)
289  le->SetTextColor(kGreen + 2);
290  else if (ks_test >= .05)
291  le->SetTextColor(kYellow + 1);
292  else
293  le->SetTextColor(kRed + 1);
294  legend->Draw();
295  }
296  else
297  {
298  TLegendEntry *le = legend->AddEntry(hnew, Form("New Result"), "lpe");
299  legend->Draw();
300  }
301  }
302  else
303  {
304  gPad->SetTopMargin(.07);
305  TLegend *legend = new TLegend(0, .93, 0, 1, hnew->GetTitle(), "NB NDC");
306  legend->Draw();
307  }
308 
309 
310  return ks_test;
311 }
312 
315 void DrawReference(TGraph *hnew, TGraph *href, bool draw_href_error = true)
316 {
317  hnew->SetLineColor(kBlue + 3);
318  hnew->SetMarkerColor(kBlue + 3);
319  hnew->SetLineWidth(2);
320  hnew->SetMarkerStyle(kFullCircle);
321  hnew->SetMarkerSize(1);
322 
323  if (href)
324  {
325  if (draw_href_error)
326  {
327  href->SetLineColor(kGreen + 1);
328  href->SetFillColor(kGreen + 1);
329  href->SetFillStyle(0);
330  href->SetLineStyle(kSolid);
331  href->SetMarkerColor(kGreen + 1);
332  href->SetLineWidth(4);
333  href->SetMarkerStyle(kDot);
334  href->SetMarkerSize(0);
335  }
336  else
337  {
338  href->SetLineColor(kGreen + 1);
339  href->SetFillColor(kGreen + 1);
340  href->SetLineStyle(0);
341  href->SetMarkerColor(kGreen + 1);
342  href->SetLineWidth(0);
343  href->SetMarkerStyle(kDot);
344  href->SetMarkerSize(0);
345  }
346  }
347 
348  if (href)
349  {
350  if (draw_href_error)
351  {
352  href->DrawClone("E2");
353  }
354  else
355  href->Draw("HIST same");
356  hnew->Draw("p e"); // over lay data points
357  }
358 
359  // ---------------------------------
360  // now, make summary header
361  // ---------------------------------
362  if (href)
363  {
364  gPad->SetTopMargin(.14);
365  TLegend *legend = new TLegend(0, .93, 0, 1, hnew->GetTitle(), "NB NDC");
366  legend->Draw();
367  legend = new TLegend(0, .86, .3, .93, NULL, "NB NDC");
368  legend->AddEntry(href, Form("Reference"), "f");
369  legend->Draw();
370  legend = new TLegend(0.3, .86, 1, .93, NULL, "NB NDC");
371  TLegendEntry *le = legend->AddEntry(hnew, Form("New"), "lpe");
372  legend->Draw();
373  }
374  else
375  {
376  gPad->SetTopMargin(.07);
377  TLegend *legend = new TLegend(0, .93, 0, 1, hnew->GetTitle(), "NB NDC");
378  legend->Draw();
379  }
380 }
381 
382 
384 void DivideCanvas(TVirtualPad *p, int npads)
385 {
386  if (!p) return;
387  if (!npads) return;
388  const double ratio = double(p->GetWw()) / p->GetWh();
389  Int_t columns = std::max(1, int(std::sqrt(npads * ratio)));
390  Int_t rows = std::max(1, int(npads / columns));
391  while (rows * columns < npads)
392  {
393  columns++;
394  if (rows * columns < npads) rows++;
395  }
396  p->Divide(rows, columns);
397 }