EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CharmTagPlots.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CharmTagPlots.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 "TFile.h"
15 
16 #include <glob.h>
17 #include <iostream>
18 #include <iomanip>
19 #include <vector>
20 
21 #include "PlotFunctions.h"
22 
23 
24 float target_xsec = LookupCrossSection("CC_DIS");
25 
26 float target_lumi = 100 * u_fb;
27 
28 float n_gen = -1;
29 
30 
31 TEfficiency DifferentialTaggingEfficiency(TTree *data, plot_config draw_config, TString xvar, TString which = "all")
32 {
33  std::cout << "DifferentialTaggingEfficiency: processing " << which.Data() << " ... " << std::endl;
34 
35  TH1F *tru_yield = static_cast<TH1F *>(draw_config.htemplate->Clone(Form("tru_yield_%d", histUID)));
36  TH1F *tag_yield = static_cast<TH1F *>(draw_config.htemplate->Clone(Form("tag_yield_%d", histUID)));
37 
38  TCut *tru_selection = nullptr;
39  TCut *tag_selection = nullptr;
40 
41  if (which == "light") {
42  tru_selection = new TCut("jet_flavor < 4 || jet_flavor == 21");
43  tag_selection = new TCut(*tru_selection && TCut(draw_config.cuts));
44  } else if (which == "charm") {
45  tru_selection = new TCut("jet_flavor == 4");
46  tag_selection = new TCut(*tru_selection && TCut(draw_config.cuts));
47  }
48 
49 
50  data->Project(tru_yield->GetName(), xvar.Data(), *tru_selection);
51  data->Project(tag_yield->GetName(), xvar.Data(), *tag_selection);
52 
53  std::cout << " True Yield: " << tru_yield->Integral() << std::endl;
54  std::cout << " Tag Yield: " << tag_yield->Integral() << std::endl;
55 
56  TEfficiency eff(*tag_yield, *tru_yield);
57 
58  if (tru_yield)
59  delete tru_yield;
60 
61  if (tag_yield)
62  delete tag_yield;
63 
64  histUID++;
65 
66  return eff;
67 }
68 
69 TH1F* DifferentialTaggingYield(TTree *data, plot_config draw_config, TString xvar, TString which = "all")
70 {
71  std::cout << "DifferentialTaggingYield: processing " << which.Data() << " ... " << std::endl;
72 
73  TH1F *tag_yield = static_cast<TH1F *>(draw_config.htemplate->Clone(Form("%s_%s_tag_yield_%d", which.Data(), xvar.Data(), histUID)));
74  tag_yield->Sumw2();
75 
76  TCut *tag_selection = nullptr;
77 
78  if (which == "light") {
79  tag_selection = new TCut(TCut("(jet_flavor < 4 || jet_flavor == 21)") && TCut(draw_config.cuts));
80  } else if (which == "charm") {
81  tag_selection = new TCut(TCut("jet_flavor == 4") && TCut(draw_config.cuts));
82  }
83 
84  // Add other analysis-level cuts to the tag selection
85  *tag_selection = (*tag_selection) && TCut("met_et > 10.0");
86 
87 
88  data->Project(tag_yield->GetName(), xvar.Data(), *tag_selection);
89 
90  std::cout << " Unnormalized Tag Yield: " << tag_yield->GetEntries() << std::endl;
91  tag_yield->Scale(target_xsec * target_lumi / data->GetEntries());
92  std::cout << " Normalized Tag Yield: " << tag_yield->Integral() << std::endl;
93 
94 
95  // auto tag_graph = new TGraphErrors(tag_yield);
96 
97  histUID++;
98  return tag_yield;
99 }
100 
101 void DrawTagEfficiencyPlot(TCanvas *pad, TTree *data, plot_config draw_config, TString xvar)
102 {
103  auto light_eff = DifferentialTaggingEfficiency(data, draw_config, xvar, "light");
104  auto charm_eff = DifferentialTaggingEfficiency(data, draw_config, xvar, "charm");
105 
106 
107  // make a template histogram to fine-tune layout of the final plot
108  TH1F *htemplate = new TH1F("htemplate", "", 1, draw_config.xlimits[0], draw_config.xlimits[1]);
109 
110  set_axis_range(htemplate, draw_config.ylimits[0], draw_config.ylimits[1], "Y");
111  set_axis_range(htemplate, draw_config.xlimits[0], draw_config.xlimits[1], "X");
112 
113  // Configure plots
114  configure_plot<TEfficiency>(&charm_eff, draw_config, "charm");
115  configure_plot<TEfficiency>(&light_eff, draw_config, "light");
116 
117  // Draw Layout
118  htemplate->Draw("HIST");
119  charm_eff.Draw("P E1 SAME");
120  light_eff.Draw("E1 P SAME");
121 
122  // Title
123  TLatex plot_title = make_title();
124  plot_title.Draw("SAME");
125 
126  // Configure the Pad
127  pad->SetLogy(draw_config.logy);
128  pad->SetLogx(draw_config.logx);
129  pad->SetGrid(1, 1);
130 
131 
132  // Configure the Legend
133  TLegend *legend = smart_legend("lower right");
134  legend->SetFillStyle(0);
135  legend->SetBorderSize(0);
136  legend->AddEntry(&charm_eff, "Charm Jets", "lp");
137  legend->AddEntry(&light_eff, "Light Jets", "lp");
138  legend->Draw();
139 
140  pad->Update();
141  pad->SaveAs(Form("CharmTagPlot_tagging_efficiency_%s_%s.pdf", data->GetTitle(), xvar.Data()));
142  pad->SaveAs(Form("CharmTagPlot_tagging_efficiency_%s_%s.C", data->GetTitle(), xvar.Data()));
143 
144  // Cleanup
145  if (htemplate != nullptr)
146  delete htemplate;
147 
148  if (legend != nullptr)
149  delete legend;
150 }
151 
152 void DrawTagYieldPlot(TCanvas *pad, TTree *data, plot_config draw_config, TString xvar)
153 {
154  std::cout << "Drawing tag yield plot for variable " << xvar << std::endl;
155  auto charm_yield = DifferentialTaggingYield(data, draw_config, xvar, "charm");
156  auto light_yield = DifferentialTaggingYield(data, draw_config, xvar, "light");
157 
158  // Configure plots
159 
160  configure_plot<TH1F>(charm_yield, draw_config, "charm");
161  configure_plot<TH1F>(light_yield, draw_config, "light");
162 
163  // generate template histogram
164  TH1F *htemplate = new TH1F("htemplate", "", 1, draw_config.xlimits[0], draw_config.xlimits[1]);
165 
166  set_axis_range(htemplate, draw_config.ylimits[0], draw_config.ylimits[1], "Y");
167  set_axis_range(htemplate, draw_config.xlimits[0], draw_config.xlimits[1], "X");
168  htemplate->SetXTitle(draw_config.xtitle);
169  htemplate->SetYTitle(draw_config.ytitle);
170 
171  pad->SetGrid(1, 1);
172  pad->SetLogy(draw_config.logy);
173  pad->SetLogx(draw_config.logx);
174 
175  htemplate->Draw("HIST");
176  charm_yield->Draw("E1 P SAME");
177  light_yield->Draw("E1 P SAME");
178 
179  TLatex plot_title = make_title();
180  plot_title.Draw("SAME");
181 
182 
183  // Legend
184  TLegend *legend = smart_legend("upper right");
185  legend->SetFillStyle(0);
186  legend->SetBorderSize(0);
187  legend->AddEntry(charm_yield, "Charm Jets", "lp");
188  legend->AddEntry(light_yield, "Light Jets", "lp");
189  legend->Draw();
190  pad->Update();
191 
192 
193  pad->SaveAs(Form("CharmTagPlot_tagging_yield_%s_%s.pdf", data->GetTitle(), xvar.Data()));
194  pad->SaveAs(Form("CharmTagPlot_tagging_yield_%s_%s.C", data->GetTitle(), xvar.Data()));
195 
196  // Cleanup
197  if (htemplate != nullptr)
198  delete htemplate;
199 
200  if (legend != nullptr)
201  delete legend;
202 }
203 
204 float CutEfficiency(TTree *data, TCut cut, TString cutdesc = "")
205 {
206  if (n_gen < 0) {
207  n_gen = static_cast<float>(data->GetEntries());
208  }
209 
210  float n_yield = data->GetEntries(cut.GetTitle());
211 
212  TString desc = cutdesc;
213 
214  if (desc == "") {
215  desc = cut.GetTitle();
216  }
217 
218  std::cout << setw(30) << desc << " efficiency: " << setprecision(4) << n_yield / n_gen * 100.0 << "%" << std::endl;
219 
220  return n_yield / n_gen;
221 }
222 
223 void CharmTagPlots(TString dir, TString input, TString xvar, TString filePattern = "*/out.root")
224 {
225  // Global options
226  gStyle->SetOptStat(0);
227 
228  // Create the TCanvas
229  TCanvas *pad = new TCanvas("pad", "", 800, 600);
230  TLegend *legend = nullptr;
231  TH1F *htemplate = nullptr;
232 
233  auto default_data = new TChain("tree");
234  default_data->SetTitle(input.Data());
235  auto files = fileVector(Form("%s/%s/%s", dir.Data(), input.Data(), filePattern.Data()));
236 
237  for (auto file : files)
238  {
239  default_data->Add(file.c_str());
240  }
241 
242 
243  // Some basic cut efficiency information
244  std::cout << "Inclusive Jet Efficiency Information: " << std::endl;
245  TCut selection("jet_n>0");
246  CutEfficiency(default_data, selection, "Jet Reconstructed");
247 
248  selection = selection && TCut("TMath::Abs(jet_eta) < 3.0");
249  CutEfficiency(default_data, selection, "Jet Fiducial");
250 
251  selection = selection && TCut("met_et > 10.0");
252  CutEfficiency(default_data, selection, "MET > 10 GeV");
253 
254  TCut main_preselection = selection;
255 
256  // Print charm jet MET cut efficiency
257  std::cout << "Charm Jet Efficiency Information: " << std::endl;
258  selection = TCut("jet_n > 0 && TMath::Abs(jet_eta) < 3.0 && jet_flavor == 4");
259  CutEfficiency(default_data, selection, "Charm Jet Reconstructed");
260 
261  selection = selection && TCut("met_et > 10.0");
262  CutEfficiency(default_data, selection, "MET > 10 GeV");
263 
264 
265  TTree *default_data_selected = default_data->CopyTree(main_preselection.GetTitle());
266  default_data_selected->SetTitle(input.Data());
267 
268  bool do_TagEffPlot = kFALSE;
269  bool do_TagYieldPlot = kFALSE; // kTRUE;
270  bool do_100fbProjPlot = kTRUE; // kTRUE;
271  bool do_Helicity = kFALSE; // kTRUE;
272 
273 
274  // plot configurations
275  plot_config draw_config;
276 
277  if (xvar == "jet_pt") {
278  float xbins[] = { 10, 12.5, 15, 17.5, 20, 25, 35, 60 };
279  int nbins = sizeof(xbins) / sizeof(xbins[0]) - 1;
280  draw_config.htemplate = new TH1F(xvar, "", nbins, xbins);
281  draw_config.xlimits = std::vector<float>();
282  draw_config.xlimits.push_back(0);
283  draw_config.xlimits.push_back(60);
284  draw_config.ylimits = std::vector<float>();
285  draw_config.ylimits.push_back(1e-3);
286  draw_config.ylimits.push_back(1.0);
287  draw_config.xtitle = "Reconstructed Jet p_{T} [GeV]";
288  draw_config.ytitle = "#varepsilon_{tag}";
289  draw_config.logy = kTRUE;
290  draw_config.logx = kFALSE;
291  draw_config.cuts = "jet_sip3dtag==1";
292  } else if (xvar == "bjorken_x") {
293  float xbins[] = { 0.01, 0.043333, 0.076666, 0.1, 0.25, 0.5, 1.0 };
294  int nbins = sizeof(xbins) / sizeof(xbins[0]) - 1;
295  draw_config.htemplate = new TH1F(xvar, "", nbins, xbins);
296  draw_config.xlimits = std::vector<float>();
297  draw_config.xlimits.push_back(0.01);
298  draw_config.xlimits.push_back(1.0);
299  draw_config.ylimits = std::vector<float>();
300  draw_config.ylimits.push_back(1e-3);
301  draw_config.ylimits.push_back(1.0);
302  draw_config.xtitle = "Bjorken x";
303  draw_config.ytitle = "#varepsilon_{tag}";
304  draw_config.logy = kTRUE;
305  draw_config.logx = kTRUE;
306  draw_config.cuts = "jet_sip3dtag==1";
307  } else if (xvar == "jb_x") {
308  float xbins[] = { 0.01, 0.043333, 0.076666, 0.1, 0.25, 0.5, 1.0 };
309  int nbins = sizeof(xbins) / sizeof(xbins[0]) - 1;
310  draw_config.htemplate = new TH1F(xvar, "", nbins, xbins);
311  draw_config.xlimits = std::vector<float>();
312  draw_config.xlimits.push_back(0.01);
313  draw_config.xlimits.push_back(1.0);
314  draw_config.ylimits = std::vector<float>();
315  draw_config.ylimits.push_back(1e-3);
316  draw_config.ylimits.push_back(1.0);
317  draw_config.xtitle = "Reconstructed x_{JB}";
318  draw_config.ytitle = "#varepsilon_{tag}";
319  draw_config.logy = kTRUE;
320  draw_config.logx = kTRUE;
321  draw_config.cuts = "jet_sip3dtag==1";
322  } else {
323  do_TagEffPlot = kFALSE;
324  }
325 
326 
327  //
328  // Efficiency Plot
329  //
330 
331  if (do_TagEffPlot) {
332  DrawTagEfficiencyPlot(pad, default_data_selected, draw_config, xvar);
333  }
334 
335  //
336  // Cleanup
337  //
338 
339  pad->Clear();
340 
341  //
342  // Yield plot
343  //
344  gStyle->SetErrorX(0.5);
345 
346  if (xvar == "jet_pt") {
347  draw_config.ylimits = std::vector<float>();
348 
349  // draw_config.ylimits.push_back(0.1);
350  // draw_config.ylimits.push_back(5000);
351  draw_config.ylimits.push_back(0.0);
352  draw_config.ylimits.push_back(2500);
353  draw_config.logy = kFALSE;
354 
355  draw_config.xtitle = "Reconstructed Jet p_{T} [GeV]";
356  draw_config.ytitle = "#varepsilon_{tag} #times #sigma_{CC-DIS} #times 100fb^{-1}";
357  draw_config.cuts = "jet_sip3dtag==1";
358  } else if (xvar == "bjorken_x") {
359  draw_config.ylimits = std::vector<float>();
360  draw_config.ylimits.push_back(0.1);
361  draw_config.ylimits.push_back(5000);
362  draw_config.xtitle = "Bjorken x";
363  draw_config.ytitle = "#varepsilon_{tag} #times #sigma_{CC-DIS} #times 100fb^{-1}";
364  draw_config.cuts = "jet_sip3dtag==1";
365  } else if (xvar == "jb_x") {
366  draw_config.ylimits = std::vector<float>();
367  draw_config.ylimits.push_back(0.1);
368  draw_config.ylimits.push_back(5000);
369  draw_config.xtitle = "Reconstructed x_{JB}";
370  draw_config.ytitle = "#varepsilon_{tag} #times #sigma_{CC-DIS} #times 100fb^{-1}";
371  draw_config.cuts = "jet_sip3dtag==1";
372  } else if (xvar == "jet_mlp_ktagger") {
373  draw_config.htemplate = new TH1F(xvar, "", 50, -0.00001, 1.00001);
374  draw_config.xlimits = std::vector<float>();
375  draw_config.xlimits.push_back(0);
376  draw_config.xlimits.push_back(1);
377  draw_config.ylimits = std::vector<float>();
378  draw_config.ylimits.push_back(0.1);
379  draw_config.ylimits.push_back(1e7);
380  draw_config.xtitle = "MLP Kaon-based Charm Jet Tagger Output";
381  draw_config.ytitle = "#varepsilon_{tag} #times #sigma_{CC-DIS} #times 100fb^{-1}";
382  draw_config.logy = kTRUE;
383  draw_config.logx = kFALSE;
384  draw_config.cuts = "1==1";
385  } else if (xvar == "jet_mlp_eltagger") {
386  draw_config.htemplate = new TH1F(xvar, "", 50, -0.00001, 1.00001);
387  draw_config.xlimits = std::vector<float>();
388  draw_config.xlimits.push_back(0);
389  draw_config.xlimits.push_back(1);
390  draw_config.ylimits = std::vector<float>();
391  draw_config.ylimits.push_back(0.1);
392  draw_config.ylimits.push_back(1e7);
393  draw_config.xtitle = "MLP Electron-based Charm Jet Tagger Output";
394  draw_config.ytitle = "#varepsilon_{tag} #times #sigma_{CC-DIS} #times 100fb^{-1}";
395  draw_config.logy = kTRUE;
396  draw_config.logx = kFALSE;
397  draw_config.cuts = "1==1";
398  } else if (xvar == "jet_mlp_mutagger") {
399  draw_config.htemplate = new TH1F(xvar, "", 50, -0.00001, 1.00001);
400  draw_config.xlimits = std::vector<float>();
401  draw_config.xlimits.push_back(0);
402  draw_config.xlimits.push_back(1);
403  draw_config.ylimits = std::vector<float>();
404  draw_config.ylimits.push_back(0.1);
405  draw_config.ylimits.push_back(1e7);
406  draw_config.xtitle = "MLP Muon-based Charm Jet Tagger Output";
407  draw_config.ytitle = "#varepsilon_{tag} #times #sigma_{CC-DIS} #times 100fb^{-1}";
408  draw_config.logy = kTRUE;
409  draw_config.logx = kFALSE;
410  draw_config.cuts = "1==1";
411  } else if (xvar == "jet_mlp_ip3dtagger") {
412  draw_config.htemplate = new TH1F(xvar, "", 50, -0.00001, 1.00001);
413  draw_config.xlimits = std::vector<float>();
414  draw_config.xlimits.push_back(0);
415  draw_config.xlimits.push_back(1);
416  draw_config.ylimits = std::vector<float>();
417  draw_config.ylimits.push_back(0.1);
418  draw_config.ylimits.push_back(1e7);
419  draw_config.xtitle = "MLP sIP_{3D}-based Charm Jet Tagger Output";
420  draw_config.ytitle = "#varepsilon_{tag} #times #sigma_{CC-DIS} #times 100fb^{-1}";
421  draw_config.logy = kTRUE;
422  draw_config.logx = kFALSE;
423  draw_config.cuts = "1==1";
424  } else if (xvar == "jet_mlp_globaltagger") {
425  draw_config.htemplate = new TH1F(xvar, "", 50, -0.00001, 1.00001);
426  draw_config.xlimits = std::vector<float>();
427  draw_config.xlimits.push_back(0);
428  draw_config.xlimits.push_back(1);
429  draw_config.ylimits = std::vector<float>();
430  draw_config.ylimits.push_back(0.1);
431  draw_config.ylimits.push_back(1e7);
432  draw_config.xtitle = "MLP Global Charm Jet Tagger Output";
433  draw_config.ytitle = "#varepsilon_{tag} #times #sigma_{CC-DIS} #times 100fb^{-1}";
434  draw_config.logy = kTRUE;
435  draw_config.logx = kFALSE;
436  draw_config.cuts = "1==1";
437  } else {
438  do_TagYieldPlot = kFALSE;
439  }
440 
441 
442  if (do_TagYieldPlot) {
443  DrawTagYieldPlot(pad, default_data_selected, draw_config, xvar);
444  }
445 
446  //
447  // Error Bands at 100/fb
448  //
449 
450  if (xvar == "jet_pt") {
451  draw_config.ylimits = std::vector<float>();
452  draw_config.ylimits.push_back(0.1);
453  draw_config.ylimits.push_back(5000);
454  draw_config.xtitle = "Reconstructed Jet p_{T} [GeV]";
455  draw_config.ytitle = "Relative Variation to Suppressed Strangeness";
456  } else if (xvar == "bjorken_x") {
457  draw_config.ytitle = "Relative Variation to Suppressed Strangeness";
458  } else if (xvar == "jb_x") {
459  draw_config.ytitle = "Relative Variation to Suppressed Strangeness";
460  } else {
461  do_100fbProjPlot = kFALSE;
462  }
463 
464  bool do_BkgSubtraction = kTRUE;
465 
466  if (do_100fbProjPlot == kTRUE) {
467  pad->Clear();
468  gStyle->SetErrorX(0.5);
469 
470  auto charm_yield = DifferentialTaggingYield(default_data_selected, draw_config, xvar, "charm");
471  auto charm_yield_100fb = static_cast<TH1F *>(charm_yield->Clone("charm_yield_100fb"));
472 
473  for (Int_t i = 1; i <= charm_yield_100fb->GetNbinsX(); i++) {
474  charm_yield_100fb->SetBinError(i, TMath::Sqrt(charm_yield_100fb->GetBinContent(i)));
475  }
476 
477  TH1F *light_yield = nullptr;
478  TH1F *light_yield_100fb = nullptr;
479 
480  if (do_BkgSubtraction == kTRUE) {
481  std::cout << " For Yield Estimate, performing a naive background subtraction error estimation..." << std::endl;
482  light_yield = DifferentialTaggingYield(default_data_selected, draw_config, xvar, "light");
483  light_yield_100fb = static_cast<TH1F *>(light_yield->Clone("light_yield_100fb"));
484 
485  for (Int_t i = 1; i <= light_yield_100fb->GetNbinsX(); i++) {
486  light_yield_100fb->SetBinError(i, TMath::Sqrt(light_yield_100fb->GetBinContent(i)));
487  }
488 
489  // inflate the charm yield errors to account for "background subtraction"
490  for (Int_t i = 1; i <= light_yield_100fb->GetNbinsX(); i++) {
491  float err_charm = charm_yield_100fb->GetBinError(i);
492  float err_light = light_yield_100fb->GetBinError(i);
493  float err_total = TMath::Sqrt(charm_yield_100fb->GetBinContent(i) + light_yield_100fb->GetBinContent(i));
494 
495  float new_charm_err = TMath::Sqrt(TMath::Power(err_total, 2.0) + TMath::Power(err_light, 2.0));
496 
497 
498  charm_yield_100fb->SetBinError(i, new_charm_err);
499  }
500  }
501 
502  std::map<TString, TString> alt_samples;
503  //alt_samples["CC_DIS_e10_p275_lha_21Rs2"] = "CT18ZNNLO with enhanced strangeness [R_{s}=0.863]";
504  alt_samples["CC_DIS_e10_p275_lha_21Rs2"] = "Rs-High NNLO (enhanced strange)";
505  alt_samples["CC_DIS_e10_p275_CT18ANNLO"] = "CT18A NNLO (intermediate strange)";
506  //alt_samples["CC_DIS_e10_p275_lha_22Rs2"] = "CT18ZNNLO with intermediate strangeness (22Rs2)";
507  //alt_samples["CC_DIS_e10_p275_lha_22Rs2ver3"] = "R_{s}^{INT}=0.594";
508  //alt_samples["CC_DIS_e10_p275_MMHT2014nnlo68cl"] = "MMHT2014nnlo68cl";
509  //alt_samples["CC_DIS_e10_p275_NNPDF31_nnlo_as_0118"] = "NNPDF31_nnlo_as_0118";
510 
511  std::map<TString, Int_t> alt_marker;
512  alt_marker["CC_DIS_e10_p275_lha_21Rs2"] = kOpenDiamond;
513  alt_marker["CC_DIS_e10_p275_CT18ANNLO"] = kOpenCrossX;
514  alt_marker["CC_DIS_e10_p275_lha_22Rs2ver3" ] = kOpenCross;
515 
516  std::map<TString, Int_t> alt_color;
517  alt_color["CC_DIS_e10_p275_lha_21Rs2"] = kBlue - 5;
518  alt_color["CC_DIS_e10_p275_CT18ANNLO"] = kViolet - 5;
519  alt_color["CC_DIS_e10_p275_lha_22Rs2ver3" ] = kSpring - 6;
520 
521 
522 
523 
524 
525  // Load alternative samples
526  auto ccdis_20Rs2_data = new TChain("tree");
527  files = fileVector(Form("%s/CC_DIS_e10_p275_lha_20Rs2/%s", dir.Data(), filePattern.Data()));
528  for (auto file : files)
529  {
530  ccdis_20Rs2_data->Add(file.c_str());
531  }
532 
533  std::cout << "Running pre-selection on alternative samples." << std::endl;
534  std::cout << "... suppressed strangeness ..." << std::endl;
535  TTree *ccdis_20Rs2_data_selected = ccdis_20Rs2_data->CopyTree(main_preselection.GetTitle());
536  auto ccdis_20Rs2_yield = DifferentialTaggingYield(ccdis_20Rs2_data_selected, draw_config, xvar, "charm");
537  auto ccdis_20Rs2_yield_100fb = static_cast<TH1F *>(ccdis_20Rs2_yield->Clone("20Rs2_yield_100fb"));
538 
539  for (Int_t i = 1; i <= charm_yield_100fb->GetNbinsX(); i++) {
540  // charm_yield_100fb->SetBinError(i,
541  // TMath::Sqrt(charm_yield_100fb->GetBinContent(i)) );
542  ccdis_20Rs2_yield_100fb->SetBinError(i, TMath::Sqrt(ccdis_20Rs2_yield_100fb->GetBinContent(i)));
543 
544  if (do_BkgSubtraction == kTRUE) {
545  float err_charm = ccdis_20Rs2_yield_100fb->GetBinError(i);
546  float err_light = light_yield_100fb->GetBinError(i);
547  float err_total = TMath::Sqrt(ccdis_20Rs2_yield_100fb->GetBinContent(i) + light_yield_100fb->GetBinContent(i));
548 
549  float new_charm_err = TMath::Sqrt(TMath::Power(err_total, 2.0) + TMath::Power(err_light, 2.0));
550  ccdis_20Rs2_yield_100fb->SetBinError(i, new_charm_err);
551  }
552 
553  }
554 
555  // Generate the error band histogram for the suppressed case
556  TH1F *error_band_100fb = static_cast<TH1F *>(ccdis_20Rs2_yield_100fb->Clone("error_band_100fb"));
557 
558  for (Int_t i = 1; i <= error_band_100fb->GetNbinsX(); i++) {
559  error_band_100fb->SetBinError(i, error_band_100fb->GetBinError(i) / error_band_100fb->GetBinContent(i));
560  error_band_100fb->SetBinContent(i, 1.0);
561  }
562 
563 
564  // Now handle the suppressed-only error band
565  error_band_100fb->SetLineColor(kGray);
566  error_band_100fb->SetFillColor(kGray);
567  error_band_100fb->SetMarkerSize(0.0001);
568  error_band_100fb->SetMarkerColor(kGray);
569  error_band_100fb->SetTitle("#splitline{Stat. Uncertainty}{[CT18NNLO, R_{s}=2s/(#bar{u} + #bar{d})= 0.325 (suppressed)]}");
570 
571  htemplate = new TH1F("htemplate", "", 1, draw_config.xlimits[0], draw_config.xlimits[1]);
572  htemplate->SetXTitle(draw_config.xtitle);
573  htemplate->SetYTitle(draw_config.ytitle);
574  htemplate->GetYaxis()->SetTitleSize(0.04);
575  set_axis_range(htemplate, 0.0, 2.0, "Y");
576  set_axis_range(htemplate, draw_config.xlimits[0], draw_config.xlimits[1], "X");
577 
578  htemplate->Draw("HIST");
579  error_band_100fb->Draw("E2 SAME");
580 
581 
582  // Legend
583  legend = smart_legend("lower left", 0.75, 0.28);
584  legend->SetTextSize(0.045);
585  legend->SetFillStyle(0);
586  legend->SetBorderSize(0);
587  //legend->AddEntry(error_band_100fb, "Stat. Uncertainty [CT18NNLO, R_{s}=2s/(#bar{u} + #bar{d})= 0.325 (suppressed)]", "lf");
588  legend->AddEntry(error_band_100fb, "#splitline{Stat. Uncertainty}{[Rs-Low NNLO, suppressed strange]}", "lf");
589 
590 
591  // Other alternative samples
592 
593  Int_t plot_index = 0;
594  for (auto& alt_sample : alt_samples) {
595  auto alt_tree = new TChain("tree");
596  files = fileVector(Form("%s/%s/%s", dir.Data(), alt_sample.first.Data(), filePattern.Data()));
597 
598  if (files.size() == 0) {
599  std::cout << "Sample " << alt_sample.first.Data() << " appears to contain no files - check the name of the directory, etc." << std::endl;
600  exit(0);
601  }
602 
603  for (auto file : files)
604  {
605  // Check file first
606  TFile* testfile = TFile::Open(file.c_str());
607  if (testfile != nullptr) {
608  alt_tree->Add(file.c_str());
609  delete testfile;
610  } else {
611  std::cout << "File " << file.c_str() << " is damaged!" << std::endl;
612  }
613  }
614 
615  std::cout << "Generated events in sample " << alt_sample.first.Data() << ": " << alt_tree->GetEntries() << std::endl;
616 
617  TTree* alt_selected = alt_tree->CopyTree(main_preselection.GetTitle());
618  auto yield = DifferentialTaggingYield(alt_selected, draw_config, xvar, "charm");
619 
620  delete alt_selected;
621 
622  auto yield_100fb = static_cast<TH1F *>(yield->Clone(alt_sample.second));
623 
624 
625  for (Int_t i = 1; i <= charm_yield_100fb->GetNbinsX(); i++) {
626  yield_100fb->SetBinError(i, TMath::Sqrt(yield_100fb->GetBinContent(i)));
627  }
628 
629  // Draw the enhanced-suppressed range overlay
630  std::cout << "Draw " << alt_sample.second << std::endl;
631  TH1F *range_band_100fb = static_cast<TH1F *>(ccdis_20Rs2_yield_100fb->Clone(Form("range plot %s", alt_sample.second.Data())));
632  range_band_100fb->Scale(-1.0);
633  range_band_100fb->Add(yield_100fb);
634  range_band_100fb->Divide(ccdis_20Rs2_yield_100fb);
635 
636  for (Int_t i = 1; i <= range_band_100fb->GetNbinsX(); i++) {
637  range_band_100fb->SetBinError(i, 0.0);
638  range_band_100fb->SetBinContent(i, range_band_100fb->GetBinContent(i) + 1.0);
639  }
640  TGraphErrors *range_plot_100fb = new TGraphErrors(range_band_100fb);
641  range_plot_100fb->SetName(Form("range_plot_100fb_%d", plot_index));
642  range_plot_100fb->SetMarkerColor(alt_color[alt_sample.first]);
643  range_plot_100fb->SetLineColor(alt_color[alt_sample.first]);
644  range_plot_100fb->SetMarkerSize(2.0);
645  range_plot_100fb->SetMarkerStyle(alt_marker[alt_sample.first]);
646  range_plot_100fb->SetTitle(alt_sample.second.Data());
647 
648  range_plot_100fb->Draw("E1 SAME P");
649  legend->AddEntry(range_plot_100fb, alt_sample.second.Data(), "lp");
650 
651  plot_index++;
652  }
653 
654 
655  // auto ccdis_21Rs2_data = new TChain("tree");
656  // files = fileVector(Form("%s/CC_DIS_e10_p275_lha_21Rs2/%s", dir.Data(), filePattern.Data()));
657 
658  // for (auto file : files)
659  // {
660  // ccdis_21Rs2_data->Add(file.c_str());
661  // }
662 
663  // auto ccdis_CT18ANNLO_data = new TChain("tree");
664  // files = fileVector(Form("%s/CC_DIS_e10_p275_CT18ANNLO/%s", dir.Data(), filePattern.Data()));
665 
666  // for (auto file : files)
667  // {
668  // ccdis_CT18ANNLO_data->Add(file.c_str());
669  // }
670 
671 
672  // std::cout << "... enhanced strangeness ..." << std::endl;
673  // TTree *ccdis_21Rs2_data_selected = ccdis_21Rs2_data->CopyTree(main_preselection.GetTitle());
674  // std::cout << "... CT18ANNLO ..." << std::endl;
675  // TTree *ccdis_CT18ANNLO_data_selected = ccdis_CT18ANNLO_data->CopyTree(main_preselection.GetTitle());
676 
677  // auto ccdis_21Rs2_yield = DifferentialTaggingYield(ccdis_21Rs2_data_selected, draw_config, xvar, "charm");
678  // auto ccdis_CT18ANNLO_yield = DifferentialTaggingYield(ccdis_CT18ANNLO_data_selected, draw_config, xvar, "charm");
679 
680  // auto ccdis_21Rs2_yield_100fb = static_cast<TH1F *>(ccdis_21Rs2_yield->Clone("ccdis_21Rs2_yield_100fb"));
681  // auto ccdis_CT18ANNLO_yield_100fb = static_cast<TH1F *>(ccdis_CT18ANNLO_yield->Clone("ccdis_CT18ANNLO_yield_100fb"));
682 
683 
684 
685 
686  // // Draw the CT18ANNLO-suppressed range overlay
687  // std::cout << "Draw the Rs range band for CT18ANNLO vs. Suppressed Strangeness..." << std::endl;
688  // TH1F *range_CT18ANNLO_band_100fb = static_cast<TH1F *>(ccdis_20Rs2_yield_100fb->Clone("range_CT18ANNLO_band_100fb"));
689  // range_CT18ANNLO_band_100fb->Scale(-1.0);
690  // range_CT18ANNLO_band_100fb->Add(ccdis_CT18ANNLO_yield_100fb);
691  // range_CT18ANNLO_band_100fb->Divide(ccdis_20Rs2_yield_100fb);
692 
693  // for (Int_t i = 1; i <= range_CT18ANNLO_band_100fb->GetNbinsX(); i++) {
694  // range_CT18ANNLO_band_100fb->SetBinError(i, 0.0);
695  // range_CT18ANNLO_band_100fb->SetBinContent(i, range_CT18ANNLO_band_100fb->GetBinContent(i) + 1.0);
696  // }
697  // TGraphErrors *range_CT18ANNLO_plot_100fb = new TGraphErrors(range_CT18ANNLO_band_100fb);
698  // range_CT18ANNLO_plot_100fb->SetMarkerColor(kGreen - 5);
699  // range_CT18ANNLO_plot_100fb->SetLineColor(kGreen - 5);
700  // range_CT18ANNLO_plot_100fb->SetMarkerSize(2.0);
701  // range_CT18ANNLO_plot_100fb->SetMarkerStyle(kOpenTriangleUp);
702  // range_CT18ANNLO_plot_100fb->SetTitle("CT18ANNLO");
703 
704 
705  // range_plot_100fb->Draw("E1 SAME P");
706  // range_CT18ANNLO_plot_100fb->Draw("E1 SAME P");
707 
708  TLine OneLine(draw_config.xlimits[0], 1.0, draw_config.xlimits[1], 1.0);
709  OneLine.SetLineWidth(2);
710  OneLine.SetLineColor(kBlack);
711  OneLine.Draw("SAME");
712 
713  TLatex plot_title = make_title();
714  plot_title.Draw("SAME");
715 
716 
717  pad->SetLogy(kFALSE);
718 
719  if (xvar.Contains("_x"))
720  pad->SetLogx(kTRUE);
721 
722  pad->SetGrid(1, 1);
723 
724  legend->Draw();
725 
726  pad->Update();
727 
728  pad->SaveAs(Form("CharmTagPlot_tagging_yield_100fb_%s_%s.pdf", input.Data(), xvar.Data()));
729  pad->SaveAs(Form("CharmTagPlot_tagging_yield_100fb_%s_%s.C", input.Data(), xvar.Data()));
730 
731 
732  if (htemplate != nullptr)
733  delete htemplate;
734  }
735 
736 
737  //
738  // Polarized strangeness sensitivity estimate plot
739  //
740 
741  if (xvar == "jet_pt") {
742  draw_config.ylimits = std::vector<float>();
743  draw_config.ylimits.push_back(-100);
744  draw_config.ylimits.push_back(100);
745  draw_config.xtitle = "Reconstructed Jet p_{T} [GeV]";
746  draw_config.ytitle = "Asymmetry Uncertainty [%]";
747  draw_config.logy = kFALSE;
748  draw_config.logx = kFALSE;
749  } else {
750  do_Helicity = kFALSE;
751  }
752 
753  if (do_Helicity == kTRUE) {
754  pad->Clear();
755  gStyle->SetErrorX(0.5);
756 
757  auto charm_yield = DifferentialTaggingYield(default_data_selected, draw_config, xvar, "charm");
758 
759  // auto light_yield = DifferentialTaggingYield(default_data_selected,
760  // draw_config, xvar, "light");
761  Float_t polarization_e = 0.70;
762  Float_t polarization_p = 0.70;
763  Float_t polarization = 0.70; // # single beam polarization
764 
765 
766  // Generate the error band histogram for the suppressed case
767  TH1F *dA_band_100fb = static_cast<TH1F *>(charm_yield->Clone("dA_band_100fb"));
768  Float_t n_gen = static_cast<Float_t>(default_data_selected->GetEntries());
769  std::cout << "N_gen: " << n_gen << std::endl;
770 
771  for (Int_t i = 1; i <= dA_band_100fb->GetNbinsX(); i++) {
772  Float_t n_events = dA_band_100fb->GetBinContent(i);
773  Float_t err_n = TMath::Sqrt(n_events * (1.0 + polarization_e) / 2.0);
774  Float_t dA = 0.0;
775 
776  if (err_n > 0.0) {
777  dA = 1.0 / (err_n * polarization_p) * 100.0;
778  }
779  dA_band_100fb->SetBinError(i, dA);
780  dA_band_100fb->SetBinContent(i, 0.0);
781  }
782 
783  htemplate = new TH1F("htemplate", "", 1, draw_config.xlimits[0], draw_config.xlimits[1]);
784  htemplate->SetXTitle(draw_config.xtitle);
785  htemplate->SetYTitle(draw_config.ytitle);
786  set_axis_range(htemplate, draw_config.ylimits[0], draw_config.ylimits[1], "Y");
787  set_axis_range(htemplate, draw_config.xlimits[0], draw_config.xlimits[1], "X");
788 
789  dA_band_100fb->SetFillColor(kBlue + 2);
790  dA_band_100fb->SetLineColor(kBlue + 2);
791  dA_band_100fb->SetMarkerSize(0.001);
792  dA_band_100fb->SetMarkerColor(kBlue + 2);
793  dA_band_100fb->SetFillColorAlpha(kBlue + 0.2, 0.5);
794 
795  htemplate->Draw("HIST");
796  dA_band_100fb->Draw("E2 SAME");
797 
798  TLine ZeroLine(draw_config.xlimits[0], 0.0, draw_config.xlimits[1], 0.0);
799  ZeroLine.SetLineWidth(2);
800  ZeroLine.SetLineColor(kBlack);
801  ZeroLine.Draw("SAME");
802 
803  TLatex plot_title = make_title();
804  plot_title.Draw("SAME");
805 
806 
807  pad->SetLogy(draw_config.logx);
808  pad->SetLogx(draw_config.logy);
809 
810  pad->SetGrid(1, 1);
811 
812  // Legend
813  TLegend *legend = smart_legend("upper right", 0.25, 0.10);
814  legend->SetFillStyle(0);
815  legend->SetBorderSize(0);
816  legend->AddEntry(dA_band_100fb, "p=70%", "f");
817  legend->Draw();
818 
819  pad->Update();
820 
821  pad->SaveAs(Form("CharmTagPlot_dA_strange_helicity_100fb_%s_%s.pdf", input.Data(), xvar.Data()));
822  pad->SaveAs(Form("CharmTagPlot_dA_strange_helicity_100fb_%s_%s.C", input.Data(), xvar.Data()));
823  }
824 
825 
826  //
827  // Generator-level plots
828  //
829 
830  pad->Clear();
831 
832  bool do_GenJetPlot = kTRUE;
833 
834  if (xvar == "GenJet.PT") {
835  float xbins[] = { 0, 2.5, 5.0, 7.5, 10, 12.5, 15, 20, 25, 35, 60 };
836  int nbins = sizeof(xbins) / sizeof(xbins[0]) - 1;
837  draw_config.htemplate = new TH1F(xvar, "", nbins, xbins);
838 
839  draw_config.xlimits = std::vector<float>();
840  draw_config.xlimits.push_back(0);
841  draw_config.xlimits.push_back(60);
842 
843  draw_config.ylimits = std::vector<float>();
844  draw_config.ylimits.push_back(1);
845  draw_config.ylimits.push_back(2e5);
846 
847  draw_config.logy = kTRUE;
848  draw_config.logx = kFALSE;
849 
850  draw_config.xtitle = "Generated Jet p_{T} [GeV]";
851  draw_config.ytitle = "d#sigma/dp_{T} #times 100fb^{-1} [GeV^{-1}]";
852  } else {
853  do_GenJetPlot = kFALSE;
854  }
855 
856 
857  if (do_GenJetPlot == kTRUE) {
858  pad = new TCanvas("pad", "", 900, 1200);
859 
860  // Load Delphes samples for this plot
861  auto delphes_data = new TChain("Delphes");
862  files = fileVector(Form("%s/../%s/%s", dir.Data(), input.Data(), filePattern.Data()));
863 
864  for (auto file : files)
865  {
866  delphes_data->Add(file.c_str());
867  }
868 
869  TH1F *light_genjets = static_cast<TH1F *>(draw_config.htemplate->Clone(Form("GenJet_light_%s_%d", xvar.Data(), getHistUID())));
870  light_genjets->Sumw2();
871 
872  TH1F *charm_genjets = static_cast<TH1F *>(draw_config.htemplate->Clone(Form("GenJet_charm_%s_%d", xvar.Data(), getHistUID())));
873  charm_genjets->Sumw2();
874 
875  delphes_data->Project(light_genjets->GetName(), xvar.Data(), "");
876  delphes_data->Project(charm_genjets->GetName(), xvar.Data(), "GenJet.Flavor == 4");
877 
878  // Transform these into differential cross-section distributions
879  Float_t n_gen = static_cast<float>(delphes_data->GetEntries());
880 
881  for (Int_t ibin = 1; ibin <= light_genjets->GetNbinsX(); ibin++) {
882  Float_t dydx = light_genjets->GetBinContent(ibin) * target_lumi * target_xsec / n_gen / light_genjets->GetBinWidth(ibin);
883  Float_t dydx_relerr = light_genjets->GetBinError(ibin) / light_genjets->GetBinContent(ibin);
884  light_genjets->SetBinContent(ibin, dydx);
885  light_genjets->SetBinError(ibin, dydx * dydx_relerr);
886 
887  dydx = charm_genjets->GetBinContent(ibin) * target_lumi * target_xsec / n_gen / charm_genjets->GetBinWidth(ibin);
888  dydx_relerr = charm_genjets->GetBinError(ibin) / charm_genjets->GetBinContent(ibin);
889  charm_genjets->SetBinContent(ibin, dydx);
890  charm_genjets->SetBinError(ibin, dydx * dydx_relerr);
891  }
892 
893  configure_plot<TH1F>(charm_genjets, draw_config, "charm");
894  configure_plot<TH1F>(light_genjets, draw_config, "light");
895 
896  // generate template histogram
897  htemplate = new TH1F("htemplate", "", 1, draw_config.xlimits[0], draw_config.xlimits[1]);
898 
899  set_axis_range(htemplate, draw_config.ylimits[0], draw_config.ylimits[1], "Y");
900  set_axis_range(htemplate, draw_config.xlimits[0], draw_config.xlimits[1], "X");
901  htemplate->SetXTitle(draw_config.xtitle);
902  htemplate->SetYTitle(draw_config.ytitle);
903 
904  set_axis_range(charm_genjets, draw_config.ylimits[0], draw_config.ylimits[1], "Y");
905  set_axis_range(charm_genjets, draw_config.xlimits[0], draw_config.xlimits[1], "X");
906  charm_genjets->SetXTitle(draw_config.xtitle);
907  charm_genjets->SetYTitle(draw_config.ytitle);
908 
909  set_axis_range(light_genjets, draw_config.ylimits[0], draw_config.ylimits[1], "Y");
910  set_axis_range(light_genjets, draw_config.xlimits[0], draw_config.xlimits[1], "X");
911  light_genjets->SetXTitle(draw_config.xtitle);
912  light_genjets->SetYTitle(draw_config.ytitle);
913 
914 
915  pad->cd();
916  pad->SetLogy(draw_config.logy);
917  pad->SetLogx(draw_config.logx);
918  pad->SetGrid(1, 1);
919  pad->Update();
920 
921  auto rp = new TRatioPlot(charm_genjets, light_genjets);
922 
923  rp->SetH1DrawOpt("E1 P");
924  rp->SetH2DrawOpt("E1 P");
925 
926  rp->Draw();
927  rp->GetLowerRefYaxis()->SetTitle("Charm-to-All Ratio");
928  rp->GetLowerRefYaxis()->SetLabelSize(0.022);
929  rp->GetLowerRefYaxis()->SetTitleSize(0.027);
930  rp->GetLowerRefYaxis()->SetTitleOffset(2);
931  rp->SetLeftMargin(0.15);
932  rp->SetRightMargin(0.05);
933  rp->SetLowBottomMargin(0.55);
934 
935  rp->GetLowerRefYaxis()->SetRangeUser(0.0, 0.05);
936  rp->GetLowerRefYaxis()->SetNdivisions(5, 1, 0, kFALSE);
937  rp->GetLowYaxis()->SetRangeUser(-0.01, 0.051);
938  rp->GetLowYaxis()->SetNdivisions(5, 1, 0, kTRUE);
939 
940  TLatex plot_title = make_title();
941  plot_title.Draw("SAME");
942 
943 
944  // Legend
945  rp->GetUpperPad()->cd();
946  TLegend *legend = smart_legend("lower left", 0.55, 0.15);
947  legend->SetFillStyle(0);
948  legend->SetBorderSize(0);
949  legend->AddEntry(light_genjets, "All Jets [CT18NNLO]", "lp");
950  legend->AddEntry(charm_genjets, "Charm Jets [CT18NNLO]", "lp");
951  legend->Draw();
952 
953  rp->GetUpperPad()->SetGrid(1, 1);
954 
955  pad->Update();
956  pad->SaveAs(Form("CharmTagPlot_differential_xs_%s_%s.pdf", delphes_data->GetTitle(), xvar.Data()));
957  pad->SaveAs(Form("CharmTagPlot_differential_xs_%s_%s.C", delphes_data->GetTitle(), xvar.Data()));
958 
959  // Cleanup
960  if (htemplate != nullptr)
961  delete htemplate;
962 
963  if (legend != nullptr)
964  delete legend;
965  }
966 }