EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KaonIDStudy.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KaonIDStudy.C
1 #include "TROOT.h"
2 #include "TChain.h"
3 #include "TEfficiency.h"
4 #include "TH1.h"
5 #include "TGraphErrors.h"
6 #include "TCut.h"
7 #include "TCanvas.h"
8 #include "TStyle.h"
9 #include "TLegend.h"
10 #include "TMath.h"
11 #include "TLine.h"
12 #include "TLatex.h"
13 #include "TRatioPlot.h"
14 #include "TString.h"
15 
16 #include <glob.h>
17 #include <iostream>
18 #include <iomanip>
19 #include <vector>
20 
21 #include "PlotFunctions.h"
22 
23 void ConfigHistogram(TH1 *h, std::string which = "charm")
24 {
25  if (which == "charm") {
26  h->SetLineColor(kBlue);
27  h->SetMarkerColor(kBlue);
28  h->SetMarkerStyle(kOpenSquare);
29  h->SetMarkerColor(kBlue);
30  } else if (which == "light") {
31  h->SetLineColor(kRed);
32  h->SetMarkerColor(kRed);
33  h->SetMarkerStyle(kFullCrossX);
34  h->SetMarkerColor(kRed);
35  }
36 }
37 
38 void KaonIDStudy(TString dir,
39  TString input,
40  TString trackname = "barrelDircTrack",
41  TString filePattern = "*/out.root")
42 {
43  auto data = new TChain("tree");
44 
45  data->SetTitle(input.Data());
46  auto files = fileVector(Form("%s/%s/%s", dir.Data(), input.Data(), filePattern.Data()));
47 
48  for (auto file : files)
49  {
50  data->Add(file.c_str());
51  }
52 
53 
54  // TString trackname = "CF4RICHTrack";
55  // TString trackname = "mRICHTrack";
56  // TString trackname = "barrelDircTrack";
57  // TString trackname = "tofBarrelTrack";
58  // TString trackname = "dRICHagTrack";
59 
60  Float_t etaMin = -8;
61  Float_t etaMax = 8;
62 
63  TCut *cut_fiducial = nullptr;
64  Float_t xmin = 0.0;
65  Float_t xmax = 55.0;
66  Bool_t logplot = kFALSE;
67 
68  if (trackname == "mRICHTrack") {
69  // mRICH
70  etaMin = -3.5;
71  etaMax = -1.0;
72  xmax = 12.0;
73  } else if (trackname == "barrelDIRCTrack") {
74  // barrel DIRC
75  etaMin = -1.0;
76  etaMax = 1.0;
77  xmax = 55.0;
78  } else if (trackname == "CF4RICHTrack") {
79  // CF4RICH
80  etaMin = 1.0;
81  etaMax = 4.0;
82  } else if (trackname == "tofBarrelTrack") {
83  // TOF Barrel Detector
84  etaMin = -2.0;
85  etaMax = 2.0;
86  } else if (trackname == "dualRICHTrack") {
87  // dualRICH, aerogel-based
88  etaMin = 1.00;
89  etaMax = 3.50;
90 
91  // etaMin = 1.48;
92  // etaMax = 3.91;
93  xmax = 50.0;
94  logplot = kTRUE;
95 
96  }
97  cut_fiducial = new TCut(Form("%0.1f < pid_track_eta && pid_track_eta < %0.1f && pid_track_pt>0.1", etaMin, etaMax));
98 
99  Float_t xbinsize = 0.5;
100  Int_t nbinsx = TMath::Ceil((xmax - xmin) / xbinsize);
101 
102  TCut cut_truekaon("TMath::Abs(pid_track_true_pid) == 321");
103  TCut cut_recokaon("TMath::Abs(pid_track_reco_pid) == 321");
104  TCut cut_truepion("TMath::Abs(pid_track_true_pid) == 211");
105  TCut cut_recopion("TMath::Abs(pid_track_reco_pid) == 211");
106 
107  plot_config draw_config;
108  draw_config.htemplate = new TH1F("TrackPT",
109  "",
110  nbinsx,
111  xmin,
112  xmax);
113 
114  // draw_config.xlimits = std::vector<float>();
115  // draw_config.xlimits.push_back(0);
116  // draw_config.xlimits.push_back(20);
117  draw_config.ylimits = std::vector < float > ();
118  draw_config.ylimits.push_back(0.0);
119  draw_config.ylimits.push_back(1e4);
120  draw_config.xtitle = "Track Momentum [GeV]";
121  draw_config.ytitle = "Frequency";
122  draw_config.logy = kFALSE;
123  draw_config.logx = kFALSE;
124  draw_config.cuts = "";
125 
126  // Kaon -> Kaon ID
127  auto true_kaon_pt = GeneratePlot(draw_config, data, "True Kaon PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && cut_truekaon);
128  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);
129 
130 
131  TH1F *h_template = static_cast < TH1F * > (reco_kaon_pt->Clone("h_template"));
132  h_template->SetAxisRange(0.0, 1.12, "Y");
133  if (logplot)
134  h_template->SetAxisRange(1.0e-5, 10.0, "Y");
135 
136  h_template->GetYaxis()->SetTitle("Efficiency");
137 
138  auto pid_efficiency = new TEfficiency(*reco_kaon_pt,
139  *true_kaon_pt);
140  pid_efficiency->SetMarkerColor(kBlue);
141  pid_efficiency->SetFillColor(kBlue);
142  pid_efficiency->SetLineColor(kBlue);
143  pid_efficiency->SetMarkerStyle(kOpenSquare);
144 
145  // Pion -> Kaon ID
146  auto true_pion_pt = GeneratePlot(draw_config, data, "True Pion PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial && cut_truepion);
147  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);
148 
149 
150  auto misid_efficiency = new TEfficiency(*reco_pion_pt,
151  *true_pion_pt);
152  misid_efficiency->SetMarkerColor(kRed);
153  misid_efficiency->SetFillColor(kRed);
154  misid_efficiency->SetLineColor(kRed);
155  misid_efficiency->SetMarkerStyle(kFullCrossX);
156 
157  // Draw!
158  TCanvas c1("c1",
159  "",
160  800,
161  600);
162 
163  c1.cd();
164  c1.SetLogy(logplot);
165 
166  // true_kaon_pt->Draw("HIST");
167  // reco_kaon_pt->Draw("E1 SAME");
168  h_template->Draw("AXIS");
169  pid_efficiency->Draw("E1 SAME");
170  misid_efficiency->Draw("E1 SAME");
171 
172  TLatex etaRange((xmin + (xmax - xmin) * 0.67),
173  1.05,
174  Form("%0.1f < #eta < %0.1f", etaMin, etaMax));
175  etaRange.Draw();
176 
177  // legend
178  auto *legend = smart_legend("lower left");
179  legend->AddEntry(pid_efficiency, "K^{#pm} #rightarrow K^{#pm}", "lp");
180  legend->AddEntry(misid_efficiency, "#pi^{#pm} #rightarrow K^{#pm}", "lp");
181  legend->Draw();
182 
183  c1.SaveAs(Form("%s_KaonIDStudy.pdf", trackname.Data()));
184 
185 
186  return;
187 
188  // Plot the pT spectrum for true and reconstructed candidates from charm and
189  // light jets
190 
191  // Kaon -> Kaon ID
192  TCut true_charm_mother("pid_track_jetmother == 4");
193  TCut true_light_mother("pid_track_jetmother > 0 && pid_track_jetmother < 4");
194 
195  if (cut_fiducial != nullptr) delete cut_fiducial;
196  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"));
197 
198 
199  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);
200  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);
201 
202  c1.Clear();
203 
204  charm_true_kaon_pt->SetYTitle("Probability");
205 
206  ConfigHistogram(charm_true_kaon_pt, "charm");
207  ConfigHistogram(light_true_kaon_pt, "light");
208 
209  charm_true_kaon_pt->Rebin(2);
210  light_true_kaon_pt->Rebin(2);
211 
212  TH1 *charm_true_kaon = charm_true_kaon_pt->DrawNormalized("E1");
213  TH1 *light_true_kaon = light_true_kaon_pt->DrawNormalized("E1 SAME");
214 
215  charm_true_kaon->SetAxisRange(1e-5, 1.0, "Y");
216 
217  c1.SetLogy();
218  c1.SaveAs(Form("True_Kaon_Momentum.pdf"));
219 
220 
221  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);
222  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);
223 
224  c1.Clear();
225 
226  charm_reco_kaon_pt->SetYTitle("Probability");
227 
228  ConfigHistogram(charm_reco_kaon_pt, "charm");
229  ConfigHistogram(light_reco_kaon_pt, "light");
230 
231  charm_reco_kaon_pt->Rebin(2);
232  light_reco_kaon_pt->Rebin(2);
233 
234  TH1 *charm_reco_kaon = charm_reco_kaon_pt->DrawNormalized("E1");
235  TH1 *light_reco_kaon = light_reco_kaon_pt->DrawNormalized("E1 SAME");
236 
237  charm_reco_kaon->SetAxisRange(1e-5, 1.0, "Y");
238 
239  c1.SetLogy();
240  c1.SaveAs(Form("Reco_Kaon_Momentum.pdf"));
241 
242  //
243  // Zoom in on the low-momentum region, where most of the true kaons in CC DIS
244  // live.
245  //
246 
247  plot_config kaon_low_p_config;
248  kaon_low_p_config.htemplate = new TH1F("TrackMomentum_Low",
249  "",
250  100,
251  0.0,
252  10.0);
253 
254  // kaon_low_p_config.xlimits = std::vector < float > ();
255  // kaon_low_p_config.xlimits.push_back(0);
256  // kaon_low_p_config.xlimits.push_back(10);
257  kaon_low_p_config.ylimits = std::vector < float > ();
258  kaon_low_p_config.ylimits.push_back(0.0);
259  kaon_low_p_config.ylimits.push_back(1.0e6);
260  kaon_low_p_config.xtitle = "Track Momentum [GeV]";
261  kaon_low_p_config.ytitle = "Frequency";
262  kaon_low_p_config.logy = kFALSE;
263  kaon_low_p_config.logx = kFALSE;
264  kaon_low_p_config.cuts = "";
265 
266  // auto reco_kaon_pt = GeneratePlot(draw_config, data, "Reconstructed Kaon
267  // PT", "pid_track_pt*TMath::CosH(pid_track_eta)", *cut_fiducial &&
268  // cut_truekaon && cut_recokaon);
269 
270  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);
271  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);
272 
273  std::cout << charm_true_kaon_p_low->GetNbinsX() << std::endl;
274  std::cout << light_true_kaon_p_low->GetNbinsX() << std::endl;
275 
276  c1.Clear();
277  c1.SetLogy(kaon_low_p_config.logy);
278  c1.SetLogx(kaon_low_p_config.logx);
279 
280  charm_true_kaon_p_low->SetYTitle("Probability");
281 
282  ConfigHistogram(charm_true_kaon_p_low, "charm");
283  ConfigHistogram(light_true_kaon_p_low, "light");
284 
285  charm_true_kaon = charm_true_kaon_p_low->DrawNormalized("E1");
286  light_true_kaon = light_true_kaon_p_low->DrawNormalized("E1 SAME");
287 
288  charm_true_kaon->SetAxisRange(0.0, 0.04, "Y");
289 
290  c1.SaveAs(Form("True_Kaon_Momentum_Low.pdf"));
291 
292  // Zoomed, reconstructed momentum
293  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);
294  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);
295 
296  std::cout << charm_reco_kaon_p_low->GetNbinsX() << std::endl;
297  std::cout << light_reco_kaon_p_low->GetNbinsX() << std::endl;
298 
299  c1.Clear();
300  c1.SetLogy(kaon_low_p_config.logy);
301  c1.SetLogx(kaon_low_p_config.logx);
302 
303  charm_reco_kaon_p_low->SetYTitle("Probability");
304 
305  ConfigHistogram(charm_reco_kaon_p_low, "charm");
306  ConfigHistogram(light_reco_kaon_p_low, "light");
307 
308  charm_reco_kaon = charm_reco_kaon_p_low->DrawNormalized("E1");
309  light_reco_kaon = light_reco_kaon_p_low->DrawNormalized("E1 SAME");
310 
311  charm_reco_kaon->SetAxisRange(0.0, 0.04, "Y");
312 
313  c1.SaveAs(Form("Reco_Kaon_Momentum_Low.pdf"));
314 }