Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrackingQA.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrackingQA.C
1 #include "../CommonTools.h"
2 #include <sPhenixStyle.C>
3 
4 namespace
5 {
6  // Normalization
7  double Nevent_new = 1;
8  double Nevent_ref = 1;
9 
10  void GetNormalization(TFile *qa_file_new, TFile *qa_file_ref, const TString &prefix, const TString &tag)
11  {
12  if (qa_file_new)
13  {
14  TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(prefix + TString("Normalization"), "TH1");
15  assert(h_norm);
16  Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin(tag));
17  }
18 
19  if (qa_file_ref)
20  {
21  TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(prefix + TString("Normalization"), "TH1");
22  assert(h_norm);
23  Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin(tag));
24  }
25  }
26 
27  void Draw(TFile *qa_file_new, TFile *qa_file_ref, const TString &prefix, const TString &tag)
28  {
29  auto h_new = static_cast<TH1 *>(qa_file_new->GetObjectChecked(prefix + tag, "TH1"));
30  assert(h_new);
31  //h_new->Sumw2();
32  h_new->Scale(1. / Nevent_new);
33 
34  TH1 *h_ref = nullptr;
35  if (qa_file_ref)
36  {
37  h_ref = static_cast<TH1 *>(qa_file_ref->GetObjectChecked(prefix + tag, "TH1"));
38  assert(h_ref);
39  //h_ref->Sumw2();
40  h_ref->Scale(1.0 / Nevent_ref);
41  }
42 
43  DrawReference(h_new, h_ref);
44  HorizontalLine(gPad, 1)->Draw();
45  }
46 
47 
48  // square
49 inline double square( const double& x ) { return x*x; }
50 
51 // christal ball function for momentum resolution fit
52 Double_t CBcalc(Double_t *xx, Double_t *par)
53 {
54  // Crystal Ball fit to one state
55  double f;
56  const double x = xx[0];
57 
58  // The four parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
59 
60  const double alpha = par[0];
61  const double n = par[1];
62  const double x_mean = par[2];
63  const double sigma = par[3];
64  const double N = par[4];
65 
66  // dimensionless variable
67  const double t = (x-x_mean)/sigma;
68 
69  // The Crystal Ball function is:
70 
71  if( t > -alpha)
72  {
73  return N * exp( -square(t)/2 );
74  }
75  else
76  {
77  const double A = pow( (n/TMath::Abs(alpha)),n) * exp(-square(alpha)/2.0);
78  const double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
79  return N * A * pow(B - t, -n);
80  }
81 
82 }
83 
84 // christal ball function for Upsilon fits
85 Double_t CBcalc2(Double_t *xx, Double_t *par)
86 {
87  // Crystal Ball fit to one state
88  const double x = xx[0];
89 
90  /* The parameters are:
91  * N the normalization
92  * x_mean, sigma
93  * alpha_left, n_left
94  * alpha_right, n_right
95  */
96  const double N = par[0];
97  const double x_mean = par[1];
98  const double sigma = par[2];
99  const double t = (x-x_mean)/sigma;
100 
101  // left tail
102  const double alpha_left = std::abs(par[3]);
103  const double n_left = par[4];
104 
105  // right tail
106  const double alpha_right = std::abs(par[5]);
107  const double n_right = par[6];
108 
109  // The Crystal Ball function is:
110  if( t < -alpha_left )
111  {
112  const double A = pow( (n_left/TMath::Abs(alpha_left)),n_left) * exp(-square(alpha_left)/2.0);
113  const double B = n_left/std::abs(alpha_left) - std::abs(alpha_left);
114  return N * A * pow(B - t, -n_left);
115  } else if( t > alpha_right ) {
116  const double A = pow( (n_right/TMath::Abs(alpha_right)),n_right) * exp(-square(alpha_right)/2.0);
117  const double B = n_right/std::abs(alpha_right) - std::abs(alpha_right);
118  return N * A * pow(B + t, -n_right);
119  } else {
120  return N * exp( -square(t)/2);
121  }
122 }
123 
124 } // namespace
125 
126 
127 void TrackingQA(std::string reffile, std::string newfile, std::string outfile)
128 {
129  SetsPhenixStyle();
130  TFile *qa_file_new = TFile::Open(newfile.c_str());
131  TFile *qa_file_ref = TFile::Open(reffile.c_str());
132 
133  TFile *outfilef = new TFile(outfile.c_str(), "recreate");
134 
135  {
136  const char *hist_name_prefix = "QAG4SimulationTracking";
137  TString prefix = TString("h_") + hist_name_prefix + TString("_");
138 
139  TCanvas *rawResCan = new TCanvas(TString("QA_Draw_pTResolution") + TString("_") + hist_name_prefix,
140  TString("QA_Draw_pTResolution") + TString("_") + hist_name_prefix,
141  1800,1000);
142  rawResCan->Divide(2,2);
143 
144  TH2F *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new =
145  (TH2F *)qa_file_new->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",
146  "TH2");
147  assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new);
148 
149  h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new->SetDirectory(nullptr);
150 
151  TGraphErrors *newScale = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new, false, 1);
152  TGraphErrors *newRes = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new, false, 2);
153  newScale->SetMarkerStyle(20);
154  newScale->SetMarkerColor(kBlack);
155  newScale->SetMarkerSize(1.0);
156  newRes->SetMarkerStyle(20);
157  newRes->SetMarkerColor(kBlack);
158  newRes->SetMarkerSize(1.0);
159 
160  rawResCan->cd(3);
161  newScale->GetXaxis()->SetTitle("Truth p_{T} [GeV]");
162  newScale->GetYaxis()->SetTitle("#LTp_{T}^{reco}/p_{T}^{true}#GT");
163  newScale->GetYaxis()->SetRangeUser(0.9,1.1);
164  newScale->Draw("ap");
165 
166  rawResCan->cd(4);
167  newRes->GetXaxis()->SetTitle("Truth p_{T} [GeV]");
168  newRes->GetYaxis()->SetTitle("#sigma(p_{T}^{reco}/p_{T}^{true})");
169  newRes->GetYaxis()->SetRangeUser(0,0.08);
170  newRes->Draw("ap");
171 
172  rawResCan->cd(2);
173  gPad->SetRightMargin(0.2);
174  gPad->SetLogx();
175  gPad->SetLogz();
176  h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new->SetTitle("New pT Spectrum");
177  h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new->Draw("colz");
178 
179  TLatex newl;
180  newl.SetTextSize(0.05);
181  newl.SetNDC();
182  newl.SetTextColor(kBlack);
183  newl.DrawLatex(0.45,0.96,"New");
184 
185  if (qa_file_ref) {
186  TH2F *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref =
187  (TH2F *)qa_file_ref->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",
188  "TH2");
189  assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref);
190  h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref->SetDirectory(nullptr);
191 
192  TGraphErrors *refScale = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref, false, 1);
193  TGraphErrors *refRes = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref, false, 2);
194  refScale->SetMarkerStyle(24);
195  refScale->SetMarkerColor(kRed);
196  refScale->SetMarkerSize(1.0);
197  refRes->SetMarkerStyle(24);
198  refRes->SetMarkerColor(kRed);
199  refRes->SetMarkerSize(1.0);
200 
201  TLegend *ptresleg = new TLegend(0.4,0.77,0.5,0.9);
202  ptresleg->AddEntry(refScale,"Reference","P");
203  ptresleg->AddEntry(newScale,"New","P");
204 
205  rawResCan->cd(3);
206  refScale->Draw("psame");
207  ptresleg->Draw("same");
208 
209  rawResCan->cd(4);
210  refRes->Draw("psame");
211  ptresleg->Draw("same");
212 
213  rawResCan->cd(1);
214 
215  gPad->SetRightMargin(0.2);
216  gPad->SetLogx();
217  gPad->SetLogz();
218  h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref->SetTitle("Reference pT Spectrum");
219  h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref->Draw("colz");
220  TLatex refl;
221  refl.SetTextSize(0.05);
222  refl.SetNDC();
223  refl.SetTextColor(kBlack);
224  refl.DrawLatex(0.45,0.96,"Reference");
225 
226  }
227 
228  rawResCan->Draw();
229 
230  outfilef->cd();
231  rawResCan->Write();
232 
233  }
234 
235 
236 
237  {
238  //base histogram from the reco module name
239  const char *hist_name_prefix = "QAG4SimulationTracking";
240  TString prefix = TString("h_") + hist_name_prefix + TString("_");
241 
242  // obtain normalization
243  double Nevent_new = 1;
244  double Nevent_ref = 1;
245 
246  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_TruthMatchingOverview") +
247  TString("_") + hist_name_prefix,
248  TString("QA_Draw_Tracking_TruthMatchingOverview") +
249  TString("_") + hist_name_prefix,
250  1800, 1000);
251  c1->Divide(3, 1);
252  int idx = 1;
253  TPad *p;
254 
255  {
256  static const int nrebin = 5;
257 
258  p = (TPad *)c1->cd(idx++);
259  c1->Update();
260  p->SetLogx();
261  p->SetGridy();
262 
263  TH1 *h_pass =
264  (TH1 *)qa_file_new->GetObjectChecked(prefix + "nReco_pTGen", "TH1");
265  TH1 *h_norm =
266  (TH1 *)qa_file_new->GetObjectChecked(prefix + "nGen_pTGen", "TH1");
267  assert(h_norm);
268  assert(h_pass);
269 
270  h_norm->SetDirectory(nullptr);
271  h_pass->SetDirectory(nullptr);
272 
273  h_norm->Rebin(nrebin);
274  h_pass->Rebin(nrebin);
275 
276  TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);
277 
278  // h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);
279  h_ratio->GetYaxis()->SetTitle("Reco efficiency");
280  h_ratio->GetYaxis()->SetRangeUser(-0, 1.);
281 
282  TH1 *h_ratio_ref = NULL;
283  if (qa_file_ref) {
284  TH1 *h_pass =
285  (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nReco_pTGen", "TH1");
286  TH1 *h_norm =
287  (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nGen_pTGen", "TH1");
288  assert(h_norm);
289  assert(h_pass);
290  h_norm->SetDirectory(nullptr);
291  h_pass->SetDirectory(nullptr);
292  h_norm->Rebin(nrebin);
293  h_pass->Rebin(nrebin);
294  h_ratio_ref = GetBinominalRatio(h_pass, h_norm);
295  }
296 
297  h_ratio->SetTitle(TString(hist_name_prefix) + ": Tracking Efficiency");
298 
299  DrawReference(h_ratio, h_ratio_ref, false);
300 
301  }
302 
303  {
304  static const int nrebin = 4;
305 
306  p = (TPad *)c1->cd(idx++);
307  c1->Update();
308  // p->SetLogx();
309  p->SetGridy();
310 
311  TH1 *h_pass =
312  (TH1 *)qa_file_new->GetObjectChecked(prefix + "nReco_etaGen", "TH1");
313  TH1 *h_norm =
314  (TH1 *)qa_file_new->GetObjectChecked(prefix + "nGen_etaGen", "TH1");
315  assert(h_norm);
316  assert(h_pass);
317 
318  h_norm->SetDirectory(nullptr);
319  h_pass->SetDirectory(nullptr);
320  h_norm->Rebin(nrebin);
321  h_pass->Rebin(nrebin);
322 
323  TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);
324 
325  h_ratio->GetXaxis()->SetRangeUser(-1.1, 1.1);
326  h_ratio->GetYaxis()->SetTitle("Reco efficiency");
327  h_ratio->GetYaxis()->SetRangeUser(-0, 1.);
328 
329  TH1 *h_ratio_ref = NULL;
330  if (qa_file_ref) {
331  TH1 *h_pass =
332  (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nReco_etaGen", "TH1");
333  TH1 *h_norm =
334  (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nGen_etaGen", "TH1");
335  assert(h_norm);
336  assert(h_pass);
337  h_norm->SetDirectory(nullptr);
338  h_pass->SetDirectory(nullptr);
339  h_norm->Rebin(nrebin);
340  h_pass->Rebin(nrebin);
341  h_ratio_ref = GetBinominalRatio(h_pass, h_norm);
342  }
343 
344  h_ratio->SetTitle(TString(hist_name_prefix) + ": Tracking Efficiency");
345 
346  DrawReference(h_ratio, h_ratio_ref, false);
347  }
348 
349  {
350  p = (TPad *)c1->cd(idx++);
351  c1->Update();
352  // p->SetLogx();
353  TH1 *frame = p->DrawFrame(0, .9, 50, 1.1,
354  "Mean and sigma, p_{T,reco}/p_{T,truth};Truth p_{T} [GeV/c];<p_{T,reco}/p_{T,truth}> #pm #sigma(p_{T,reco}/p_{T,truth})");
355  //gPad->SetLeftMargin(.2);
356  gPad->SetTopMargin(-1);
357  frame->GetYaxis()->SetTitleOffset(1.7);
358  //TLine *l = new TLine(0, 1, 50, 1);
359  //l->SetLineColor(kGray);
360  //l->Draw();
361  HorizontalLine( gPad, 1 )->Draw();
362 
363  TH2 *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen =
364  (TH2 *)qa_file_new->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",
365  "TH2");
366  assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);
367 
368  h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->SetDirectory(nullptr);
369  h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Rebin2D(16, 1);
370  TGraphErrors *ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen =
371  FitProfile(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);
372  ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Draw("pe");
373  ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen->SetTitle(
374  "Mean and sigma, p_{T,reco}/p_{T,truth}");
375 
376  TGraphErrors *h_ratio_ref = NULL;
377  if (qa_file_ref) {
378  TH2 *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen =
379  (TH2 *)qa_file_ref->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",
380  "TH2");
381  assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);
382 
383  h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->SetDirectory(nullptr);
384  h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Rebin2D(16, 1);
385 
386  h_ratio_ref = FitProfile(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);
387  ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Draw("pe");
388  }
389 
390  DrawReference(ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen, h_ratio_ref,
391  true);
392 
393  }
394 
395  //SaveCanvas(c1,
396  // TString(qa_file_name_new) + TString("_") + TString(c1->GetName()),
397  // true);
398 
399  c1->Draw();
400  outfilef->cd();
401  c1->Write();
402 }
403 
404 
405 {
406  const char *hist_name_prefix = "QAG4SimulationTracking";
407  TString prefix = TString("h_") + hist_name_prefix + TString("_");
408 
409  // obtain normalization
410  double Nevent_new = 1;
411  double Nevent_ref = 1;
412 
413 
414  TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(
415  prefix + TString("pTRecoGenRatio_pTGen"), "TH2");
416  assert(h_new);
417 
418  // h_new->Rebin(1, 2);
419  //h_new->Sumw2();
420  // h_new->Scale(1. / Nevent_new);
421 
422  TH2 *h_ref = NULL;
423  if (qa_file_ref)
424  {
425  h_ref = (TH2 *) qa_file_ref->GetObjectChecked(
426  prefix + TString("pTRecoGenRatio_pTGen"), "TH2");
427  assert(h_ref);
428 
429  // h_ref->Rebin(1, 2);
430  //h_ref->Sumw2();
431  h_ref->Scale(Nevent_new / Nevent_ref);
432  }
433 
434  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_pTRatio") + TString("_") + hist_name_prefix,
435  TString("QA_Draw_Tracking_pTRatio") + TString("_") + hist_name_prefix,
436  1800, 1000);
437  c1->Divide(4, 2);
438  int idx = 1;
439  TPad *p;
440 
441  vector<pair<double, double>> gpt_ranges{
442  {0, 1},
443  {1, 5},
444  {5, 10},
445  {10, 20},
446  {20, 30},
447  {30, 40},
448  {40, 45},
449  {45, 50}};
450  TF1 *f1 = nullptr;
451  TF1 *fit = nullptr;
452  Double_t sigma = 0;
453  Double_t sigma_unc = 0;
454  char resstr[500];
455  TLatex *res = nullptr;
456  for (auto pt_range : gpt_ranges)
457  {
458  //cout << __PRETTY_FUNCTION__ << " process " << pt_range.first << " - " << pt_range.second << " GeV/c";
459 
460  p = (TPad *) c1->cd(idx++);
461  c1->Update();
462  p->SetLogy();
463 
464  const double epsilon = 1e-6;
465  const int bin_start = h_new->GetXaxis()->FindBin(pt_range.first + epsilon);
466  const int bin_end = h_new->GetXaxis()->FindBin(pt_range.second - epsilon);
467 
468  TH1 *h_proj_new = h_new->ProjectionY(
469  TString::Format(
470  "%s_New_ProjX_%d_%d",
471  h_new->GetName(), bin_start, bin_end),
472  bin_start, bin_end);
473 
474  h_proj_new->GetXaxis()->SetRangeUser(.7, 1.3);
475  h_proj_new->SetTitle(TString(hist_name_prefix) + TString::Format(
476  ": %.1f - %.1f GeV/c", pt_range.first, pt_range.second));
477  h_proj_new->GetXaxis()->SetTitle(TString::Format(
478  "Reco p_{T}/Truth p_{T}"));
479 
480  f1 = new TF1("f1", "gaus", -.85, 1.15);
481  h_proj_new->Fit(f1, "mq");
482  sigma = f1->GetParameter(2);
483  sigma_unc = f1->GetParError(2);
484 
485  TH1 *h_proj_ref = nullptr;
486  if (h_ref)
487  h_proj_ref =
488  h_ref->ProjectionY(
489  TString::Format(
490  "%s_Ref_ProjX_%d_%d",
491  h_new->GetName(), bin_start, bin_end),
492  bin_start, bin_end);
493 
494  DrawReference(h_proj_new, h_proj_ref);
495  sprintf(resstr, "#sigma = %.5f #pm %.5f", sigma, sigma_unc);
496  res = new TLatex(0.325, 0.825, resstr);
497  res->SetNDC();
498  res->SetTextSize(0.05);
499  res->SetTextAlign(13);
500  res->Draw();
501  }
502 
503  // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);
504  c1->Draw();
505  outfilef->cd();
506  c1->Write();
507 }
508 
509 
510  {
511 const char *hist_name_prefix = "QAG4SimulationTracking";
512  TString prefix = TString("h_") + hist_name_prefix + TString("_");
513 
514 
515  // obtain normalization
516  double Nevent_new = 1;
517  double Nevent_ref = 1;
518 
519  if (qa_file_new)
520  {
521  //cout << "Open new QA file " << qa_file_new->GetName() << endl;
522 
523  TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(
524  prefix + TString("Normalization"), "TH1");
525  assert(h_norm);
526 
527  Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
528  }
529  if (qa_file_ref)
530  {
531  // cout << "Open ref QA file " << qa_file_ref->GetName() << endl;
532  TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(
533  prefix + TString("Normalization"), "TH1");
534  assert(h_norm);
535 
536  Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
537  }
538 
539 
540  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_RecoTruthMatching") +
541  TString("_") + hist_name_prefix,
542  TString("QA_Draw_Tracking_RecoTruthMatching") +
543  TString("_") + hist_name_prefix,
544  1800, 1000);
545  c1->Divide(2, 1);
546  int idx = 1;
547  TPad *p;
548 
549  {
550  static const int nrebin = 5;
551 
552  p = (TPad *) c1->cd(idx++);
553  c1->Update();
554  p->SetLogx();
555  p->SetGridy();
556 
557  TH1 *h_pass =
558  (TH1 *) qa_file_new->GetObjectChecked(prefix + "nGen_pTReco", "TH1");
559  TH1 *h_norm =
560  (TH1 *) qa_file_new->GetObjectChecked(prefix + "nReco_pTReco", "TH1");
561  assert(h_norm);
562  assert(h_pass);
563 
564  h_norm->SetDirectory(nullptr);
565  h_pass->SetDirectory(nullptr);
566 
567  h_norm->Rebin(nrebin);
568  h_pass->Rebin(nrebin);
569 
570  TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);
571 
572  // h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);
573  h_ratio->GetYaxis()->SetTitle("Tracking Purity");
574  h_ratio->GetYaxis()->SetRangeUser(-0, 1.1);
575 
576  TH1 *h_ratio_ref = NULL;
577  if (qa_file_ref)
578  {
579  TH1 *h_pass =
580  (TH1 *) qa_file_ref->GetObjectChecked(prefix + "nGen_pTReco", "TH1");
581  TH1 *h_norm =
582  (TH1 *) qa_file_ref->GetObjectChecked(prefix + "nReco_pTReco", "TH1");
583  assert(h_norm);
584  assert(h_pass);
585  h_norm->SetDirectory(nullptr);
586  h_pass->SetDirectory(nullptr);
587  h_norm->Rebin(nrebin);
588  h_pass->Rebin(nrebin);
589  h_ratio_ref = GetBinominalRatio(h_pass, h_norm);
590  }
591 
592  h_ratio->SetTitle("Tracking Purity (matched truth-reco pairs)");
593 
594  DrawReference(h_ratio, h_ratio_ref, false);
595  }
596 
597  {
598  p = (TPad *) c1->cd(idx++);
599  c1->Update();
600  // p->SetLogx();
601  TH1 *frame = p->DrawFrame(0, .9, 50, 1.1,
602  "Mean and sigma p_{Tmatched}/p_{Treco};Reco p_{T} [GeV/c];<p_{T,matched}/p_{T,reco}> #pm #sigma(p_{T,matched}/p_{T,reco})");
603  // gPad->SetLeftMargin(.2);
604  gPad->SetTopMargin(-1);
605  frame->GetYaxis()->SetTitleOffset(1.7);
606  // TLine *l = new TLine(0, 1, 50, 1);
607  // l->SetLineColor(kGray);
608  // l->Draw();
609  HorizontalLine(gPad, 1)->Draw();
610 
611  TH2 *h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =
612  (TH2 *) qa_file_new->GetObjectChecked(
613  prefix + "pTRecoTruthMatchedRatio_pTReco", "TH2");
614  assert(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);
615 
616  h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetDirectory(nullptr);
617  h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Rebin2D(16, 1);
618 
619  TGraphErrors *ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =
620  FitProfile(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);
621  ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Draw("pe");
622  ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetTitle(
623  "Mean and sigma p_{Tmatched}/p_{Treco}");
624 
625  TGraphErrors *h_ratio_ref = NULL;
626  if (qa_file_ref)
627  {
628  TH2 *h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =
629  (TH2 *) qa_file_ref->GetObjectChecked(
630  prefix + "pTRecoTruthMatchedRatio_pTReco", "TH2");
631  assert(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);
632 
633  h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetDirectory(nullptr);
634  h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Rebin2D(16, 1);
635 
636  h_ratio_ref =
637  FitProfile(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);
638  ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Draw("pe");
639  }
640 
641  DrawReference(ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco,
642  h_ratio_ref, true);
643  }
644 
645  c1->Draw();
646  outfilef->cd();
647  c1->Write();
648 }
649 
650 
651 {
652  const char *hist_name_prefix = "QAG4SimulationTracking";
653  TString prefix = TString("h_") + hist_name_prefix + TString("_");
654 
655  // obtain normalization
656  double Nevent_new = 1;
657  double Nevent_ref = 1;
658 
659 
660 
661  if (qa_file_new)
662  {
663  TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(
664  prefix + TString("Normalization"), "TH1");
665  assert(h_norm);
666 
667  Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
668  }
669  if (qa_file_ref)
670  {
671  TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(
672  prefix + TString("Normalization"), "TH1");
673  assert(h_norm);
674 
675  Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
676  }
677 
678  //MVTX, INTT, TPC
679  vector<TString> detectors{"MVTX", "INTT", "TPC"};
680  vector<int> eff_ncluster_cuts{2, 2, 40};
681  vector<double> ncluster_spectrum_pt_cuts{2, 2, 2};
682  vector<TH2 *> h_pass_detectors(3, nullptr);
683  static const int nrebin = 5;
684 
685  h_pass_detectors[0] = (TH2 *) qa_file_new->GetObjectChecked(
686  prefix + "nMVTX_nReco_pTGen", "TH1") ;
687  h_pass_detectors[1] = (TH2 *) qa_file_new->GetObjectChecked(
688  prefix + "nINTT_nReco_pTGen", "TH1") ;
689  h_pass_detectors[2] = (TH2 *) qa_file_new->GetObjectChecked(
690  prefix + "nTPC_nReco_pTGen", "TH1") ;
691 
692  TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(
693  prefix + "nGen_pTGen", "TH1") ;
694  assert(h_norm);
695  h_norm->SetDirectory(nullptr);
696  h_norm->Rebin(nrebin);
697 
698  vector<TH2 *> h_pass_detectors_ref(3, nullptr);
699  TH1 *h_norm_ref = nullptr;
700  if (qa_file_ref)
701  {
702  h_pass_detectors_ref[0] = (TH2 *) qa_file_ref->GetObjectChecked(
703  prefix + "nMVTX_nReco_pTGen", "TH1") ;
704  h_pass_detectors_ref[1] = (TH2 *) qa_file_ref->GetObjectChecked(
705  prefix + "nINTT_nReco_pTGen", "TH1") ;
706  h_pass_detectors_ref[2] = (TH2 *) qa_file_ref->GetObjectChecked(
707  prefix + "nTPC_nReco_pTGen", "TH1") ;
708 
709  h_norm_ref = (TH1 *) qa_file_ref->GetObjectChecked(
710  prefix + "nGen_pTGen", "TH1") ;
711  h_norm_ref->SetDirectory(nullptr);
712  h_norm_ref->Rebin(nrebin);
713 
714  }
715 
716  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_TruthMatching_NumOfClusters") + TString("_") + hist_name_prefix,
717  TString("QA_Draw_Tracking_TruthMatching_NumOfClusters") + TString("_") + hist_name_prefix,
718  1800, 1000);
719  c1->Divide(3, 2);
720  TPad *p;
721 
722  for (int i = 0; i < 3; ++i)
723  {
724  TString detector = detectors[i];
725  TH2 *h_pass_detector = h_pass_detectors[i];
726  TH2 *h_pass_detector_ref = h_pass_detectors_ref[i];
727  assert(h_pass_detector);
728 
729  {
730  p = (TPad *) c1->cd(i + 1);
731  c1->Update();
732  p->SetLogy();
733 
734  const int bin_start = h_pass_detector->GetXaxis()->FindBin(ncluster_spectrum_pt_cuts[i]);
735 
736  TH1 *h_pass_detector_ncluster = h_pass_detector->ProjectionY(
737  TString(h_pass_detector->GetName()) + "_nCluster_new",
738  bin_start);
739  TH1 *h_pass_detector_ncluster_ref = nullptr;
740  if (h_pass_detector_ref)
741  {
742  h_pass_detector_ncluster_ref = h_pass_detector_ref->ProjectionY(
743  TString(h_pass_detector_ref->GetName()) + "_nCluster_ref",
744  bin_start);
745  }
746 
747  h_pass_detector_ncluster->SetTitle(TString(hist_name_prefix) + ": " + detector + Form(" n_{Cluster} | p_{T} #geq %.1fGeV/c", ncluster_spectrum_pt_cuts[i]));
748  h_pass_detector_ncluster->SetYTitle("# of reconstructed track");
749  DrawReference(h_pass_detector_ncluster, h_pass_detector_ncluster_ref, false);
750  }
751 
752  {
753  p = (TPad *) c1->cd(i + 3 + 1);
754  c1->Update();
755  p->SetLogx();
756  p->SetGridy();
757 
758  const int bin_start = h_pass_detector->GetYaxis()->FindBin(eff_ncluster_cuts[i]);
759  TH1 *h_pass = h_pass_detector->ProjectionX(
760  TString(h_pass_detector->GetName()) + "_nReco_new",
761  bin_start);
762 
763  assert(h_pass);
764  h_pass->SetDirectory(nullptr);
765  h_pass->Rebin(nrebin);
766 
767  TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);
768  h_ratio->GetYaxis()->SetTitle("Reco efficiency | " + detector + Form(" n_{Cluster} #geq %d", eff_ncluster_cuts[i]));
769  h_ratio->GetYaxis()->SetRangeUser(-0, 1.);
770  //
771  TH1 *h_ratio_ref = NULL;
772  if (h_pass_detector_ref)
773  {
774  TH1 *h_pass = h_pass_detector_ref->ProjectionX(
775  TString(h_pass_detector->GetName()) + "_nReco_ref",
776  bin_start);
777 
778  assert(h_pass);
779  h_pass->SetDirectory(nullptr);
780  h_pass->Rebin(nrebin);
781 
782  h_ratio_ref = GetBinominalRatio(h_pass, h_norm_ref);
783  }
784  //
785  h_ratio->SetTitle("Tracking efficiency | " + detector + Form(" n_{Cluster} #geq %d", eff_ncluster_cuts[i]));
786  DrawReference(h_ratio, h_ratio_ref, false);
787  }
788  }
789 
790  // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);
791  c1->Draw();
792  outfilef->cd();
793  c1->Write();
794 }
795 
796 
797 
798  {
799  const char *hist_name_prefix = "QAG4SimulationTracking";
800  TString prefix = TString("h_") + hist_name_prefix + TString("_");
801 
802  auto c1 = new TCanvas(TString("QA_Draw_Tracking_nClus_Layer") + TString("_") + hist_name_prefix,
803  TString("QA_Draw_Tracking_nClus_Layer") + TString("_") + hist_name_prefix,
804  1800, 1000);
805 
806  c1->Divide(2, 1);
807  c1->cd(1);
808  GetNormalization(qa_file_new, qa_file_ref, prefix, "Truth Track");
809  Draw(qa_file_new, qa_file_ref, prefix, "nClus_layerGen");
810 
811  c1->cd(2);
812  GetNormalization(qa_file_new, qa_file_ref, prefix, "Reco Track");
813  Draw(qa_file_new, qa_file_ref, prefix, "nClus_layer");
814 
815  // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);
816  c1->Draw();
817  outfilef->cd();
818  c1->Write();
819 }
820 
821 
822 {
823  const char *hist_name_prefix = "QAG4SimulationUpsilon";
824  TString prefix = TString("h_") + hist_name_prefix + TString("_");
825 
826  // obtain normalization
827  double Nevent_new = 1;
828  double Nevent_ref = 1;
829 
830  if ( qa_file_new->GetObjectChecked(
831  prefix + TString("pTRecoGenRatio_pTGen"), "TH2")
832  == nullptr )
833  {
834  cout <<"QAG4SimulationUpsilon is not enabled. Skip...."<<endl;
835  }
836  else
837  {
838 
839  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_UpsilonOverview") + TString("_") + hist_name_prefix,
840  TString("QA_Draw_Tracking_UpsilonOverview") + TString("_") + hist_name_prefix,
841  1800, 1000);
842  c1->Divide(2, 1);
843  int idx = 1;
844  TPad *p;
845 
846  {
847  p = (TPad *) c1->cd(idx++);
848  c1->Update();
849  p->SetLogy();
850 
851  TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(
852  prefix + TString("pTRecoGenRatio_pTGen"), "TH2");
853  assert(h_new);
854 
855  // h_new->Rebin(1, 2);
856  //h_new->Sumw2();
857  // h_new->Scale(1. / Nevent_new);
858 
859  TH2 *h_ref = NULL;
860  if (qa_file_ref)
861  {
862  h_ref = (TH2 *) qa_file_ref->GetObjectChecked(
863  prefix + TString("pTRecoGenRatio_pTGen"), "TH2");
864  assert(h_ref);
865 
866  // h_ref->Rebin(1, 2);
867  //h_ref->Sumw2();
868  h_ref->Scale(Nevent_new / Nevent_ref);
869  }
870 
871  TH1 *h_proj_new = h_new->ProjectionY(
872  TString::Format(
873  "%s_New_ProjX",
874  h_new->GetName()));
875 
876  h_proj_new->GetXaxis()->SetRangeUser(0, 1.3);
877  h_proj_new->SetTitle(TString(hist_name_prefix) + TString::Format(
878  ": Electron lineshape"));
879  h_proj_new->GetXaxis()->SetTitle(TString::Format(
880  "Reco p_{T}/Truth p_{T}"));
881 
882  TF1 *f_eLineshape = new TF1("f_eLineshape", CBcalc, 7, 11, 5);
883  f_eLineshape->SetParameter(0, 1.0);
884  f_eLineshape->SetParameter(1, 1.0);
885  f_eLineshape->SetParameter(2, 0.95);
886  f_eLineshape->SetParameter(3, 0.08);
887  f_eLineshape->SetParameter(4, 20.0);
888 
889  f_eLineshape->SetParNames("alpha","n","m","sigma","N");
890  f_eLineshape->SetLineColor(kRed);
891  f_eLineshape->SetLineWidth(3);
892  f_eLineshape->SetLineStyle(kSolid);
893  f_eLineshape->SetNpx(1000);
894 
895  h_proj_new->Fit(f_eLineshape);
896 
897  TH1 *h_proj_ref = nullptr;
898  if (h_ref)
899  {
900  h_proj_ref =
901  h_ref->ProjectionY(
902  TString::Format(
903  "%s_Ref_ProjX",
904  h_new->GetName()));
905  }
906  TF1 *f_eLineshape_ref = new TF1("f_eLineshape_ref", CBcalc, 7, 11, 5);
907  f_eLineshape_ref->SetParameter(0, 1.0);
908  f_eLineshape_ref->SetParameter(1, 1.0);
909  f_eLineshape_ref->SetParameter(2, 0.95);
910  f_eLineshape_ref->SetParameter(3, 0.08);
911  f_eLineshape_ref->SetParameter(4, 20.0);
912 
913  f_eLineshape_ref->SetParNames("alpha","n","m","sigma","N");
914  f_eLineshape_ref->SetLineColor(kRed);
915  f_eLineshape_ref->SetLineWidth(3);
916  f_eLineshape_ref->SetLineStyle(kSolid);
917 
918  h_proj_ref->Fit(f_eLineshape_ref);
919 
920 
921  DrawReference(h_proj_new, h_proj_ref);
922  f_eLineshape->Draw("same");
923 
924  char resstr_1[500];
925  sprintf(resstr_1,"#sigma_{dp/p} = %.2f #pm %.2f %%", f_eLineshape->GetParameter(3)*100, f_eLineshape->GetParError(3)*100);
926  TLatex *res_1 = new TLatex(0.2,0.75,resstr_1);
927  res_1->SetNDC();
928  res_1->SetTextSize(0.05);
929  res_1->SetTextAlign(13);
930  res_1->Draw();
931 
932  char resstr_2[500];
933  sprintf(resstr_2,"#sigma_{dp/p,ref} = %.2f #pm %.2f %%", f_eLineshape_ref->GetParameter(3)*100, f_eLineshape_ref->GetParError(3)*100);
934  TLatex *res_2 = new TLatex(0.2,0.7,resstr_2);
935  res_2->SetNDC();
936  res_2->SetTextSize(0.05);
937  res_2->SetTextAlign(13);
938  res_2->Draw();
939  }
940 
941  {
942  p = (TPad *) c1->cd(idx++);
943  c1->Update();
944  // p->SetLogy();
945 
946  TH1 *h_new = (TH1 *) qa_file_new->GetObjectChecked(
947  prefix + TString("nReco_Pair_InvMassReco"), "TH1");
948  assert(h_new);
949 
950  // h_new->Rebin(2);
951  // h_new->Sumw2();
952  // h_new->Scale(1. / Nevent_new);
953 
954  TF1 *f1S = new TF1("f1S", CBcalc2, 7, 11, 7);
955  f1S->SetParameter(0, 50.0);
956  f1S->SetParameter(1, 9.46);
957  f1S->SetParameter(2, 0.08);
958  f1S->SetParameter(3, 1.0);
959  f1S->SetParameter(4, 3.0);
960  f1S->SetParameter(5, 1.0);
961  f1S->SetParLimits(3, 0.120, 10);
962  f1S->SetParLimits(4, 1.05, 10);
963  f1S->SetParameter(6, 3.0);
964  f1S->SetParLimits(5, 0.120, 10);
965  f1S->SetParLimits(6, 1.05, 10);
966 
967  f1S->SetParNames("N", "m", "#sigma", "#alpha_{left}","n_{left}","#alpha_{right}","#sigma_{right}");
968  f1S->SetLineColor(kRed);
969  f1S->SetLineWidth(3);
970  f1S->SetLineStyle(kSolid);
971  f1S->SetNpx(1000);
972 
973  h_new->Fit(f1S);
974 
975  TH1 *h_ref = NULL;
976  if (qa_file_ref)
977  {
978  h_ref = (TH1 *) qa_file_ref->GetObjectChecked(
979  prefix + TString("nReco_Pair_InvMassReco"), "TH1");
980  assert(h_ref);
981 
982  // h_ref->Rebin(2);
983  // h_ref->Sumw2();
984  // h_ref->Scale(Nevent_new / Nevent_ref);
985  }
986 
987  h_new->SetTitle(TString(hist_name_prefix) + TString::Format(
988  ": #Upsilon #rightarrow e^{+}e^{-} lineshape"));
989  h_new->GetXaxis()->SetRangeUser(7, 10);
990 
991  TF1 *f1S_ref = new TF1("f1S_ref", CBcalc2, 7, 11, 7);
992  f1S_ref->SetParameter(0, 50.0);
993  f1S_ref->SetParameter(1, 9.46);
994  f1S_ref->SetParameter(2, 0.08);
995  f1S_ref->SetParameter(3, 1.0);
996  f1S_ref->SetParameter(4, 3.0);
997  f1S_ref->SetParameter(5, 1.0);
998  f1S_ref->SetParLimits(3, 0.120, 10);
999  f1S_ref->SetParLimits(4, 1.05, 10);
1000  f1S_ref->SetParameter(6, 3.0);
1001  f1S_ref->SetParLimits(5, 0.120, 10);
1002  f1S_ref->SetParLimits(6, 1.05, 10);
1003 
1004  f1S_ref->SetParNames("N", "m", "#sigma", "#alpha_{left}","n_{left}","#alpha_{right}","#sigma_{right}");
1005  f1S_ref->SetLineColor(kRed);
1006  f1S_ref->SetLineWidth(3);
1007  f1S_ref->SetLineStyle(kSolid);
1008 
1009  h_ref->Fit(f1S_ref);
1010 
1011  DrawReference(h_new, h_ref, false);
1012  f1S->Draw("same");
1013 
1014  // cout << "f1S pars " << f1S->GetParameter(3) << " " << f1S->GetParError(3) << endl;
1015 
1016  char resstr_3[500];
1017  sprintf(resstr_3,"#sigma_{1S} = %.1f #pm %.1f MeV", f1S->GetParameter(2)*1000, f1S->GetParError(2)*1000);
1018  TLatex *res_3 = new TLatex(0.2,0.75,resstr_3);
1019  res_3->SetNDC();
1020  res_3->SetTextSize(0.05);
1021  res_3->SetTextAlign(13);
1022  res_3->Draw();
1023 
1024  char resstr_4[500];
1025  sprintf(resstr_4,"#sigma_{1S,ref} = %.1f #pm %.1f MeV", f1S_ref->GetParameter(2)*1000, f1S_ref->GetParError(2)*1000);
1026  TLatex *res_4 = new TLatex(0.2,0.7,resstr_4);
1027  res_4->SetNDC();
1028  res_4->SetTextSize(0.05);
1029  res_4->SetTextAlign(13);
1030  res_4->Draw();
1031  }
1032 
1033  // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);
1034 
1035  c1 -> Draw();
1036  outfilef->cd();
1037  c1->Write();
1038  }// if checks
1039  }
1040 
1041  outfilef->Close();
1042 
1043 }