EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CharmJetMVAOptimization.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CharmJetMVAOptimization.C
1 #include "TROOT.h"
2 #include "TChain.h"
3 #include "TFile.h"
4 #include "TEfficiency.h"
5 #include "TH1.h"
6 #include "TGraphErrors.h"
7 #include "TCut.h"
8 #include "TCanvas.h"
9 #include "TStyle.h"
10 #include "TLegend.h"
11 #include "TMath.h"
12 #include "TLine.h"
13 #include "TLatex.h"
14 #include "TRatioPlot.h"
15 
16 
17 #include <glob.h>
18 #include <iostream>
19 #include <iomanip>
20 #include <vector>
21 
22 #include "PlotFunctions.h"
23 
24 
25 void CharmJetMVAOptimization(TString filename = "CharmJetClassification_Results.root",
26  TString foldername = "dataset",
27  TString classifiername = "MLP")
28 {
29  // Global options
30  gStyle->SetOptStat(0);
31 
32 
33  // Scan the output for the optimal performance point and quote the results
34  float target_xsec = LookupCrossSection("CC_DIS");
35  float eff_xsec_light = target_xsec * 0.983;
36  float eff_xsec_charm = target_xsec * 0.020;
37  float target_lumi = 100*u_fb;
38 
39  auto result_file = TFile::Open(filename.Data());
40  auto result_tree = static_cast<TTree*>(result_file->Get(Form("%s/TestTree", foldername.Data())));
41 
42  float max_charm_yield_signifiance = 0.0;
43  float optimal_mva_cut = 0.0;
44  float optimal_charm_efficiency = 0.0;
45  float optimal_light_efficiency = 0.0;
46  float optimal_charm_yield = 0.0;
47  float optimal_light_yield = 0.0;
48 
49 
50  float n_light_all = static_cast<float>(result_tree->GetEntries("(jet_flavor==21 || jet_flavor < 4)"));
51  float n_charm_all = static_cast<float>(result_tree->GetEntries("(jet_flavor == 4)"));
52 
53  float w_light = target_lumi * eff_xsec_light / n_light_all;
54  float w_charm = target_lumi * eff_xsec_charm / n_charm_all;
55 
56  TH1F* h_light = new TH1F("h_light","",200,0,1.0);
57  h_light->Sumw2();
58  auto h_charm = static_cast<TH1F*>(h_light->Clone("h_charm"));
59 
60  result_tree->Project(h_light->GetName(), classifiername.Data(), "jet_pt>5.0 && TMath::Abs(jet_eta) < 3.0 && (jet_flavor==21 || jet_flavor < 4) && met_et > 10");
61  result_tree->Project(h_charm->GetName(), classifiername.Data(), "jet_pt>5.0 && TMath::Abs(jet_eta) < 3.0 && (jet_flavor==4) && met_et > 10");
62 
63 
64  Int_t nbins = h_light->GetNbinsX();
65 
66  for (Int_t i = 1; i <= h_light->GetNbinsX(); i++) {
67  float min = h_light->GetBinLowEdge(i);
68  std::cout << min << std::endl;
69 
70  float n_light = h_light->Integral(i, nbins);
71  float n_charm = h_charm->Integral(i, nbins);
72 
73  float eff_light = n_light/n_light_all;
74  float eff_charm = n_charm/n_charm_all;
75 
76 
77  // Scale the number of light and charm yields to target_lumi
78  n_light = n_light * w_light;
79  n_charm = n_charm * w_charm;
80  float n_total = n_light + n_charm;
81 
82  float err_light = TMath::Sqrt(n_light);
83  float err_charm = TMath::Sqrt(n_charm);
84 
85 
86  float err_yield = TMath::Sqrt(n_total + n_light);
87  float yield_significance = n_charm / err_yield;
88 
89  if (yield_significance > max_charm_yield_signifiance) {
90  max_charm_yield_signifiance = yield_significance;
91  optimal_mva_cut = min;
92  optimal_charm_efficiency = eff_charm;
93  optimal_light_efficiency = eff_light;
94  optimal_charm_yield = n_charm;
95  optimal_light_yield = n_light;
96  }
97 
98  }
99 
100 
101  // for (float min = 0.0; min <= 1.0; min += 0.02) {
102  // std::cout << min << std::endl;
103 
104  // float n_light = static_cast<float>(result_tree->GetEntries(Form("jet_pt>5.0 && TMath::Abs(jet_eta) < 3.0 && (jet_flavor==21 || jet_flavor < 4) && met_et > 10 && MLP > %.5f", min)));
105  // float n_charm = static_cast<float>(result_tree->GetEntries(Form("jet_pt>5.0 && TMath::Abs(jet_eta) < 3.0 && (jet_flavor == 4) && met_et > 10 && MLP > %.5f", min)));
106 
107 
108  // float eff_light = n_light/n_light_all;
109  // float eff_charm = n_charm/n_charm_all;
110 
111  // // Scale the number of light and charm yields to target_lumi
112  // n_light = n_light * w_light;
113  // n_charm = n_charm * w_charm;
114  // float n_total = n_light + n_charm;
115 
116  // float err_light = TMath::Sqrt(n_light);
117  // float err_charm = TMath::Sqrt(n_charm);
118 
119 
120  // float err_yield = TMath::Sqrt(n_total + n_light);
121  // float yield_significance = n_charm / err_yield;
122 
123  // if (yield_significance > max_charm_yield_signifiance) {
124  // max_charm_yield_signifiance = yield_significance;
125  // optimal_mva_cut = min;
126  // optimal_charm_efficiency = eff_charm;
127  // optimal_light_efficiency = eff_light;
128  // }
129 
130 
131  // }
132 
133  std::cout << "Optimal MVA cut: " << optimal_mva_cut
134  << " with charm yield significance of " << max_charm_yield_signifiance
135  << " charm (light) efficiency " << optimal_charm_efficiency << " (" << optimal_light_efficiency << ")"
136  << " in 100/fb" << std::endl;
137  std::cout << "Estimated charm (light) yield: " << optimal_charm_yield << "(" << optimal_light_yield << ")" << std::endl;
138 
139 }