SimpleAnalysis/scripts/KaonIDStudy.C

0001 #include "TROOT.h"
0002 #include "TChain.h"
0003 #include "TEfficiency.h"
0004 #include "TH1.h"
0005 #include "TGraphErrors.h"
0006 #include "TCut.h"
0007 #include "TCanvas.h"
0008 #include "TStyle.h"
0009 #include "TLegend.h"
0010 #include "TMath.h"
0011 #include "TLine.h"
0012 #include "TLatex.h"
0013 #include "TRatioPlot.h"
0014 #include "TString.h"
0015
0016 #include <glob.h>
0017 #include <iostream>
0018 #include <iomanip>
0019 #include <vector>
0020
0021 #include "PlotFunctions.h"
0022
0023 void ConfigHistogram(TH1 *h, std::string which = "charm")
0024 {
0025   if (which == "charm") {
0026     h->SetLineColor(kBlue);
0027     h->SetMarkerColor(kBlue);
0028     h->SetMarkerStyle(kOpenSquare);
0029     h->SetMarkerColor(kBlue);
0030   } else if (which == "light") {
0031     h->SetLineColor(kRed);
0032     h->SetMarkerColor(kRed);
0033     h->SetMarkerStyle(kFullCrossX);
0034     h->SetMarkerColor(kRed);
0035   }
0036 }
0037
0038 void KaonIDStudy(TString dir,
0039                  TString input,
0040                  TString trackname   = "barrelDircTrack",
0041                  TString filePattern = "*/out.root")
0042 {
0043   auto data = new TChain("tree");
0044
0045   data->SetTitle(input.Data());
0046   auto files = fileVector(Form("%s/%s/%s", dir.Data(), input.Data(), filePattern.Data()));
0047
0048   for (auto file : files)
0049   {
0050     data->Add(file.c_str());
0051   }
0052
0053
0054   // TString trackname = "CF4RICHTrack";
0055   // TString trackname = "mRICHTrack";
0056   // TString trackname = "barrelDircTrack";
0057   // TString trackname = "tofBarrelTrack";
0058   // TString trackname = "dRICHagTrack";
0059
0060   Float_t etaMin = -8;
0061   Float_t etaMax = 8;
0062
0063   TCut   *cut_fiducial = nullptr;
0064   Float_t xmin         = 0.0;
0065   Float_t xmax         = 55.0;
0066   Bool_t logplot = kFALSE;
0067
0068   if (trackname == "mRICHTrack") {
0069     // mRICH
0070     etaMin = -3.5;
0071     etaMax = -1.0;
0072     xmax   = 12.0;
0073   } else if (trackname == "barrelDIRCTrack") {
0074     // barrel DIRC
0075     etaMin = -1.0;
0076     etaMax = 1.0;
0077     xmax   = 55.0;
0078   } else if (trackname == "CF4RICHTrack") {
0079     // CF4RICH
0080     etaMin = 1.0;
0081     etaMax = 4.0;
0082   } else if (trackname == "tofBarrelTrack") {
0083     // TOF Barrel Detector
0084     etaMin = -2.0;
0085     etaMax = 2.0;
0086   } else if (trackname == "dualRICHTrack") {
0087     // dualRICH, aerogel-based
0088     etaMin = 1.00;
0089     etaMax = 3.50;
0090
0091     // etaMin       = 1.48;
0092     // etaMax       = 3.91;
0093     xmax = 50.0;
0094     logplot = kTRUE;
0095
0096   }
0097   cut_fiducial = new TCut(Form("%0.1f < pid_track_eta && pid_track_eta < %0.1f && pid_track_pt>0.1", etaMin, etaMax));
0098
0099   Float_t xbinsize = 0.5;
0100   Int_t   nbinsx   = TMath::Ceil((xmax - xmin) / xbinsize);
0101
0102   TCut                        cut_truekaon("TMath::Abs(pid_track_true_pid) == 321");
0103   TCut                        cut_recokaon("TMath::Abs(pid_track_reco_pid) == 321");
0104   TCut                        cut_truepion("TMath::Abs(pid_track_true_pid) == 211");
0105   TCut                        cut_recopion("TMath::Abs(pid_track_reco_pid) == 211");
0106
0107   plot_config draw_config;
0108   draw_config.htemplate = new TH1F("TrackPT",
0109                                    "",
0110                                    nbinsx,
0111                                    xmin,
0112                                    xmax);
0113
0114   // draw_config.xlimits   = std::vector<float>();
0115   // draw_config.xlimits.push_back(0);
0116   // draw_config.xlimits.push_back(20);
0117   draw_config.ylimits = std::vector < float > ();
0118   draw_config.ylimits.push_back(0.0);
0119   draw_config.ylimits.push_back(1e4);
0120   draw_config.xtitle = "Track Momentum [GeV]";
0121   draw_config.ytitle = "Frequency";
0122   draw_config.logy   = kFALSE;
0123   draw_config.logx   = kFALSE;
0124   draw_config.cuts   = "";
0125
0126   // Kaon -> Kaon ID
0127   auto true_kaon_pt = GeneratePlot(draw_config, data, "True Kaon PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && cut_truekaon);
0128   auto reco_kaon_pt = GeneratePlot(draw_config, data, "Reconstructed Kaon PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && cut_truekaon && cut_recokaon);
0129
0130
0131   TH1F *h_template = static_cast < TH1F * > (reco_kaon_pt->Clone("h_template"));
0132   h_template->SetAxisRange(0.0, 1.12, "Y");
0133   if (logplot)
0134     h_template->SetAxisRange(1.0e-5, 10.0, "Y");
0135
0136   h_template->GetYaxis()->SetTitle("Efficiency");
0137
0138   auto pid_efficiency = new TEfficiency(*reco_kaon_pt,
0139                                         *true_kaon_pt);
0140   pid_efficiency->SetMarkerColor(kBlue);
0141   pid_efficiency->SetFillColor(kBlue);
0142   pid_efficiency->SetLineColor(kBlue);
0143   pid_efficiency->SetMarkerStyle(kOpenSquare);
0144
0145   // Pion -> Kaon ID
0146   auto true_pion_pt = GeneratePlot(draw_config, data, "True Pion PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && cut_truepion);
0147   auto reco_pion_pt = GeneratePlot(draw_config, data, "Reconstructed Pion PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && cut_truepion && cut_recokaon);
0148
0149
0150   auto misid_efficiency = new TEfficiency(*reco_pion_pt,
0151                                           *true_pion_pt);
0152   misid_efficiency->SetMarkerColor(kRed);
0153   misid_efficiency->SetFillColor(kRed);
0154   misid_efficiency->SetLineColor(kRed);
0155   misid_efficiency->SetMarkerStyle(kFullCrossX);
0156
0157   // Draw!
0158   TCanvas c1("c1",
0159              "",
0160              800,
0161              600);
0162
0163   c1.cd();
0164   c1.SetLogy(logplot);
0165
0166   // true_kaon_pt->Draw("HIST");
0167   // reco_kaon_pt->Draw("E1 SAME");
0168   h_template->Draw("AXIS");
0169   pid_efficiency->Draw("E1 SAME");
0170   misid_efficiency->Draw("E1 SAME");
0171
0172   TLatex etaRange((xmin + (xmax - xmin) * 0.67),
0173                   1.05,
0174                   Form("%0.1f < #eta < %0.1f", etaMin, etaMax));
0175   etaRange.Draw();
0176
0177   // legend
0178   auto *legend = smart_legend("lower left");
0179   legend->AddEntry(pid_efficiency,   "K^{#pm} #rightarrow K^{#pm}",   "lp");
0180   legend->AddEntry(misid_efficiency, "#pi^{#pm} #rightarrow K^{#pm}", "lp");
0181   legend->Draw();
0182
0183   c1.SaveAs(Form("%s_KaonIDStudy.pdf", trackname.Data()));
0184
0185
0186   return;
0187
0188   // Plot the pT spectrum for true and reconstructed candidates from charm and
0189   // light jets
0190
0191   // Kaon -> Kaon ID
0192   TCut               true_charm_mother("pid_track_jetmother == 4");
0193   TCut               true_light_mother("pid_track_jetmother > 0 && pid_track_jetmother < 4");
0194
0195   if (cut_fiducial != nullptr) delete cut_fiducial;
0196   cut_fiducial = new TCut(Form("-4.0 < pid_track_eta && pid_track_eta < 4.0 && pid_track_pt>0.1 && pid_track_jet_pt > 10.0 && TMath::Abs(pid_track_jet_eta)<3.0"));
0197
0198
0199   auto charm_true_kaon_pt = GeneratePlot(draw_config, data, "True Kaon PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && true_charm_mother && cut_truekaon);
0200   auto light_true_kaon_pt = GeneratePlot(draw_config, data, "True Kaon PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && true_light_mother && cut_truekaon);
0201
0202   c1.Clear();
0203
0204   charm_true_kaon_pt->SetYTitle("Probability");
0205
0206   ConfigHistogram(charm_true_kaon_pt, "charm");
0207   ConfigHistogram(light_true_kaon_pt, "light");
0208
0209   charm_true_kaon_pt->Rebin(2);
0210   light_true_kaon_pt->Rebin(2);
0211
0212   TH1 *charm_true_kaon = charm_true_kaon_pt->DrawNormalized("E1");
0213   TH1 *light_true_kaon = light_true_kaon_pt->DrawNormalized("E1 SAME");
0214
0215   charm_true_kaon->SetAxisRange(1e-5, 1.0, "Y");
0216
0217   c1.SetLogy();
0218   c1.SaveAs(Form("True_Kaon_Momentum.pdf"));
0219
0220
0221   auto charm_reco_kaon_pt = GeneratePlot(draw_config, data, "reco Kaon PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && true_charm_mother && cut_truekaon && cut_recokaon);
0222   auto light_reco_kaon_pt = GeneratePlot(draw_config, data, "True Kaon PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && true_light_mother && cut_truekaon && cut_recokaon);
0223
0224   c1.Clear();
0225
0226   charm_reco_kaon_pt->SetYTitle("Probability");
0227
0228   ConfigHistogram(charm_reco_kaon_pt, "charm");
0229   ConfigHistogram(light_reco_kaon_pt, "light");
0230
0231   charm_reco_kaon_pt->Rebin(2);
0232   light_reco_kaon_pt->Rebin(2);
0233
0234   TH1 *charm_reco_kaon = charm_reco_kaon_pt->DrawNormalized("E1");
0235   TH1 *light_reco_kaon = light_reco_kaon_pt->DrawNormalized("E1 SAME");
0236
0237   charm_reco_kaon->SetAxisRange(1e-5, 1.0, "Y");
0238
0239   c1.SetLogy();
0240   c1.SaveAs(Form("Reco_Kaon_Momentum.pdf"));
0241
0242   //
0243   // Zoom in on the low-momentum region, where most of the true kaons in CC DIS
0244   // live.
0245   //
0246
0247   plot_config kaon_low_p_config;
0248   kaon_low_p_config.htemplate = new TH1F("TrackMomentum_Low",
0249                                          "",
0250                                          100,
0251                                          0.0,
0252                                          10.0);
0253
0254   // kaon_low_p_config.xlimits = std::vector < float > ();
0255   // kaon_low_p_config.xlimits.push_back(0);
0256   // kaon_low_p_config.xlimits.push_back(10);
0257   kaon_low_p_config.ylimits = std::vector < float > ();
0258   kaon_low_p_config.ylimits.push_back(0.0);
0259   kaon_low_p_config.ylimits.push_back(1.0e6);
0260   kaon_low_p_config.xtitle = "Track Momentum [GeV]";
0261   kaon_low_p_config.ytitle = "Frequency";
0262   kaon_low_p_config.logy   = kFALSE;
0263   kaon_low_p_config.logx   = kFALSE;
0264   kaon_low_p_config.cuts   = "";
0265
0266   // auto reco_kaon_pt = GeneratePlot(draw_config, data, "Reconstructed Kaon
0267   // PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial &&
0268   // cut_truekaon && cut_recokaon);
0269
0270   auto charm_true_kaon_p_low = GeneratePlot(kaon_low_p_config, data, "Charm True Kaon P", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && true_charm_mother && cut_truekaon);
0271   auto light_true_kaon_p_low = GeneratePlot(kaon_low_p_config, data, "Light True Kaon P", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && true_light_mother && cut_truekaon);
0272
0273   std::cout << charm_true_kaon_p_low->GetNbinsX() << std::endl;
0274   std::cout << light_true_kaon_p_low->GetNbinsX() << std::endl;
0275
0276   c1.Clear();
0277   c1.SetLogy(kaon_low_p_config.logy);
0278   c1.SetLogx(kaon_low_p_config.logx);
0279
0280   charm_true_kaon_p_low->SetYTitle("Probability");
0281
0282   ConfigHistogram(charm_true_kaon_p_low, "charm");
0283   ConfigHistogram(light_true_kaon_p_low, "light");
0284
0285   charm_true_kaon = charm_true_kaon_p_low->DrawNormalized("E1");
0286   light_true_kaon = light_true_kaon_p_low->DrawNormalized("E1 SAME");
0287
0288   charm_true_kaon->SetAxisRange(0.0, 0.04, "Y");
0289
0290   c1.SaveAs(Form("True_Kaon_Momentum_Low.pdf"));
0291
0292   // Zoomed, reconstructed momentum
0293   auto charm_reco_kaon_p_low = GeneratePlot(kaon_low_p_config, data, "Charm True Kaon P", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && true_charm_mother && cut_truekaon && cut_recokaon);
0294   auto light_reco_kaon_p_low = GeneratePlot(kaon_low_p_config, data, "Light True Kaon P", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && true_light_mother && cut_truekaon && cut_recokaon);
0295
0296   std::cout << charm_reco_kaon_p_low->GetNbinsX() << std::endl;
0297   std::cout << light_reco_kaon_p_low->GetNbinsX() << std::endl;
0298
0299   c1.Clear();
0300   c1.SetLogy(kaon_low_p_config.logy);
0301   c1.SetLogx(kaon_low_p_config.logx);
0302
0303   charm_reco_kaon_p_low->SetYTitle("Probability");
0304
0305   ConfigHistogram(charm_reco_kaon_p_low, "charm");
0306   ConfigHistogram(light_reco_kaon_p_low, "light");
0307
0308   charm_reco_kaon = charm_reco_kaon_p_low->DrawNormalized("E1");
0309   light_reco_kaon = light_reco_kaon_p_low->DrawNormalized("E1 SAME");
0310
0311   charm_reco_kaon->SetAxisRange(0.0, 0.04, "Y");
0312
0313   c1.SaveAs(Form("Reco_Kaon_Momentum_Low.pdf"));
0314 }