EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EventSelectionModule.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EventSelectionModule.cc
1 #include "EventSelectionModule.h"
2 
3 #include "TClonesArray.h"
4 #include "TRandom3.h"
5 #include "TDatabasePDG.h"
6 #include "Math/PdfFuncMathCore.h"
7 
8 #include "AnalysisFunctions.cc"
9 #include "TreeHandler.h"
10 #include "TSystem.h"
11 
12 #include <iostream>
13 #include <iomanip>
14 #include <fstream>
15 #include <algorithm>
16 
17 EventSelectionModule::EventSelectionModule(ExRootTreeReader *data)
18  : Module(data)
19 {}
20 
21 EventSelectionModule::~EventSelectionModule()
22 {}
23 
24 void EventSelectionModule::initialize()
25 {
26  // Get the pointer to the mRICH tracks
27  _branch_mRICHTracks = getData()->UseBranch("mRICHTrack");
28  _branch_barrelDircTracks = getData()->UseBranch("barrelDIRCTrack");
29  _branch_dualRICHagTracks = getData()->UseBranch("dualRICHagTrack");
30  _branch_dualRICHcfTracks = getData()->UseBranch("dualRICHcfTrack");
31 
32 
33  // Tree handler initialization
34  TreeHandler *tree_handler = tree_handler->getInstance();
35 
36  if (tree_handler->getTree() != nullptr) {
37  _jet_n = 0;
38  _jet_pt = std::vector<float>();
39  _jet_eta = std::vector<float>();
40  _jet_flavor = std::vector<int>();
41  _jet_nconstituents = std::vector<int>();
42  _jet_sip3dtag = std::vector<int>();
43 
44  _jet_ktag = std::vector<int>();
45  _jet_k1_pt = std::vector<float>();
46  _jet_k1_sIP3D = std::vector<float>();
47  _jet_k1_pid = std::vector<int>();
48  _jet_k2_pt = std::vector<float>();
49  _jet_k2_sIP3D = std::vector<float>();
50  _jet_k2_pid = std::vector<int>();
51 
52  _jet_e1_pt = std::vector<float>();
53  _jet_e1_sIP3D = std::vector<float>();
54  _jet_e2_pt = std::vector<float>();
55  _jet_e2_sIP3D = std::vector<float>();
56 
57  _jet_mu1_pt = std::vector<float>();
58  _jet_mu1_sIP3D = std::vector<float>();
59  _jet_mu2_pt = std::vector<float>();
60  _jet_mu2_sIP3D = std::vector<float>();
61 
62  _jet_t1_pt = std::vector<float>();
63  _jet_t1_sIP3D = std::vector<float>();
64  _jet_t2_pt = std::vector<float>();
65  _jet_t2_sIP3D = std::vector<float>();
66  _jet_t3_pt = std::vector<float>();
67  _jet_t3_sIP3D = std::vector<float>();
68  _jet_t4_pt = std::vector<float>();
69  _jet_t4_sIP3D = std::vector<float>();
70 
71 
72  _jet_etag = std::vector<int>();
73  _jet_mutag = std::vector<int>();
74 
75  _jet_mlp_ktagger = std::vector<float>();
76  _jet_mlp_eltagger = std::vector<float>();
77  _jet_mlp_mutagger = std::vector<float>();
78  _jet_mlp_ip3dtagger = std::vector<float>();
79 
80 
81  _jet_charge = std::vector<float>();
82 
83 
84  tree_handler->getTree()->Branch("jet_n", &_jet_n, "jet_n/I");
85  tree_handler->getTree()->Branch("jet_pt", "std::vector<float>", &_jet_pt);
86  tree_handler->getTree()->Branch("jet_eta", "std::vector<float>", &_jet_eta);
87  tree_handler->getTree()->Branch("jet_flavor", "std::vector<int>", &_jet_flavor);
88  tree_handler->getTree()->Branch("jet_nconstituents", "std::vector<int>", &_jet_nconstituents);
89  tree_handler->getTree()->Branch("jet_sip3dtag", "std::vector<int>", &_jet_sip3dtag);
90 
91  tree_handler->getTree()->Branch("jet_ktag", "std::vector<int>", &_jet_ktag);
92  tree_handler->getTree()->Branch("jet_k1_pt", "std::vector<float>", &_jet_k1_pt);
93  tree_handler->getTree()->Branch("jet_k1_pid", "std::vector<int>", &_jet_k1_pid);
94  tree_handler->getTree()->Branch("jet_k1_sIP3D", "std::vector<float>", &_jet_k1_sIP3D);
95  tree_handler->getTree()->Branch("jet_k2_pt", "std::vector<float>", &_jet_k2_pt);
96  tree_handler->getTree()->Branch("jet_k2_pid", "std::vector<int>", &_jet_k2_pid);
97  tree_handler->getTree()->Branch("jet_k2_sIP3D", "std::vector<float>", &_jet_k2_sIP3D);
98 
99  tree_handler->getTree()->Branch("jet_e1_pt", "std::vector<float>", &_jet_e1_pt);
100  tree_handler->getTree()->Branch("jet_e1_sIP3D", "std::vector<float>", &_jet_e1_sIP3D);
101  tree_handler->getTree()->Branch("jet_e2_pt", "std::vector<float>", &_jet_e2_pt);
102  tree_handler->getTree()->Branch("jet_e2_sIP3D", "std::vector<float>", &_jet_e2_sIP3D);
103 
104  tree_handler->getTree()->Branch("jet_mu1_pt", "std::vector<float>", &_jet_mu1_pt);
105  tree_handler->getTree()->Branch("jet_mu1_sIP3D", "std::vector<float>", &_jet_mu1_sIP3D);
106  tree_handler->getTree()->Branch("jet_mu2_pt", "std::vector<float>", &_jet_mu2_pt);
107  tree_handler->getTree()->Branch("jet_mu2_sIP3D", "std::vector<float>", &_jet_mu2_sIP3D);
108 
109  tree_handler->getTree()->Branch("jet_t1_sIP3D", "std::vector<float>", &_jet_t1_sIP3D);
110  tree_handler->getTree()->Branch("jet_t2_sIP3D", "std::vector<float>", &_jet_t2_sIP3D);
111  tree_handler->getTree()->Branch("jet_t3_sIP3D", "std::vector<float>", &_jet_t3_sIP3D);
112  tree_handler->getTree()->Branch("jet_t4_sIP3D", "std::vector<float>", &_jet_t4_sIP3D);
113 
114  tree_handler->getTree()->Branch("jet_t1_pt", "std::vector<float>", &_jet_t1_pt);
115  tree_handler->getTree()->Branch("jet_t2_pt", "std::vector<float>", &_jet_t2_pt);
116  tree_handler->getTree()->Branch("jet_t3_pt", "std::vector<float>", &_jet_t3_pt);
117  tree_handler->getTree()->Branch("jet_t4_pt", "std::vector<float>", &_jet_t4_pt);
118 
119  tree_handler->getTree()->Branch("jet_mlp_ip3dtagger", "std::vector<float>", &_jet_mlp_ip3dtagger);
120  tree_handler->getTree()->Branch("jet_mlp_ktagger", "std::vector<float>", &_jet_mlp_ktagger);
121  tree_handler->getTree()->Branch("jet_mlp_eltagger", "std::vector<float>", &_jet_mlp_eltagger);
122  tree_handler->getTree()->Branch("jet_mlp_mutagger", "std::vector<float>", &_jet_mlp_mutagger);
123  tree_handler->getTree()->Branch("jet_mlp_globaltagger", "std::vector<float>", &_jet_mlp_globaltagger);
124 
125  tree_handler->getTree()->Branch("jet_etag", "std::vector<int>", &_jet_etag);
126  tree_handler->getTree()->Branch("jet_mutag", "std::vector<int>", &_jet_mutag);
127  tree_handler->getTree()->Branch("jet_charge", "std::vector<float>", &_jet_charge);
128 
129  _jet_Ks_mass = std::vector<std::vector<float> >();
130  _jet_Ks_p = std::vector<std::vector<float> >();
131  _jet_Ks_flightlength = std::vector<std::vector<float> >();
132  _jet_Ks_sumpt = std::vector<float>();
133  _jet_K_sumpt = std::vector<float>();
134  _jet_Ks_zhadron = std::vector<std::vector<float> >();
135  _jet_K_zhadron = std::vector<std::vector<float> >();
136  _jet_Ks_leading_zhadron = std::vector<float>();
137  _jet_K_leading_zhadron = std::vector<float>();
138  _jet_ehadoveremratio = std::vector<float>();
139 
140  tree_handler->getTree()->Branch("jet_Ks_mass", "std::vector<std::vector<float>>", &_jet_Ks_mass);
141  tree_handler->getTree()->Branch("jet_Ks_p", "std::vector<std::vector<float>>", &_jet_Ks_p);
142  tree_handler->getTree()->Branch("jet_Ks_flightlength", "std::vector<std::vector<float>>", &_jet_Ks_flightlength);
143  tree_handler->getTree()->Branch("jet_Ks_sumpt", "std::vector<float>", &_jet_Ks_sumpt);
144  tree_handler->getTree()->Branch("jet_K_sumpt", "std::vector<float>", &_jet_K_sumpt);
145  tree_handler->getTree()->Branch("jet_Ks_zhadron", "std::vector<std::vector<float>>", &_jet_Ks_zhadron);
146  tree_handler->getTree()->Branch("jet_K_zhadron", "std::vector<std::vector<float>>", &_jet_K_zhadron);
147  tree_handler->getTree()->Branch("jet_Ks_leading_zhadron", "std::vector<float>", &_jet_Ks_leading_zhadron);
148  tree_handler->getTree()->Branch("jet_K_leading_zhadron", "std::vector<float>", &_jet_K_leading_zhadron);
149  tree_handler->getTree()->Branch("jet_ehadoveremratio", "std::vector<float>", &_jet_ehadoveremratio);
150 
151  _charmjet_n = 0;
152  _charmjet_pt = std::vector<float>();
153  _charmjet_eta = std::vector<float>();
154  tree_handler->getTree()->Branch("charmjet_n", &_charmjet_n, "charmjet_n/I");
155  tree_handler->getTree()->Branch("charmjet_pt", "std::vector<float>", &_charmjet_pt);
156  tree_handler->getTree()->Branch("charmjet_eta", "std::vector<float>", &_charmjet_eta);
157 
158  _met_et = 0.0;
159  tree_handler->getTree()->Branch("met_et", &_met_et, "met_et/F");
160 
161  _bjorken_x = 0.0;
162  _bjorken_Q2 = 0.0;
163  _bjorken_y = 0.0;
164  _jb_x = 0.0;
165  _jb_Q2 = 0.0;
166  tree_handler->getTree()->Branch("bjorken_x", &_bjorken_x, "bjorken_x/F");
167  tree_handler->getTree()->Branch("bjorken_Q2", &_bjorken_Q2, "bjorken_Q2/F");
168  tree_handler->getTree()->Branch("bjorken_y", &_bjorken_y, "bjorken_y/F");
169  tree_handler->getTree()->Branch("jb_x", &_jb_x, "jb_x/F");
170  tree_handler->getTree()->Branch("jb_Q2", &_jb_Q2, "jb_Q2/F");
171 
172  // pid tracks
173  _pid_track_n = 0;
174  _pid_track_jetmother = std::vector<int>();
175  _pid_track_jet_pt = std::vector<float>();
176  _pid_track_jet_eta = std::vector<float>();
177  _pid_track_true_pid = std::vector<int>();
178  _pid_track_reco_pid = std::vector<int>();
179  _pid_track_eta = std::vector<float>();
180  _pid_track_pt = std::vector<float>();
181 
182  tree_handler->getTree()->Branch("pid_track_n", &_pid_track_n, "pid_track_n/I");
183  tree_handler->getTree()->Branch("pid_track_true_pid", "std::vector<int>", &_pid_track_true_pid);
184  tree_handler->getTree()->Branch("pid_track_reco_pid", "std::vector<int>", &_pid_track_reco_pid);
185  tree_handler->getTree()->Branch("pid_track_pt", "std::vector<float>", &_pid_track_pt);
186  tree_handler->getTree()->Branch("pid_track_eta", "std::vector<float>", &_pid_track_eta);
187  tree_handler->getTree()->Branch("pid_track_jetmother", "std::vector<int>", &_pid_track_jetmother);
188  tree_handler->getTree()->Branch("pid_track_jet_pt", "std::vector<float>", &_pid_track_jet_pt);
189  tree_handler->getTree()->Branch("pid_track_jet_eta", "std::vector<float>", &_pid_track_jet_eta);
190  }
191 
192 
193  // Initialize the cut flow
194 
195  _cut_flow["1: All events"] = 0;
196  _cut_flow["2: MET > 10 GeV"] = 0;
197  _cut_flow["3: Fiducial Jets >= 1"] = 0;
198  _cut_flow["4: Charm Jet == 1"] = 0;
199 
200  // Global variables
201  _mpi = TDatabasePDG().GetParticle(211)->Mass();
202  _mK = TDatabasePDG().GetParticle(321)->Mass();
203 
204 
205  // MVA Algorithms
206  _mva_reader_ip3dtagger = new TMVA::Reader("Silent");
207  _mva_reader_ktagger = new TMVA::Reader("Silent");
208  _mva_reader_eltagger = new TMVA::Reader("Silent");
209  _mva_reader_mutagger = new TMVA::Reader("Silent");
210  _mva_reader_globaltagger = new TMVA::Reader("Silent");
211 
212 
213  _mva_inputs_float = std::map<TString, Float_t>();
214  _mva_inputs_int = std::map<TString, Int_t>();
215  _mva_inputs_float["jet_pt"] = 0.0;
216  _mva_inputs_float["jet_eta"] = 0.0;
217  _mva_inputs_int["jet_flavor"] = 0;
218  _mva_inputs_float["met_et"] = 0.0;
219 
220 
221  _mva_inputs_float["jet_e1_pt"] = 0.0;
222  _mva_inputs_float["jet_e1_sIP3D"] = 0.0;
223  _mva_inputs_float["jet_e2_pt"] = 0.0;
224  _mva_inputs_float["jet_e2_sIP3D"] = 0.0;
225 
226  _mva_inputs_float["jet_t1_pt"] = 0.0;
227  _mva_inputs_float["jet_t1_sIP3D"] = 0.0;
228  _mva_inputs_float["jet_t2_pt"] = 0.0;
229  _mva_inputs_float["jet_t2_sIP3D"] = 0.0;
230  _mva_inputs_float["jet_t3_pt"] = 0.0;
231  _mva_inputs_float["jet_t3_sIP3D"] = 0.0;
232  _mva_inputs_float["jet_t4_pt"] = 0.0;
233  _mva_inputs_float["jet_t4_sIP3D"] = 0.0;
234 
235  _mva_inputs_float["jet_k1_pt"] = 0.0;
236  _mva_inputs_float["jet_k1_sIP3D"] = 0.0;
237  _mva_inputs_float["jet_k1_pid"] = 0;
238  _mva_inputs_float["jet_k2_pt"] = 0.0;
239  _mva_inputs_float["jet_k2_sIP3D"] = 0.0;
240  _mva_inputs_float["jet_k2_pid"] = 0;
241 
242  _mva_inputs_float["jet_mu1_pt"] = 0.0;
243  _mva_inputs_float["jet_mu1_sIP3D"] = 0.0;
244  _mva_inputs_float["jet_mu2_pt"] = 0.0;
245  _mva_inputs_float["jet_mu2_sIP3D"] = 0.0;
246 
247 
248  _mva_inputs_float["jet_mlp_ip3dtagger"] = 0.0;
249  _mva_inputs_float["jet_mlp_ktagger"] = 0.0;
250  _mva_inputs_float["jet_mlp_eltagger"] = 0.0;
251  _mva_inputs_float["jet_mlp_mutagger"] = 0.0;
252 
253 
254  _mva_reader_ip3dtagger->AddSpectator("jet_pt", &(_mva_inputs_float["jet_pt"]));
255  _mva_reader_ip3dtagger->AddSpectator("jet_eta", &(_mva_inputs_float["jet_eta"]));
256  _mva_reader_ip3dtagger->AddSpectator("jet_flavor", &(_mva_inputs_int["jet_flavor"]));
257  _mva_reader_ip3dtagger->AddSpectator("met_et", &(_mva_inputs_float["met_et"]));
258 
259  _mva_reader_ktagger->AddSpectator("jet_pt", &(_mva_inputs_float["jet_pt"]));
260  _mva_reader_ktagger->AddSpectator("jet_eta", &(_mva_inputs_float["jet_eta"]));
261  _mva_reader_ktagger->AddSpectator("jet_flavor", &(_mva_inputs_int["jet_flavor"]));
262  _mva_reader_ktagger->AddSpectator("met_et", &(_mva_inputs_float["met_et"]));
263 
264  _mva_reader_eltagger->AddSpectator("jet_pt", &(_mva_inputs_float["jet_pt"]));
265  _mva_reader_eltagger->AddSpectator("jet_eta", &(_mva_inputs_float["jet_eta"]));
266  _mva_reader_eltagger->AddSpectator("jet_flavor", &(_mva_inputs_int["jet_flavor"]));
267  _mva_reader_eltagger->AddSpectator("met_et", &(_mva_inputs_float["met_et"]));
268 
269  _mva_reader_mutagger->AddSpectator("jet_pt", &(_mva_inputs_float["jet_pt"]));
270  _mva_reader_mutagger->AddSpectator("jet_eta", &(_mva_inputs_float["jet_eta"]));
271  _mva_reader_mutagger->AddSpectator("jet_flavor", &(_mva_inputs_int["jet_flavor"]));
272  _mva_reader_mutagger->AddSpectator("met_et", &(_mva_inputs_float["met_et"]));
273 
274  _mva_reader_globaltagger->AddSpectator("jet_pt", &(_mva_inputs_float["jet_pt"]));
275  _mva_reader_globaltagger->AddSpectator("jet_eta", &(_mva_inputs_float["jet_eta"]));
276  _mva_reader_globaltagger->AddSpectator("jet_flavor", &(_mva_inputs_int["jet_flavor"]));
277  _mva_reader_globaltagger->AddSpectator("met_et", &(_mva_inputs_float["met_et"]));
278 
279 
280  _mva_reader_eltagger->AddVariable("jet_e1_pt", &(_mva_inputs_float["jet_e1_pt"]));
281  _mva_reader_eltagger->AddVariable("jet_e1_sIP3D", &(_mva_inputs_float["jet_e1_sIP3D"]));
282  _mva_reader_eltagger->AddVariable("jet_e2_pt", &(_mva_inputs_float["jet_e2_pt"]));
283  _mva_reader_eltagger->AddVariable("jet_e2_sIP3D", &(_mva_inputs_float["jet_e2_sIP3D"]));
284 
285  _mva_reader_ip3dtagger->AddVariable("jet_t1_pt", &(_mva_inputs_float["jet_t1_pt"]));
286  _mva_reader_ip3dtagger->AddVariable("jet_t1_sIP3D", &(_mva_inputs_float["jet_t1_sIP3D"]));
287  _mva_reader_ip3dtagger->AddVariable("jet_t2_pt", &(_mva_inputs_float["jet_t2_pt"]));
288  _mva_reader_ip3dtagger->AddVariable("jet_t2_sIP3D", &(_mva_inputs_float["jet_t2_sIP3D"]));
289  _mva_reader_ip3dtagger->AddVariable("jet_t3_pt", &(_mva_inputs_float["jet_t3_pt"]));
290  _mva_reader_ip3dtagger->AddVariable("jet_t3_sIP3D", &(_mva_inputs_float["jet_t3_sIP3D"]));
291  _mva_reader_ip3dtagger->AddVariable("jet_t4_pt", &(_mva_inputs_float["jet_t4_pt"]));
292  _mva_reader_ip3dtagger->AddVariable("jet_t4_sIP3D", &(_mva_inputs_float["jet_t4_sIP3D"]));
293 
294  _mva_reader_ktagger->AddVariable("jet_k1_pt", &(_mva_inputs_float["jet_k1_pt"]));
295  _mva_reader_ktagger->AddVariable("jet_k1_sIP3D", &(_mva_inputs_float["jet_k1_sIP3D"]));
296  _mva_reader_ktagger->AddVariable("jet_k2_pt", &(_mva_inputs_float["jet_k2_pt"]));
297  _mva_reader_ktagger->AddVariable("jet_k2_sIP3D", &(_mva_inputs_float["jet_k2_sIP3D"]));
298 
299  _mva_reader_mutagger->AddVariable("jet_mu1_pt", &(_mva_inputs_float["jet_mu1_pt"]));
300  _mva_reader_mutagger->AddVariable("jet_mu1_sIP3D", &(_mva_inputs_float["jet_mu1_sIP3D"]));
301  _mva_reader_mutagger->AddVariable("jet_mu2_pt", &(_mva_inputs_float["jet_mu2_pt"]));
302  _mva_reader_mutagger->AddVariable("jet_mu2_sIP3D", &(_mva_inputs_float["jet_mu2_sIP3D"]));
303 
304  _mva_reader_globaltagger->AddVariable("jet_mlp_ip3dtagger", &(_mva_inputs_float["jet_mlp_ip3dtagger"]));
305  _mva_reader_globaltagger->AddVariable("jet_mlp_ktagger", &(_mva_inputs_float["jet_mlp_ktagger"]));
306  _mva_reader_globaltagger->AddVariable("jet_mlp_eltagger", &(_mva_inputs_float["jet_mlp_eltagger"]));
307  _mva_reader_globaltagger->AddVariable("jet_mlp_mutagger", &(_mva_inputs_float["jet_mlp_mutagger"]));
308 
309  _mva_reader_ip3dtagger->BookMVA("CharmIP3DTagger", "mva_taggers/TMVAClassification_CharmIP3DTagger.weights.xml");
310  _mva_reader_eltagger->BookMVA("CharmETagger", "mva_taggers/TMVAClassification_CharmETagger.weights.xml");
311  _mva_reader_ktagger->BookMVA("CharmKTagger", "mva_taggers/TMVAClassification_CharmKTagger.weights.xml");
312  _mva_reader_mutagger->BookMVA("CharmMuTagger", "mva_taggers/TMVAClassification_CharmMuTagger.weights.xml");
313  _mva_reader_globaltagger->BookMVA("CharmGlobalTagger", "mva_taggers/TMVAClassification_CharmGlobalTagger.weights.xml");
314 }
315 
316 void EventSelectionModule::finalize()
317 {
318  ofstream csvfile;
319 
320  csvfile.open("cut_flow.csv");
321 
322  csvfile << "Cut,Yield" << std::endl;
323 
324  for (auto&[cut, yield] : _cut_flow) {
325  csvfile << "\"" << cut << "\"," << int(yield) << std::endl;
326  }
327 
328  csvfile.close();
329 
330  if (_mva_reader_ip3dtagger) {
331  delete _mva_reader_ip3dtagger;
332  }
333 
334  if (_mva_reader_ktagger) {
335  delete _mva_reader_ktagger;
336  }
337 
338  if (_mva_reader_eltagger) {
339  delete _mva_reader_eltagger;
340  }
341 
342  if (_mva_reader_mutagger) {
343  delete _mva_reader_mutagger;
344  }
345 
346  if (_mva_reader_globaltagger) {
347  delete _mva_reader_globaltagger;
348  }
349 }
350 
351 bool EventSelectionModule::execute(std::map<std::string, std::any> *DataStore)
352 {
353  auto data = getData();
354 
355 
356  // Compute global DIS variables
357  auto dis_variables = DISVariables(getGenParticles());
358 
359  _bjorken_x = dis_variables["x"];
360  _bjorken_Q2 = dis_variables["Q2"];
361  _bjorken_y = dis_variables["y"];
362 
363  auto jb_variables = DISJacquetBlondel(getEFlowTracks(), getElectrons(), getPhotons(), getNeutralHadrons());
364  _jb_x = jb_variables["x_JB"];
365  _jb_Q2 = jb_variables["Q2_JB"];
366 
367 
368  // Initialize output variables
369  // _charmjet_pt = std::vector<float>();
370  // _charmjet_eta = std::vector<float>();
371  _charmjet_pt.clear();
372  _charmjet_eta.clear();
373  _charmjet_n = _charmjet_pt.size();
374 
375  _pid_track_jetmother.clear();
376  _pid_track_jet_pt.clear();
377  _pid_track_jet_eta.clear();
378  _pid_track_true_pid.clear();
379  _pid_track_reco_pid.clear();
380  _pid_track_pt.clear();
381  _pid_track_eta.clear();
382  _pid_track_n = _pid_track_pt.size();
383 
384 
385  // Cut flow
386  _cut_flow["1: All events"] += 1;
387  bool passed = true;
388 
389  // Get the MET object and use it
390  MissingET *MET = nullptr;
391 
392  for (int imet = 0; imet < getMET()->GetEntries(); imet++) {
393  MET = static_cast<MissingET *>(getMET()->At(imet));
394  }
395 
396  if (MET == nullptr) {
397  passed = false;
398  }
399 
400  _met_et = MET->MET;
401 
402  if ((passed == true) && (MET->MET > 10.0)) {
403  _cut_flow["2: MET > 10 GeV"] += 1;
404  } else {
405  passed = false;
406  }
407 
408 
409  // If event contains at least 1 jet
410 
411 
412  _jet_n = 0;
413  _jet_pt.clear();
414  _jet_eta.clear();
415  _jet_flavor.clear();
416  _jet_nconstituents.clear();
417  _jet_sip3dtag.clear();
418  _jet_ktag.clear();
419  _jet_k1_pt.clear();
420  _jet_k1_pid.clear();
421  _jet_k1_sIP3D.clear();
422  _jet_k2_pt.clear();
423  _jet_k2_pid.clear();
424  _jet_k2_sIP3D.clear();
425 
426  _jet_e1_pt.clear();
427  _jet_e1_sIP3D.clear();
428  _jet_e2_pt.clear();
429  _jet_e2_sIP3D.clear();
430 
431  _jet_mu1_pt.clear();
432  _jet_mu1_sIP3D.clear();
433  _jet_mu2_pt.clear();
434  _jet_mu2_sIP3D.clear();
435 
436  _jet_t1_sIP3D.clear();
437  _jet_t2_sIP3D.clear();
438  _jet_t3_sIP3D.clear();
439  _jet_t4_sIP3D.clear();
440 
441  _jet_t1_pt.clear();
442  _jet_t2_pt.clear();
443  _jet_t3_pt.clear();
444  _jet_t4_pt.clear();
445 
446 
447  _jet_mlp_ktagger.clear();
448  _jet_mlp_eltagger.clear();
449  _jet_mlp_mutagger.clear();
450  _jet_mlp_ip3dtagger.clear();
451  _jet_mlp_globaltagger.clear();
452 
453 
454  _jet_etag.clear();
455  _jet_mutag.clear();
456  _jet_charge.clear();
457 
458  _jet_Ks_mass.clear();
459  _jet_Ks_p.clear();
460  _jet_Ks_flightlength.clear();
461  _jet_Ks_sumpt.clear();
462  _jet_K_sumpt.clear();
463  _jet_Ks_zhadron.clear();
464  _jet_K_zhadron.clear();
465  _jet_Ks_leading_zhadron.clear();
466  _jet_K_leading_zhadron.clear();
467  _jet_ehadoveremratio.clear();
468 
469  auto tracks = getTracks();
470 
471 
472  bool use_kaons = true; // Kaons, for now, are constructed from PID maps in this code.
473 
474  bool use_electrons = false;
475 
476  if (DataStore->find("Electrons") != DataStore->end()) {
477  // store the number of electrons in the jets
478  use_electrons = true;
479  }
480 
481  bool use_muons = false;
482 
483  if (DataStore->find("Muons") != DataStore->end()) {
484  // store the number of muons in the jets
485  use_muons = true;
486  }
487 
488 
489  // Loop over tracks in the event; make opposite-sign pairs; compute mass and
490  // call them Ks if within some window
491  // of the Ks0 mass.
492 
493  // auto particles = getGenParticles();
494 
495  std::vector<Track *> all_tracks;
496 
497  // std::vector<GenParticle*> all_tracks;
498  // std::vector<TLorentzVector> all_tracks;
499  for (auto obj_track : *tracks)
500  {
501  auto track = static_cast<Track *>(obj_track);
502 
503  if (track->PT < 0.1) continue;
504  all_tracks.push_back(track);
505 
506 
507  // if (TMath::Abs(track->PID) == 321 && track->P4().Vect().Mag() > 3.5 && track->P4().Vect().Mag() < 4.50 && track->Eta < -1.0 ) {
508 
509  // std::cout << " True Kaon" << std::endl;
510 
511  // if (_branch_mRICHTracks != nullptr) {
512  // for (int itrk = 0; itrk < _branch_mRICHTracks->GetEntries(); itrk++) {
513  // Track *pid_track = (Track *)_branch_mRICHTracks->At(itrk);
514  // if (pid_track->Particle == track->Particle) {
515  // std::cout << " ID'd as " << pid_track->PID << std::endl;
516  // }
517  // }
518  // }
519  // }
520 
521 
522 
523  }
524 
525 
526  // Build all Ks candidates
527  std::vector<Candidate> Ks_candidates;
528 
529  for (int i = 0; i < all_tracks.size(); i++)
530  {
531  for (int j = i + 1; j < all_tracks.size(); j++)
532  {
533  auto track1 = all_tracks[i];
534  auto track2 = all_tracks[j];
535 
536  if (track1->Charge * track2->Charge != -1) continue;
537 
538  // Treat the tracks under the charged pion hypothesis
539  auto track1P4 = TLorentzVector();
540  track1P4.SetPtEtaPhiM(track1->PT, track1->Eta, track1->Phi, _mpi);
541  auto track2P4 = TLorentzVector();
542  track2P4.SetPtEtaPhiM(track2->PT, track2->Eta, track2->Phi, _mpi);
543 
544  // std::cout <<
545  // "===========================================================" <<
546  // std::endl;
547  // track1->P4().Print();
548  // track2->P4().Print();
549 
550  // TLorentzVector Ks_candidate = track1->P4() + track2->P4();
551  TLorentzVector Ks_candidate = track1P4 + track2P4;
552 
553  if (TMath::Abs(Ks_candidate.M() - 0.497) > 0.250) continue;
554 
555  // Spatial coincidence requirement
556  TVector3 track1_POCA(track1->X, track1->Y, track1->Z);
557  TVector3 track2_POCA(track2->X, track2->Y, track2->Z);
558 
559  auto intertrack_displacement = track1_POCA - track2_POCA;
560  float intertrack_distance = intertrack_displacement.Mag();
561 
562  float track1_d0err = track1->ErrorD0;
563  float track1_z0err = track1->ErrorDZ;
564  float track2_d0err = track2->ErrorD0;
565  float track2_z0err = track2->ErrorDZ;
566 
567  float err_3D = TMath::Sqrt(TMath::Power(track1_d0err, 2.0) + TMath::Power(track1_z0err, 2.0) +
568  TMath::Power(track2_d0err, 2.0) + TMath::Power(track2_z0err, 2.0));
569 
570  float intertrack_distance_signif = intertrack_distance / err_3D;
571 
572  // std::cout << " Intertrack distance, error, significance: " <<
573  // intertrack_distance << " mm"
574  // << ", " << err_3D
575  // << ", " << intertrack_distance_signif << std::endl;
576 
577  // To be from a common decay, their displacement significance should be
578  // small
579 
580  if (intertrack_distance_signif > 1.5) continue;
581 
582  // build a new Candidate
583  Candidate Ks;
584  Ks.PID = 310;
585  Ks.Mass = Ks_candidate.M();
586  Ks.Momentum = Ks_candidate;
587 
588  // midpoint between the pocas of the two tracks
589  Ks.Position = TLorentzVector(track2_POCA + 0.5 * intertrack_displacement, 0.0);
590 
591  Ks_candidates.push_back(Ks);
592  }
593  }
594 
595 
596  // Get Kaons from the EIC Detector PID Systems
597  std::vector<Track *> kaon_candidates;
598  std::vector<Track> pid_candidates;
599 
600  // The mRICH and barrel DIRC are physicall distinct systems covering different
601  // eta ranges
602  // The dualRICH system has two components - an aerogel detector for low-momentum
603  // performance
604  // and a CF6 system for high-momentum performance.
605 
606 
607  if (_branch_mRICHTracks != nullptr) {
608  for (int itrk = 0; itrk < _branch_mRICHTracks->GetEntries(); itrk++) {
609  Track *track = (Track *)_branch_mRICHTracks->At(itrk);
610 
611  if ((track->Eta < -3.5) || (-1.0 < track->Eta)) continue;
612 
613  pid_candidates.push_back(*track);
614 
615  Int_t reco_pid = track->PID;
616 
617  if (TMath::Abs(reco_pid) == 321) {
618  kaon_candidates.push_back(track);
619  }
620  }
621  }
622 
623  if (_branch_barrelDircTracks != nullptr) {
624  for (int itrk = 0; itrk < _branch_barrelDircTracks->GetEntries(); itrk++) {
625  Track *track = (Track *)_branch_barrelDircTracks->At(itrk);
626 
627  if ((track->Eta < -1.0) || (1.0 < track->Eta)) continue;
628 
629  pid_candidates.push_back(*track);
630 
631  Int_t reco_pid = track->PID;
632 
633  if (TMath::Abs(reco_pid) == 321) {
634  kaon_candidates.push_back(track);
635  }
636  }
637  }
638 
639 
640  // Handle tracks in the forward direction (dualRICH)
641  for (auto track : all_tracks) {
642  if (1.0 <= track->Eta && track->Eta <= 3.5) {
643  Int_t final_pid = 0;
644  Double_t p_track = track->P4().Vect().Mag();
645  Double_t ag_p_threshold = 12.0;
646 
647  if (p_track < ag_p_threshold) {
648  // region of sensitivity for Aerogel
649  if (_branch_dualRICHagTracks != nullptr) {
650  for (int itrk = 0; itrk < _branch_dualRICHagTracks->GetEntries(); itrk++) {
651  Track *track_ag = (Track *)_branch_dualRICHagTracks->At(itrk);
652  if (track_ag->Particle == track->Particle) {
653  final_pid = track_ag->PID;
654  break;
655  }
656  }
657  }
658  } else {
659  if (_branch_dualRICHcfTracks != nullptr) {
660  for (int itrk = 0; itrk < _branch_dualRICHcfTracks->GetEntries(); itrk++) {
661  Track *track_cf = (Track *)_branch_dualRICHcfTracks->At(itrk);
662  if (track_cf->Particle == track->Particle) {
663  final_pid = track_cf->PID;
664  break;
665  }
666  }
667  }
668  }
669 
670  Track drich_track = *track;
671  drich_track.PID = final_pid;
672  pid_candidates.push_back(drich_track);
673 
674  if (TMath::Abs(final_pid) == 321)
675  kaon_candidates.push_back(track);
676 
677  }
678  }
679 
680  _pid_track_n = pid_candidates.size();
681 
682  for (auto track : all_tracks) {
683  // try to identify the jet that might contain this track and save its
684  // mother's identity (e.g. 1, 4, 5, etc.)
685  float jet_min_dR = 1e9;
686  int jet_mother_id = 0;
687  float jet_mother_pt = 0;
688  float jet_mother_eta = 1e6;
689 
690  for (int ijet = 0; ijet < getJets()->GetEntries(); ijet++)
691  {
692  // Take first jet
693  Jet *jet = (Jet *)getJets()->At(ijet);
694  float this_dR = jet->P4().DeltaR(track->P4());
695 
696  if ((this_dR < 1.0) && (this_dR < jet_min_dR)) {
697  jet_min_dR = this_dR;
698  jet_mother_id = jet->Flavor;
699  jet_mother_pt = jet->PT;
700  jet_mother_eta = jet->Eta;
701  }
702  }
703  _pid_track_jetmother.push_back(jet_mother_id);
704  _pid_track_jet_pt.push_back(jet_mother_pt);
705  _pid_track_jet_eta.push_back(jet_mother_eta);
706 
707  Int_t pid_id = 0;
708  for (auto pid_track : pid_candidates) {
709  if ( pid_track.Particle == track->Particle ) {
710  pid_id = pid_track.PID;
711  break;
712  }
713  }
714 
715  // if (TMath::Abs(track->PID) == 321 && track->P4().Vect().Mag() > 3.5 && track->P4().Vect().Mag() < 4.50 && track->Eta < -1.0 ) {
716 
717  // std::cout << " > True Kaon" << std::endl;
718 
719  // if (_branch_mRICHTracks != nullptr) {
720  // for (int itrk = 0; itrk < _branch_mRICHTracks->GetEntries(); itrk++) {
721  // Track *pid_track = (Track *)_branch_mRICHTracks->At(itrk);
722  // if (pid_track->Particle == track->Particle) {
723  // std::cout << " ID'd as " << pid_track->PID << std::endl;
724  // }
725  // }
726  // }
727  // std::cout << " Alt. ID'd as " << pid_id << std::endl;
728  // }
729 
730  _pid_track_true_pid.push_back(track->PID);
731  _pid_track_reco_pid.push_back( pid_id );
732  _pid_track_pt.push_back(track->PT);
733  _pid_track_eta.push_back(track->Eta);
734  }
735 
736 
737  std::vector<Jet *> all_jets;
738 
739  for (int ijet = 0; ijet < getJets()->GetEntries(); ijet++)
740  {
741  // Take first jet
742  Jet *jet = (Jet *)getJets()->At(ijet);
743  all_jets.push_back(jet);
744 
745  _jet_pt.push_back(jet->PT);
746  _jet_eta.push_back(jet->Eta);
747  _jet_flavor.push_back(jet->Flavor);
748  _jet_nconstituents.push_back(jet->Constituents.GetEntries());
749 
750  // pre-long-lived particle decay optimization
751  // _jet_sip3dtag.push_back( Tagged_sIP3D(jet, *getEFlowTracks(), 3.75, 1.00,
752  // 2.0) );
753  // after including long-lived particle decay, optimization
754  // _jet_sip3dtag.push_back( Tagged_sIP3D(jet, *getEFlowTracks(), 2.75, 0.25,
755  // 2.0) );
756  // _jet_sip3dtag.push_back( Tagged_sIP3D(jet, *tagging_tracks, 2.75, 0.25,
757  // 2.0) );
758  // _jet_sip3dtag.push_back( Tagged_sIP3D(jet, *tagging_tracks, 3.75, 1.00,
759  // 2.0) );
760 
761  // After using the significance of the charm tag yield in 100/fb to optimize
762  // the tagging,
763  // and after fixing the PYTHIA ctau_max bug:
764  _jet_sip3dtag.push_back(Tagged_sIP3D(jet, *getEFlowTracks(), 3.00, 0.25, 2.0));
765 
766 
767  // Find all the good flavor-tagging tracks in this jet
768  std::vector<Track *> jet_tracks;
769 
770  for (int itrk = 0; itrk < getEFlowTracks()->GetEntries(); itrk++) {
771  auto track = (Track *)getEFlowTracks()->At(itrk);
772 
773  if ((track->P4().DeltaR(jet->P4()) < 0.5) && IsTaggingTrack(track)) {
774  jet_tracks.push_back(track);
775  }
776  }
777 
778 
779  // Rank by sIP3D
780  std::sort(jet_tracks.begin(), jet_tracks.end(), [jet](auto& lhs, const auto& rhs)
781  {
782  return sIP3D(jet, lhs) > sIP3D(jet, rhs);
783  });
784 
785  // store the HIP track information
786  if (jet_tracks.size() > 0) {
787  _jet_t1_sIP3D.push_back(sIP3D(jet, jet_tracks[0]));
788  _jet_t1_pt.push_back(jet_tracks[0]->PT);
789  } else {
790  _jet_t1_sIP3D.push_back(-199.0);
791  _jet_t1_pt.push_back(-1.0);
792  }
793 
794  if (jet_tracks.size() > 1) {
795  _jet_t2_sIP3D.push_back(sIP3D(jet, jet_tracks[1]));
796  _jet_t2_pt.push_back(jet_tracks[1]->PT);
797  } else {
798  _jet_t2_sIP3D.push_back(-199.0);
799  _jet_t2_pt.push_back(-1.0);
800  }
801 
802  if (jet_tracks.size() > 2) {
803  _jet_t3_sIP3D.push_back(sIP3D(jet, jet_tracks[2]));
804  _jet_t3_pt.push_back(jet_tracks[2]->PT);
805  } else {
806  _jet_t3_sIP3D.push_back(-199.0);
807  _jet_t3_pt.push_back(-1.0);
808  }
809 
810  if (jet_tracks.size() > 3) {
811  _jet_t4_sIP3D.push_back(sIP3D(jet, jet_tracks[3]));
812  _jet_t4_pt.push_back(jet_tracks[3]->PT);
813  } else {
814  _jet_t4_sIP3D.push_back(-199.0);
815  _jet_t4_pt.push_back(-1.0);
816  }
817 
818 
819  if (use_electrons) {
820  _jet_etag.push_back(Tagged_Electron(jet, std::any_cast<std::vector<Track *> >((*DataStore)["Electrons"]), 3.0, 1.0, 1));
821  } else {
822  _jet_etag.push_back(0.0);
823  }
824 
825 
826  if (use_muons) {
827  _jet_mutag.push_back(Tagged_Muon(jet, std::any_cast<std::vector<Track *> >((*DataStore)["Muons"]), 3.0, 1.0, 1));
828  } else {
829  _jet_mutag.push_back(0.0);
830  }
831 
832  if (use_kaons) {
833  _jet_ktag.push_back(Tagged_Kaon(jet, std::any_cast<std::vector<Track *> >(kaon_candidates), 3.0, 1.0, 1));
834  } else {
835  _jet_ktag.push_back(0.0);
836  }
837 
838  // PID-system kaons
839 
840  // Sort by PT
841  std::sort(kaon_candidates.begin(), kaon_candidates.end(), [](auto& lhs, const auto& rhs)
842  {
843  return lhs->PT > rhs->PT;
844  });
845 
846  auto kaons_list = kaon_candidates;
847  auto electrons_list = std::any_cast<std::vector<Track *> >((*DataStore)["Electrons"]);
848  auto muons_list = std::any_cast<std::vector<Track *> >((*DataStore)["Muons"]);
849 
850  // if (kaon_candidates.size() > 0) {
851  if (kaons_list.size() > 0) {
852  auto k1 = kaons_list.at(0);
853  _jet_k1_pt.push_back(k1->PT);
854  _jet_k1_sIP3D.push_back(sIP3D(jet, k1));
855  _jet_k1_pid.push_back(k1->PID);
856  } else {
857  _jet_k1_pt.push_back(-1.0);
858  _jet_k1_sIP3D.push_back(-199.0);
859  _jet_k1_pid.push_back(0);
860  }
861 
862  // if (kaon_candidates.size() > 1) {
863  if (kaons_list.size() > 1) {
864  auto k2 = kaons_list.at(1);
865  _jet_k2_pt.push_back(k2->PT);
866  _jet_k2_sIP3D.push_back(sIP3D(jet, k2));
867  _jet_k2_pid.push_back(k2->PID);
868  } else {
869  _jet_k2_pt.push_back(-1.0);
870  _jet_k2_sIP3D.push_back(-199.0);
871  _jet_k2_pid.push_back(0);
872  }
873 
874  TVector3 Ks_sumpt;
875  _jet_Ks_mass.push_back(std::vector<float>());
876  _jet_Ks_p.push_back(std::vector<float>());
877  _jet_Ks_flightlength.push_back(std::vector<float>());
878  _jet_Ks_zhadron.push_back(std::vector<float>());
879 
880  for (auto Ks_candidate : Ks_candidates) {
881  if (Ks_candidate.Position.Rho() < 5) // 5mm minimum displacement from IP
882  continue;
883 
884  if (Ks_candidate.Momentum.DeltaR(jet->P4()) < 0.5) {
885  Ks_sumpt += Ks_candidate.Momentum.Vect();
886  _jet_Ks_mass[ijet].push_back(Ks_candidate.Mass);
887  _jet_Ks_p[ijet].push_back(Ks_candidate.Momentum.Rho());
888  _jet_Ks_flightlength[ijet].push_back(Ks_candidate.Position.Rho());
889 
890  _jet_Ks_zhadron[ijet].push_back(TMath::Abs((Ks_candidate.Momentum.Vect() * jet->P4().Vect()) / jet->P4().Vect().Mag2()));
891  }
892  }
893 
894  if (_jet_Ks_zhadron[ijet].size() == 0) {
895  _jet_Ks_zhadron[ijet].push_back(-1.0); // TMVA hates empty vector
896  // branches!
897  }
898 
899  _jet_Ks_leading_zhadron.push_back(*max_element(_jet_Ks_zhadron[ijet].begin(), _jet_Ks_zhadron[ijet].end()));
900 
901 
902  _jet_Ks_sumpt.push_back(Ks_sumpt.Perp());
903 
904  // handle charged kaons
905  _jet_K_zhadron.push_back(std::vector<float>());
906  TVector3 K_sumpt;
907 
908  if (use_kaons) {
909  auto kaons = kaon_candidates;
910 
911  for (auto kaon : kaons) {
912  if (kaon->P4().DeltaR(jet->P4()) < 0.5) {
913  K_sumpt += kaon->P4().Vect();
914  _jet_K_zhadron[ijet].push_back(TMath::Abs((kaon->P4().Vect() * jet->P4().Vect()) / jet->P4().Vect().Mag2()));
915  }
916  }
917  }
918 
919  if (_jet_K_zhadron[ijet].size() == 0) {
920  _jet_K_zhadron[ijet].push_back(-1.0); // TMVA hates empty vector branches!
921  }
922 
923  _jet_K_leading_zhadron.push_back(*max_element(_jet_K_zhadron[ijet].begin(), _jet_K_zhadron[ijet].end()));
924 
925  _jet_K_sumpt.push_back(K_sumpt.Perp());
926 
927  // Electrons
928  if (electrons_list.size() > 0) {
929  auto e1 = electrons_list.at(0);
930  _jet_e1_pt.push_back(e1->PT);
931  _jet_e1_sIP3D.push_back(sIP3D(jet, e1));
932  } else {
933  _jet_e1_pt.push_back(-1.0);
934  _jet_e1_sIP3D.push_back(-199.0);
935  }
936 
937  if (electrons_list.size() > 1) {
938  auto e2 = electrons_list.at(1);
939  _jet_e2_pt.push_back(e2->PT);
940  _jet_e2_sIP3D.push_back(sIP3D(jet, e2));
941  } else {
942  _jet_e2_pt.push_back(-1.0);
943  _jet_e2_sIP3D.push_back(-199.0);
944  }
945 
946  // Muons
947  if (muons_list.size() > 0) {
948  auto mu1 = muons_list.at(0);
949  _jet_mu1_pt.push_back(mu1->PT);
950  _jet_mu1_sIP3D.push_back(sIP3D(jet, mu1));
951  } else {
952  _jet_mu1_pt.push_back(-1.0);
953  _jet_mu1_sIP3D.push_back(-199.0);
954  }
955 
956  if (muons_list.size() > 1) {
957  auto mu2 = muons_list.at(1);
958  _jet_mu2_pt.push_back(mu2->PT);
959  _jet_mu2_sIP3D.push_back(sIP3D(jet, mu2));
960  } else {
961  _jet_mu2_pt.push_back(-1.0);
962  _jet_mu2_sIP3D.push_back(-199.0);
963  }
964 
965  // MLP Tagger Decisions
966  _mva_inputs_int["jet_flavor"] = static_cast<Int_t>(jet->Flavor);
967  _mva_inputs_float["jet_pt"] = jet->PT;
968  _mva_inputs_float["jet_eta"] = jet->Eta;
969  _mva_inputs_float["met_et"] = _met_et;
970 
971 
972  _mva_inputs_float["jet_t1_pt"] = _jet_t1_pt.back();
973  _mva_inputs_float["jet_t1_sIP3D"] = _jet_t1_sIP3D.back();
974  _mva_inputs_float["jet_t2_pt"] = _jet_t2_pt.back();
975  _mva_inputs_float["jet_t2_sIP3D"] = _jet_t2_sIP3D.back();
976  _mva_inputs_float["jet_t3_pt"] = _jet_t3_pt.back();
977  _mva_inputs_float["jet_t3_sIP3D"] = _jet_t3_sIP3D.back();
978  _mva_inputs_float["jet_t4_pt"] = _jet_t4_pt.back();
979  _mva_inputs_float["jet_t4_sIP3D"] = _jet_t4_sIP3D.back();
980 
981  _mva_inputs_float["jet_k1_pt"] = _jet_k1_pt.back();
982  _mva_inputs_float["jet_k2_pt"] = _jet_k2_pt.back();
983  _mva_inputs_float["jet_k1_sIP3D"] = _jet_k1_sIP3D.back();
984  _mva_inputs_float["jet_k2_sIP3D"] = _jet_k2_sIP3D.back();
985 
986  _mva_inputs_float["jet_e1_pt"] = _jet_e1_pt.back();
987  _mva_inputs_float["jet_e2_pt"] = _jet_e2_pt.back();
988  _mva_inputs_float["jet_e1_sIP3D"] = _jet_e1_sIP3D.back();
989  _mva_inputs_float["jet_e2_sIP3D"] = _jet_e2_sIP3D.back();
990 
991  _mva_inputs_float["jet_mu1_pt"] = _jet_mu1_pt.back();
992  _mva_inputs_float["jet_mu2_pt"] = _jet_mu2_pt.back();
993  _mva_inputs_float["jet_mu1_sIP3D"] = _jet_mu1_sIP3D.back();
994  _mva_inputs_float["jet_mu2_sIP3D"] = _jet_mu2_sIP3D.back();
995 
996 
997  _jet_mlp_ip3dtagger.push_back(_mva_reader_ip3dtagger->EvaluateMVA("CharmIP3DTagger"));
998  _jet_mlp_ktagger.push_back(_mva_reader_ktagger->EvaluateMVA("CharmKTagger"));
999  _jet_mlp_eltagger.push_back(_mva_reader_eltagger->EvaluateMVA("CharmETagger"));
1000  _jet_mlp_mutagger.push_back(_mva_reader_mutagger->EvaluateMVA("CharmMuTagger"));
1001 
1002  _mva_inputs_float["jet_mlp_ip3dtagger"] = _mva_reader_ip3dtagger->EvaluateMVA("CharmIP3DTagger");
1003  _mva_inputs_float["jet_mlp_ktagger"] = _mva_reader_ktagger->EvaluateMVA("CharmKTagger");
1004  _mva_inputs_float["jet_mlp_eltagger"] = _mva_reader_eltagger->EvaluateMVA("CharmETagger");
1005  _mva_inputs_float["jet_mlp_mutagger"] = _mva_reader_mutagger->EvaluateMVA("CharmMuTagger");
1006 
1007  _jet_mlp_globaltagger.push_back(_mva_reader_globaltagger->EvaluateMVA("CharmGlobalTagger"));
1008 
1009  // Calorimeter energy ratios
1010  _jet_ehadoveremratio.push_back(jet->EhadOverEem);
1011 
1012  // Jet charge!
1013  float jet_Q = -10.0;
1014  float kappa = 0.5;
1015 
1016  for (auto obj_track : *getEFlowTracks()) {
1017  auto track = static_cast<Track *>(obj_track);
1018 
1019  if (track->P4().DeltaR(jet->P4()) < 0.5) {
1020  if (jet_Q == -10.0) {
1021  jet_Q = 0.0;
1022  }
1023  jet_Q += track->Charge * TMath::Power(track->PT, kappa);
1024  }
1025  }
1026  jet_Q /= TMath::Power(jet->PT, kappa);
1027  _jet_charge.push_back(jet_Q);
1028  }
1029 
1030  _jet_n = _jet_pt.size();
1031 
1032  std::vector<Jet *> fiducial_jets = SelectorFcn<Jet>(all_jets, [](Jet *j) {
1033  return TMath::Abs(j->Eta) < 3.0 && j->PT > 5.0;
1034  });
1035 
1036  if ((passed == true) && (fiducial_jets.size() > 0)) {
1037  _cut_flow["3: Fiducial Jets >= 1"] += 1;
1038  } else {
1039  passed = false;
1040  }
1041 
1042  // std::vector<Jet*> charmJets = SelectorFcn<Jet>(fiducial_jets, [](Jet* j){
1043  // return (j->Flavor == 4); });
1044 
1045  std::vector<Jet *> charmJets;
1046 
1047  if (DataStore->find("CharmJets") != DataStore->end()) {
1048  charmJets = std::any_cast<std::vector<Jet *> >((*DataStore)["CharmJets"]);
1049  }
1050 
1051  if ((passed == true) && (charmJets.size() > 0)) {
1052  _cut_flow["4: Charm Jet == 1"] += 1;
1053  } else {
1054  passed = false;
1055  }
1056 
1057 
1058  // Store charm jet information
1059  if (passed) {
1060  for (auto jet : charmJets) {
1061  _charmjet_pt.push_back(jet->PT);
1062  _charmjet_eta.push_back(jet->Eta);
1063  }
1064  _charmjet_n = _charmjet_pt.size();
1065  }
1066 
1067 
1068  return true;
1069 }