BenchmarkingTools/PhysicsPerformanceTools/TrackingQA.C

0001 #include "../CommonTools.h"
0002 #include <sPhenixStyle.C>
0003
0004 namespace
0005 {
0006   // Normalization
0007   double Nevent_new = 1;
0008   double Nevent_ref = 1;
0009
0010   void GetNormalization(TFile *qa_file_new, TFile *qa_file_ref, const TString &prefix, const TString &tag)
0011   {
0012     if (qa_file_new)
0013     {
0014       TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(prefix + TString("Normalization"), "TH1");
0015       assert(h_norm);
0016       Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin(tag));
0017     }
0018
0019     if (qa_file_ref)
0020     {
0021       TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(prefix + TString("Normalization"), "TH1");
0022       assert(h_norm);
0023       Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin(tag));
0024     }
0025   }
0026
0027   void Draw(TFile *qa_file_new, TFile *qa_file_ref, const TString &prefix, const TString &tag)
0028   {
0029     auto h_new = static_cast<TH1 *>(qa_file_new->GetObjectChecked(prefix + tag, "TH1"));
0030     assert(h_new);
0031     //h_new->Sumw2();
0032     h_new->Scale(1. / Nevent_new);
0033
0034     TH1 *h_ref = nullptr;
0035     if (qa_file_ref)
0036     {
0037       h_ref = static_cast<TH1 *>(qa_file_ref->GetObjectChecked(prefix + tag, "TH1"));
0038       assert(h_ref);
0039       //h_ref->Sumw2();
0040       h_ref->Scale(1.0 / Nevent_ref);
0041     }
0042
0043     DrawReference(h_new, h_ref);
0044     HorizontalLine(gPad, 1)->Draw();
0045   }
0046
0047
0048   // square
0049 inline double square( const double& x ) { return x*x; }
0050
0051 // christal ball function for momentum resolution fit
0052 Double_t CBcalc(Double_t *xx, Double_t *par)
0053 {
0054   // Crystal Ball fit to one state
0055   double f;
0056   const double x = xx[0];
0057
0058   // The four parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
0059
0060   const double alpha = par[0];
0061   const double n = par[1];
0062   const double x_mean = par[2];
0063   const double sigma = par[3];
0064   const double N = par[4];
0065
0066   // dimensionless variable
0067   const double t = (x-x_mean)/sigma;
0068
0069   // The Crystal Ball function is:
0070
0071   if( t > -alpha)
0072   {
0073     return N * exp( -square(t)/2 );
0074   }
0075   else
0076   {
0077     const double A = pow( (n/TMath::Abs(alpha)),n) * exp(-square(alpha)/2.0);
0078     const double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
0079     return N * A * pow(B - t, -n);
0080   }
0081
0082 }
0083
0084 // christal ball function for Upsilon fits
0085 Double_t CBcalc2(Double_t *xx, Double_t *par)
0086 {
0087   // Crystal Ball fit to one state
0088   const double x = xx[0];
0089
0090   /* The parameters are:
0091    * N the normalization
0092    * x_mean, sigma
0093    * alpha_left, n_left
0094    * alpha_right, n_right
0095    */
0096   const double N = par[0];
0097   const double x_mean = par[1];
0098   const double sigma = par[2];
0099   const double t = (x-x_mean)/sigma;
0100
0101   // left tail
0102   const double alpha_left = std::abs(par[3]);
0103   const double n_left = par[4];
0104
0105   // right tail
0106   const double alpha_right = std::abs(par[5]);
0107   const double n_right = par[6];
0108
0109   // The Crystal Ball function is:
0110   if( t < -alpha_left )
0111   {
0112     const double A = pow( (n_left/TMath::Abs(alpha_left)),n_left) * exp(-square(alpha_left)/2.0);
0113     const double B = n_left/std::abs(alpha_left) - std::abs(alpha_left);
0114     return N * A * pow(B - t, -n_left);
0115   } else if( t > alpha_right ) {
0116      const double A = pow( (n_right/TMath::Abs(alpha_right)),n_right) * exp(-square(alpha_right)/2.0);
0117      const double B = n_right/std::abs(alpha_right) - std::abs(alpha_right);
0118      return N * A * pow(B + t, -n_right);
0119   } else {
0120       return N * exp( -square(t)/2);
0121   }
0122 }
0123
0124 }  // namespace
0125
0126
0127 void TrackingQA(std::string reffile, std::string newfile, std::string outfile)
0128 {
0129   SetsPhenixStyle();
0130   TFile *qa_file_new = TFile::Open(newfile.c_str());
0131   TFile *qa_file_ref = TFile::Open(reffile.c_str());
0132
0133   TFile *outfilef = new TFile(outfile.c_str(), "recreate");
0134
0135   {
0136     const char *hist_name_prefix = "QAG4SimulationTracking";
0137     TString prefix = TString("h_") + hist_name_prefix + TString("_");
0138
0139     TCanvas *rawResCan = new TCanvas(TString("QA_Draw_pTResolution") + TString("_") + hist_name_prefix,
0140                                      TString("QA_Draw_pTResolution") + TString("_") + hist_name_prefix,
0141                                      1800,1000);
0142     rawResCan->Divide(2,2);
0143
0144     TH2F *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new =
0145         (TH2F *)qa_file_new->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",
0146                                              "TH2");
0147     assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new);
0148
0149     h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new->SetDirectory(nullptr);
0150
0151     TGraphErrors *newScale = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new, false, 1);
0152     TGraphErrors *newRes = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new, false, 2);
0153     newScale->SetMarkerStyle(20);
0154     newScale->SetMarkerColor(kBlack);
0155     newScale->SetMarkerSize(1.0);
0156     newRes->SetMarkerStyle(20);
0157     newRes->SetMarkerColor(kBlack);
0158     newRes->SetMarkerSize(1.0);
0159
0160     rawResCan->cd(3);
0161     newScale->GetXaxis()->SetTitle("Truth p_{T} [GeV]");
0162     newScale->GetYaxis()->SetTitle("#LTp_{T}^{reco}/p_{T}^{true}#GT");
0163     newScale->GetYaxis()->SetRangeUser(0.9,1.1);
0164     newScale->Draw("ap");
0165
0166     rawResCan->cd(4);
0167     newRes->GetXaxis()->SetTitle("Truth p_{T} [GeV]");
0168     newRes->GetYaxis()->SetTitle("#sigma(p_{T}^{reco}/p_{T}^{true})");
0169     newRes->GetYaxis()->SetRangeUser(0,0.08);
0170     newRes->Draw("ap");
0171
0172     rawResCan->cd(2);
0173     gPad->SetRightMargin(0.2);
0174     gPad->SetLogx();
0175     gPad->SetLogz();
0176     h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new->SetTitle("New pT Spectrum");
0177     h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new->Draw("colz");
0178
0179     TLatex newl;
0180     newl.SetTextSize(0.05);
0181     newl.SetNDC();
0182     newl.SetTextColor(kBlack);
0183     newl.DrawLatex(0.45,0.96,"New");
0184
0185     if (qa_file_ref) {
0186       TH2F *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref =
0187           (TH2F *)qa_file_ref->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",
0188                                                "TH2");
0189       assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref);
0190       h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref->SetDirectory(nullptr);
0191
0192       TGraphErrors *refScale = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref, false, 1);
0193       TGraphErrors *refRes = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref, false, 2);
0194       refScale->SetMarkerStyle(24);
0195       refScale->SetMarkerColor(kRed);
0196       refScale->SetMarkerSize(1.0);
0197       refRes->SetMarkerStyle(24);
0198       refRes->SetMarkerColor(kRed);
0199       refRes->SetMarkerSize(1.0);
0200
0201       TLegend *ptresleg = new TLegend(0.4,0.77,0.5,0.9);
0202       ptresleg->AddEntry(refScale,"Reference","P");
0203       ptresleg->AddEntry(newScale,"New","P");
0204
0205       rawResCan->cd(3);
0206       refScale->Draw("psame");
0207       ptresleg->Draw("same");
0208
0209       rawResCan->cd(4);
0210       refRes->Draw("psame");
0211       ptresleg->Draw("same");
0212
0213       rawResCan->cd(1);
0214
0215       gPad->SetRightMargin(0.2);
0216       gPad->SetLogx();
0217       gPad->SetLogz();
0218       h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref->SetTitle("Reference pT Spectrum");
0219       h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref->Draw("colz");
0220       TLatex refl;
0221       refl.SetTextSize(0.05);
0222       refl.SetNDC();
0223       refl.SetTextColor(kBlack);
0224       refl.DrawLatex(0.45,0.96,"Reference");
0225
0226     }
0227
0228     rawResCan->Draw();
0229
0230     outfilef->cd();
0231     rawResCan->Write();
0232
0233   }
0234
0235
0236
0237   {
0238   //base histogram from the reco module name
0239   const char *hist_name_prefix = "QAG4SimulationTracking";
0240   TString prefix = TString("h_") + hist_name_prefix + TString("_");
0241
0242   // obtain normalization
0243   double Nevent_new = 1;
0244   double Nevent_ref = 1;
0245
0246   TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_TruthMatchingOverview") +
0247                                 TString("_") + hist_name_prefix,
0248                             TString("QA_Draw_Tracking_TruthMatchingOverview") +
0249                                 TString("_") + hist_name_prefix,
0250                             1800, 1000);
0251   c1->Divide(3, 1);
0252   int idx = 1;
0253   TPad *p;
0254
0255   {
0256     static const int nrebin = 5;
0257
0258     p = (TPad *)c1->cd(idx++);
0259     c1->Update();
0260     p->SetLogx();
0261     p->SetGridy();
0262
0263     TH1 *h_pass =
0264         (TH1 *)qa_file_new->GetObjectChecked(prefix + "nReco_pTGen", "TH1");
0265     TH1 *h_norm =
0266         (TH1 *)qa_file_new->GetObjectChecked(prefix + "nGen_pTGen", "TH1");
0267     assert(h_norm);
0268     assert(h_pass);
0269
0270       h_norm->SetDirectory(nullptr);
0271       h_pass->SetDirectory(nullptr);
0272
0273     h_norm->Rebin(nrebin);
0274     h_pass->Rebin(nrebin);
0275
0276     TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);
0277
0278     //    h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);
0279     h_ratio->GetYaxis()->SetTitle("Reco efficiency");
0280     h_ratio->GetYaxis()->SetRangeUser(-0, 1.);
0281
0282     TH1 *h_ratio_ref = NULL;
0283     if (qa_file_ref) {
0284       TH1 *h_pass =
0285           (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nReco_pTGen", "TH1");
0286       TH1 *h_norm =
0287           (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nGen_pTGen", "TH1");
0288       assert(h_norm);
0289       assert(h_pass);
0290       h_norm->SetDirectory(nullptr);
0291       h_pass->SetDirectory(nullptr);
0292       h_norm->Rebin(nrebin);
0293       h_pass->Rebin(nrebin);
0294       h_ratio_ref = GetBinominalRatio(h_pass, h_norm);
0295     }
0296
0297     h_ratio->SetTitle(TString(hist_name_prefix) + ": Tracking Efficiency");
0298
0299     DrawReference(h_ratio, h_ratio_ref, false);
0300
0301   }
0302
0303   {
0304     static const int nrebin = 4;
0305
0306     p = (TPad *)c1->cd(idx++);
0307     c1->Update();
0308     // p->SetLogx();
0309     p->SetGridy();
0310
0311     TH1 *h_pass =
0312         (TH1 *)qa_file_new->GetObjectChecked(prefix + "nReco_etaGen", "TH1");
0313     TH1 *h_norm =
0314         (TH1 *)qa_file_new->GetObjectChecked(prefix + "nGen_etaGen", "TH1");
0315     assert(h_norm);
0316     assert(h_pass);
0317
0318       h_norm->SetDirectory(nullptr);
0319       h_pass->SetDirectory(nullptr);
0320     h_norm->Rebin(nrebin);
0321     h_pass->Rebin(nrebin);
0322
0323     TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);
0324
0325     h_ratio->GetXaxis()->SetRangeUser(-1.1, 1.1);
0326     h_ratio->GetYaxis()->SetTitle("Reco efficiency");
0327     h_ratio->GetYaxis()->SetRangeUser(-0, 1.);
0328
0329     TH1 *h_ratio_ref = NULL;
0330     if (qa_file_ref) {
0331       TH1 *h_pass =
0332           (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nReco_etaGen", "TH1");
0333       TH1 *h_norm =
0334           (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nGen_etaGen", "TH1");
0335       assert(h_norm);
0336       assert(h_pass);
0337       h_norm->SetDirectory(nullptr);
0338       h_pass->SetDirectory(nullptr);
0339       h_norm->Rebin(nrebin);
0340       h_pass->Rebin(nrebin);
0341       h_ratio_ref = GetBinominalRatio(h_pass, h_norm);
0342     }
0343
0344     h_ratio->SetTitle(TString(hist_name_prefix) + ": Tracking Efficiency");
0345
0346     DrawReference(h_ratio, h_ratio_ref, false);
0347   }
0348
0349   {
0350     p = (TPad *)c1->cd(idx++);
0351     c1->Update();
0352     //    p->SetLogx();
0353     TH1 *frame = p->DrawFrame(0, .9, 50, 1.1,
0354                               "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})");
0355     //gPad->SetLeftMargin(.2);
0356     gPad->SetTopMargin(-1);
0357     frame->GetYaxis()->SetTitleOffset(1.7);
0358     //TLine *l = new TLine(0, 1, 50, 1);
0359     //l->SetLineColor(kGray);
0360     //l->Draw();
0361     HorizontalLine( gPad, 1 )->Draw();
0362
0363     TH2 *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen =
0364         (TH2 *)qa_file_new->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",
0365                                              "TH2");
0366     assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);
0367
0368       h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->SetDirectory(nullptr);
0369     h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Rebin2D(16, 1);
0370     TGraphErrors *ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen =
0371         FitProfile(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);
0372     ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Draw("pe");
0373     ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen->SetTitle(
0374         "Mean and sigma, p_{T,reco}/p_{T,truth}");
0375
0376     TGraphErrors *h_ratio_ref = NULL;
0377     if (qa_file_ref) {
0378       TH2 *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen =
0379           (TH2 *)qa_file_ref->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",
0380                                                "TH2");
0381       assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);
0382
0383       h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->SetDirectory(nullptr);
0384       h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Rebin2D(16, 1);
0385
0386       h_ratio_ref = FitProfile(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);
0387       ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Draw("pe");
0388     }
0389
0390     DrawReference(ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen, h_ratio_ref,
0391                   true);
0392
0393   }
0394
0395   //SaveCanvas(c1,
0396   //           TString(qa_file_name_new) + TString("_") + TString(c1->GetName()),
0397   //           true);
0398
0399   c1->Draw();
0400   outfilef->cd();
0401   c1->Write();
0402 }
0403
0404
0405 {
0406     const char *hist_name_prefix = "QAG4SimulationTracking";
0407     TString prefix = TString("h_") + hist_name_prefix + TString("_");
0408
0409   // obtain normalization
0410   double Nevent_new = 1;
0411   double Nevent_ref = 1;
0412
0413
0414   TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(
0415       prefix + TString("pTRecoGenRatio_pTGen"), "TH2");
0416   assert(h_new);
0417
0418   //  h_new->Rebin(1, 2);
0419   //h_new->Sumw2();
0420   //  h_new->Scale(1. / Nevent_new);
0421
0422   TH2 *h_ref = NULL;
0423   if (qa_file_ref)
0424   {
0425     h_ref = (TH2 *) qa_file_ref->GetObjectChecked(
0426         prefix + TString("pTRecoGenRatio_pTGen"), "TH2");
0427     assert(h_ref);
0428
0429     //    h_ref->Rebin(1, 2);
0430     //h_ref->Sumw2();
0431     h_ref->Scale(Nevent_new / Nevent_ref);
0432   }
0433
0434   TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_pTRatio") + TString("_") + hist_name_prefix,
0435                             TString("QA_Draw_Tracking_pTRatio") + TString("_") + hist_name_prefix,
0436                             1800, 1000);
0437   c1->Divide(4, 2);
0438   int idx = 1;
0439   TPad *p;
0440
0441   vector<pair<double, double>> gpt_ranges{
0442       {0, 1},
0443       {1, 5},
0444       {5, 10},
0445       {10, 20},
0446       {20, 30},
0447       {30, 40},
0448       {40, 45},
0449       {45, 50}};
0450   TF1 *f1 = nullptr;
0451   TF1 *fit = nullptr;
0452   Double_t sigma = 0;
0453   Double_t sigma_unc = 0;
0454   char resstr[500];
0455   TLatex *res = nullptr;
0456   for (auto pt_range : gpt_ranges)
0457   {
0458     //cout << __PRETTY_FUNCTION__ << " process " << pt_range.first << " - " << pt_range.second << " GeV/c";
0459
0460     p = (TPad *) c1->cd(idx++);
0461     c1->Update();
0462     p->SetLogy();
0463
0464     const double epsilon = 1e-6;
0465     const int bin_start = h_new->GetXaxis()->FindBin(pt_range.first + epsilon);
0466     const int bin_end = h_new->GetXaxis()->FindBin(pt_range.second - epsilon);
0467
0468     TH1 *h_proj_new = h_new->ProjectionY(
0469         TString::Format(
0470             "%s_New_ProjX_%d_%d",
0471             h_new->GetName(), bin_start, bin_end),
0472         bin_start, bin_end);
0473
0474     h_proj_new->GetXaxis()->SetRangeUser(.7, 1.3);
0475     h_proj_new->SetTitle(TString(hist_name_prefix) + TString::Format(
0476                                                          ": %.1f - %.1f GeV/c", pt_range.first, pt_range.second));
0477     h_proj_new->GetXaxis()->SetTitle(TString::Format(
0478         "Reco p_{T}/Truth p_{T}"));
0479
0480     f1 = new TF1("f1", "gaus", -.85, 1.15);
0481     h_proj_new->Fit(f1, "mq");
0482     sigma = f1->GetParameter(2);
0483     sigma_unc = f1->GetParError(2);
0484
0485     TH1 *h_proj_ref = nullptr;
0486     if (h_ref)
0487       h_proj_ref =
0488           h_ref->ProjectionY(
0489               TString::Format(
0490                   "%s_Ref_ProjX_%d_%d",
0491                   h_new->GetName(), bin_start, bin_end),
0492               bin_start, bin_end);
0493
0494     DrawReference(h_proj_new, h_proj_ref);
0495     sprintf(resstr, "#sigma = %.5f #pm %.5f", sigma, sigma_unc);
0496     res = new TLatex(0.325, 0.825, resstr);
0497     res->SetNDC();
0498     res->SetTextSize(0.05);
0499     res->SetTextAlign(13);
0500     res->Draw();
0501   }
0502
0503  // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);
0504     c1->Draw();
0505     outfilef->cd();
0506     c1->Write();
0507 }
0508
0509
0510  {
0511 const char *hist_name_prefix = "QAG4SimulationTracking";
0512     TString prefix = TString("h_") + hist_name_prefix + TString("_");
0513
0514
0515   // obtain normalization
0516   double Nevent_new = 1;
0517   double Nevent_ref = 1;
0518
0519   if (qa_file_new)
0520   {
0521     //cout << "Open new QA file " << qa_file_new->GetName() << endl;
0522
0523     TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(
0524         prefix + TString("Normalization"), "TH1");
0525     assert(h_norm);
0526
0527     Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
0528   }
0529   if (qa_file_ref)
0530   {
0531    // cout << "Open ref QA file " << qa_file_ref->GetName() << endl;
0532     TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(
0533         prefix + TString("Normalization"), "TH1");
0534     assert(h_norm);
0535
0536     Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
0537   }
0538
0539
0540   TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_RecoTruthMatching") +
0541                                 TString("_") + hist_name_prefix,
0542                             TString("QA_Draw_Tracking_RecoTruthMatching") +
0543                                 TString("_") + hist_name_prefix,
0544                             1800, 1000);
0545   c1->Divide(2, 1);
0546   int idx = 1;
0547   TPad *p;
0548
0549   {
0550     static const int nrebin = 5;
0551
0552     p = (TPad *) c1->cd(idx++);
0553     c1->Update();
0554     p->SetLogx();
0555     p->SetGridy();
0556
0557     TH1 *h_pass =
0558         (TH1 *) qa_file_new->GetObjectChecked(prefix + "nGen_pTReco", "TH1");
0559     TH1 *h_norm =
0560         (TH1 *) qa_file_new->GetObjectChecked(prefix + "nReco_pTReco", "TH1");
0561     assert(h_norm);
0562     assert(h_pass);
0563
0564       h_norm->SetDirectory(nullptr);
0565       h_pass->SetDirectory(nullptr);
0566
0567     h_norm->Rebin(nrebin);
0568     h_pass->Rebin(nrebin);
0569
0570     TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);
0571
0572     //    h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);
0573     h_ratio->GetYaxis()->SetTitle("Tracking Purity");
0574     h_ratio->GetYaxis()->SetRangeUser(-0, 1.1);
0575
0576     TH1 *h_ratio_ref = NULL;
0577     if (qa_file_ref)
0578     {
0579       TH1 *h_pass =
0580           (TH1 *) qa_file_ref->GetObjectChecked(prefix + "nGen_pTReco", "TH1");
0581       TH1 *h_norm =
0582           (TH1 *) qa_file_ref->GetObjectChecked(prefix + "nReco_pTReco", "TH1");
0583       assert(h_norm);
0584       assert(h_pass);
0585       h_norm->SetDirectory(nullptr);
0586       h_pass->SetDirectory(nullptr);
0587       h_norm->Rebin(nrebin);
0588       h_pass->Rebin(nrebin);
0589       h_ratio_ref = GetBinominalRatio(h_pass, h_norm);
0590     }
0591
0592     h_ratio->SetTitle("Tracking Purity (matched truth-reco pairs)");
0593
0594     DrawReference(h_ratio, h_ratio_ref, false);
0595   }
0596
0597   {
0598     p = (TPad *) c1->cd(idx++);
0599     c1->Update();
0600     //    p->SetLogx();
0601     TH1 *frame = p->DrawFrame(0, .9, 50, 1.1,
0602                               "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})");
0603     // gPad->SetLeftMargin(.2);
0604     gPad->SetTopMargin(-1);
0605     frame->GetYaxis()->SetTitleOffset(1.7);
0606     // TLine *l = new TLine(0, 1, 50, 1);
0607     // l->SetLineColor(kGray);
0608     // l->Draw();
0609     HorizontalLine(gPad, 1)->Draw();
0610
0611     TH2 *h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =
0612         (TH2 *) qa_file_new->GetObjectChecked(
0613             prefix + "pTRecoTruthMatchedRatio_pTReco", "TH2");
0614     assert(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);
0615
0616       h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetDirectory(nullptr);
0617     h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Rebin2D(16, 1);
0618
0619     TGraphErrors *ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =
0620         FitProfile(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);
0621     ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Draw("pe");
0622     ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetTitle(
0623         "Mean and sigma p_{Tmatched}/p_{Treco}");
0624
0625     TGraphErrors *h_ratio_ref = NULL;
0626     if (qa_file_ref)
0627     {
0628       TH2 *h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =
0629           (TH2 *) qa_file_ref->GetObjectChecked(
0630               prefix + "pTRecoTruthMatchedRatio_pTReco", "TH2");
0631       assert(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);
0632
0633       h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetDirectory(nullptr);
0634       h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Rebin2D(16, 1);
0635
0636       h_ratio_ref =
0637           FitProfile(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);
0638       ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Draw("pe");
0639     }
0640
0641     DrawReference(ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco,
0642                   h_ratio_ref, true);
0643   }
0644
0645     c1->Draw();
0646     outfilef->cd();
0647     c1->Write();
0648 }
0649
0650
0651 {
0652     const char *hist_name_prefix = "QAG4SimulationTracking";
0653       TString prefix = TString("h_") + hist_name_prefix + TString("_");
0654
0655   // obtain normalization
0656   double Nevent_new = 1;
0657   double Nevent_ref = 1;
0658
0659
0660
0661   if (qa_file_new)
0662   {
0663     TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(
0664         prefix + TString("Normalization"), "TH1");
0665     assert(h_norm);
0666
0667     Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
0668   }
0669   if (qa_file_ref)
0670   {
0671     TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(
0672         prefix + TString("Normalization"), "TH1");
0673     assert(h_norm);
0674
0675     Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
0676   }
0677
0678   //MVTX, INTT, TPC
0679   vector<TString> detectors{"MVTX", "INTT", "TPC"};
0680   vector<int> eff_ncluster_cuts{2, 2, 40};
0681   vector<double> ncluster_spectrum_pt_cuts{2, 2, 2};
0682   vector<TH2 *> h_pass_detectors(3, nullptr);
0683   static const int nrebin = 5;
0684
0685   h_pass_detectors[0] = (TH2 *) qa_file_new->GetObjectChecked(
0686       prefix + "nMVTX_nReco_pTGen", "TH1") ;
0687   h_pass_detectors[1] = (TH2 *) qa_file_new->GetObjectChecked(
0688       prefix + "nINTT_nReco_pTGen", "TH1") ;
0689   h_pass_detectors[2] = (TH2 *) qa_file_new->GetObjectChecked(
0690       prefix + "nTPC_nReco_pTGen", "TH1") ;
0691
0692   TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(
0693       prefix + "nGen_pTGen", "TH1") ;
0694   assert(h_norm);
0695     h_norm->SetDirectory(nullptr);
0696   h_norm->Rebin(nrebin);
0697
0698   vector<TH2 *> h_pass_detectors_ref(3, nullptr);
0699   TH1 *h_norm_ref = nullptr;
0700   if (qa_file_ref)
0701   {
0702     h_pass_detectors_ref[0] = (TH2 *) qa_file_ref->GetObjectChecked(
0703         prefix + "nMVTX_nReco_pTGen", "TH1") ;
0704     h_pass_detectors_ref[1] = (TH2 *) qa_file_ref->GetObjectChecked(
0705         prefix + "nINTT_nReco_pTGen", "TH1") ;
0706     h_pass_detectors_ref[2] = (TH2 *) qa_file_ref->GetObjectChecked(
0707         prefix + "nTPC_nReco_pTGen", "TH1") ;
0708
0709     h_norm_ref = (TH1 *) qa_file_ref->GetObjectChecked(
0710         prefix + "nGen_pTGen", "TH1") ;
0711     h_norm_ref->SetDirectory(nullptr);
0712     h_norm_ref->Rebin(nrebin);
0713
0714   }
0715
0716   TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_TruthMatching_NumOfClusters") + TString("_") + hist_name_prefix,
0717                             TString("QA_Draw_Tracking_TruthMatching_NumOfClusters") + TString("_") + hist_name_prefix,
0718                             1800, 1000);
0719   c1->Divide(3, 2);
0720   TPad *p;
0721
0722   for (int i = 0; i < 3; ++i)
0723   {
0724     TString detector = detectors[i];
0725     TH2 *h_pass_detector = h_pass_detectors[i];
0726     TH2 *h_pass_detector_ref = h_pass_detectors_ref[i];
0727     assert(h_pass_detector);
0728
0729     {
0730       p = (TPad *) c1->cd(i + 1);
0731       c1->Update();
0732       p->SetLogy();
0733
0734       const int bin_start = h_pass_detector->GetXaxis()->FindBin(ncluster_spectrum_pt_cuts[i]);
0735
0736       TH1 *h_pass_detector_ncluster = h_pass_detector->ProjectionY(
0737           TString(h_pass_detector->GetName()) + "_nCluster_new",
0738           bin_start);
0739       TH1 *h_pass_detector_ncluster_ref = nullptr;
0740       if (h_pass_detector_ref)
0741       {
0742         h_pass_detector_ncluster_ref = h_pass_detector_ref->ProjectionY(
0743             TString(h_pass_detector_ref->GetName()) + "_nCluster_ref",
0744             bin_start);
0745       }
0746
0747       h_pass_detector_ncluster->SetTitle(TString(hist_name_prefix) + ": " + detector + Form(" n_{Cluster} | p_{T} #geq %.1fGeV/c", ncluster_spectrum_pt_cuts[i]));
0748       h_pass_detector_ncluster->SetYTitle("# of reconstructed track");
0749       DrawReference(h_pass_detector_ncluster, h_pass_detector_ncluster_ref, false);
0750     }
0751
0752     {
0753       p = (TPad *) c1->cd(i + 3 + 1);
0754       c1->Update();
0755       p->SetLogx();
0756       p->SetGridy();
0757
0758       const int bin_start = h_pass_detector->GetYaxis()->FindBin(eff_ncluster_cuts[i]);
0759       TH1 *h_pass = h_pass_detector->ProjectionX(
0760           TString(h_pass_detector->GetName()) + "_nReco_new",
0761           bin_start);
0762
0763       assert(h_pass);
0764         h_pass->SetDirectory(nullptr);
0765       h_pass->Rebin(nrebin);
0766
0767       TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);
0768       h_ratio->GetYaxis()->SetTitle("Reco efficiency | " + detector + Form(" n_{Cluster} #geq %d", eff_ncluster_cuts[i]));
0769       h_ratio->GetYaxis()->SetRangeUser(-0, 1.);
0770       //
0771       TH1 *h_ratio_ref = NULL;
0772       if (h_pass_detector_ref)
0773       {
0774         TH1 *h_pass = h_pass_detector_ref->ProjectionX(
0775             TString(h_pass_detector->GetName()) + "_nReco_ref",
0776             bin_start);
0777
0778         assert(h_pass);
0779         h_pass->SetDirectory(nullptr);
0780       h_pass->Rebin(nrebin);
0781
0782         h_ratio_ref = GetBinominalRatio(h_pass, h_norm_ref);
0783       }
0784       //
0785       h_ratio->SetTitle("Tracking efficiency | " + detector + Form(" n_{Cluster} #geq %d", eff_ncluster_cuts[i]));
0786       DrawReference(h_ratio, h_ratio_ref, false);
0787     }
0788   }
0789
0790   // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);
0791     c1->Draw();
0792     outfilef->cd();
0793     c1->Write();
0794 }
0795
0796
0797
0798  {
0799      const char *hist_name_prefix = "QAG4SimulationTracking";
0800       TString prefix = TString("h_") + hist_name_prefix + TString("_");
0801
0802   auto c1 = new TCanvas(TString("QA_Draw_Tracking_nClus_Layer") + TString("_") + hist_name_prefix,
0803                         TString("QA_Draw_Tracking_nClus_Layer") + TString("_") + hist_name_prefix,
0804                         1800, 1000);
0805
0806   c1->Divide(2, 1);
0807   c1->cd(1);
0808   GetNormalization(qa_file_new, qa_file_ref, prefix, "Truth Track");
0809   Draw(qa_file_new, qa_file_ref, prefix, "nClus_layerGen");
0810
0811   c1->cd(2);
0812   GetNormalization(qa_file_new, qa_file_ref, prefix, "Reco Track");
0813   Draw(qa_file_new, qa_file_ref, prefix, "nClus_layer");
0814
0815   // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);
0816     c1->Draw();
0817     outfilef->cd();
0818     c1->Write();
0819 }
0820
0821
0822 {
0823   const char *hist_name_prefix = "QAG4SimulationUpsilon";
0824   TString prefix = TString("h_") + hist_name_prefix + TString("_");
0825
0826   // obtain normalization
0827   double Nevent_new = 1;
0828   double Nevent_ref = 1;
0829
0830   if ( qa_file_new->GetObjectChecked(
0831         prefix + TString("pTRecoGenRatio_pTGen"), "TH2")
0832      == nullptr )
0833     {
0834         cout <<"QAG4SimulationUpsilon is not enabled. Skip...."<<endl;
0835     }
0836     else
0837     {
0838
0839       TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_UpsilonOverview") + TString("_") + hist_name_prefix,
0840                                 TString("QA_Draw_Tracking_UpsilonOverview") + TString("_") + hist_name_prefix,
0841                                 1800, 1000);
0842       c1->Divide(2, 1);
0843       int idx = 1;
0844       TPad *p;
0845
0846       {
0847         p = (TPad *) c1->cd(idx++);
0848         c1->Update();
0849         p->SetLogy();
0850
0851         TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(
0852             prefix + TString("pTRecoGenRatio_pTGen"), "TH2");
0853         assert(h_new);
0854
0855         //  h_new->Rebin(1, 2);
0856         //h_new->Sumw2();
0857         //  h_new->Scale(1. / Nevent_new);
0858
0859         TH2 *h_ref = NULL;
0860         if (qa_file_ref)
0861         {
0862           h_ref = (TH2 *) qa_file_ref->GetObjectChecked(
0863               prefix + TString("pTRecoGenRatio_pTGen"), "TH2");
0864           assert(h_ref);
0865
0866           //    h_ref->Rebin(1, 2);
0867           //h_ref->Sumw2();
0868           h_ref->Scale(Nevent_new / Nevent_ref);
0869         }
0870
0871         TH1 *h_proj_new = h_new->ProjectionY(
0872             TString::Format(
0873                 "%s_New_ProjX",
0874                 h_new->GetName()));
0875
0876         h_proj_new->GetXaxis()->SetRangeUser(0, 1.3);
0877         h_proj_new->SetTitle(TString(hist_name_prefix) + TString::Format(
0878                                                              ": Electron lineshape"));
0879         h_proj_new->GetXaxis()->SetTitle(TString::Format(
0880             "Reco p_{T}/Truth p_{T}"));
0881
0882         TF1 *f_eLineshape = new TF1("f_eLineshape", CBcalc, 7, 11, 5);
0883         f_eLineshape->SetParameter(0, 1.0);
0884         f_eLineshape->SetParameter(1, 1.0);
0885         f_eLineshape->SetParameter(2, 0.95);
0886         f_eLineshape->SetParameter(3, 0.08);
0887         f_eLineshape->SetParameter(4, 20.0);
0888
0889         f_eLineshape->SetParNames("alpha","n","m","sigma","N");
0890         f_eLineshape->SetLineColor(kRed);
0891         f_eLineshape->SetLineWidth(3);
0892         f_eLineshape->SetLineStyle(kSolid);
0893         f_eLineshape->SetNpx(1000);
0894
0895         h_proj_new->Fit(f_eLineshape);
0896
0897         TH1 *h_proj_ref = nullptr;
0898         if (h_ref)
0899         {
0900           h_proj_ref =
0901               h_ref->ProjectionY(
0902                   TString::Format(
0903                       "%s_Ref_ProjX",
0904                       h_new->GetName()));
0905         }
0906         TF1 *f_eLineshape_ref = new TF1("f_eLineshape_ref", CBcalc, 7, 11, 5);
0907         f_eLineshape_ref->SetParameter(0, 1.0);
0908         f_eLineshape_ref->SetParameter(1, 1.0);
0909         f_eLineshape_ref->SetParameter(2, 0.95);
0910         f_eLineshape_ref->SetParameter(3, 0.08);
0911         f_eLineshape_ref->SetParameter(4, 20.0);
0912
0913         f_eLineshape_ref->SetParNames("alpha","n","m","sigma","N");
0914         f_eLineshape_ref->SetLineColor(kRed);
0915         f_eLineshape_ref->SetLineWidth(3);
0916         f_eLineshape_ref->SetLineStyle(kSolid);
0917
0918         h_proj_ref->Fit(f_eLineshape_ref);
0919
0920
0921         DrawReference(h_proj_new, h_proj_ref);
0922         f_eLineshape->Draw("same");
0923
0924         char resstr_1[500];
0925         sprintf(resstr_1,"#sigma_{dp/p} = %.2f #pm %.2f %%", f_eLineshape->GetParameter(3)*100, f_eLineshape->GetParError(3)*100);
0926         TLatex *res_1 = new TLatex(0.2,0.75,resstr_1);
0927         res_1->SetNDC();
0928         res_1->SetTextSize(0.05);
0929         res_1->SetTextAlign(13);
0930         res_1->Draw();
0931
0932         char resstr_2[500];
0933         sprintf(resstr_2,"#sigma_{dp/p,ref} = %.2f #pm %.2f %%", f_eLineshape_ref->GetParameter(3)*100, f_eLineshape_ref->GetParError(3)*100);
0934         TLatex *res_2 = new TLatex(0.2,0.7,resstr_2);
0935         res_2->SetNDC();
0936         res_2->SetTextSize(0.05);
0937         res_2->SetTextAlign(13);
0938         res_2->Draw();
0939       }
0940
0941       {
0942         p = (TPad *) c1->cd(idx++);
0943         c1->Update();
0944     //    p->SetLogy();
0945
0946         TH1 *h_new = (TH1 *) qa_file_new->GetObjectChecked(
0947             prefix + TString("nReco_Pair_InvMassReco"), "TH1");
0948         assert(h_new);
0949
0950         // h_new->Rebin(2);
0951         // h_new->Sumw2();
0952         // h_new->Scale(1. / Nevent_new);
0953
0954         TF1 *f1S = new TF1("f1S", CBcalc2, 7, 11, 7);
0955         f1S->SetParameter(0, 50.0);
0956         f1S->SetParameter(1, 9.46);
0957         f1S->SetParameter(2, 0.08);
0958         f1S->SetParameter(3, 1.0);
0959         f1S->SetParameter(4, 3.0);
0960         f1S->SetParameter(5, 1.0);
0961         f1S->SetParLimits(3, 0.120, 10);
0962         f1S->SetParLimits(4, 1.05, 10);
0963         f1S->SetParameter(6, 3.0);
0964         f1S->SetParLimits(5, 0.120, 10);
0965         f1S->SetParLimits(6, 1.05, 10);
0966
0967         f1S->SetParNames("N", "m", "#sigma", "#alpha_{left}","n_{left}","#alpha_{right}","#sigma_{right}");
0968         f1S->SetLineColor(kRed);
0969         f1S->SetLineWidth(3);
0970         f1S->SetLineStyle(kSolid);
0971         f1S->SetNpx(1000);
0972
0973         h_new->Fit(f1S);
0974
0975         TH1 *h_ref = NULL;
0976         if (qa_file_ref)
0977         {
0978           h_ref = (TH1 *) qa_file_ref->GetObjectChecked(
0979               prefix + TString("nReco_Pair_InvMassReco"), "TH1");
0980           assert(h_ref);
0981
0982           // h_ref->Rebin(2);
0983           // h_ref->Sumw2();
0984           // h_ref->Scale(Nevent_new / Nevent_ref);
0985         }
0986
0987         h_new->SetTitle(TString(hist_name_prefix) + TString::Format(
0988                                                         ": #Upsilon #rightarrow e^{+}e^{-} lineshape"));
0989         h_new->GetXaxis()->SetRangeUser(7, 10);
0990
0991         TF1 *f1S_ref = new TF1("f1S_ref", CBcalc2, 7, 11, 7);
0992         f1S_ref->SetParameter(0, 50.0);
0993         f1S_ref->SetParameter(1, 9.46);
0994         f1S_ref->SetParameter(2, 0.08);
0995         f1S_ref->SetParameter(3, 1.0);
0996         f1S_ref->SetParameter(4, 3.0);
0997         f1S_ref->SetParameter(5, 1.0);
0998         f1S_ref->SetParLimits(3, 0.120, 10);
0999         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 }