TTTT Analysis  0.1
TTTT_Analysis.cpp
Go to the documentation of this file.
1 
37 #include <iostream>
38 #include <vector>
39 #include <utility>
40 #include <limits>
41 
42 #include "filval/filval.hpp"
43 #include "filval_root/filval_root.hpp"
44 
45 #include "MiniTreeDataSet.hpp"
46 #include <TSystem.h>
47 
48 #define PI 3.14159
49 #define W_MASS 80.385 // GeV/c^2
50 #define Z_MASS 91.188 // GeV/c^2
51 #define T_MASS 172.44 // GeV/c^2
52 
53 using namespace std;
54 using namespace fv;
55 using namespace fv::root;
56 
57 void enable_branches(MiniTreeDataSet& mt){
58  mt.fChain->SetBranchStatus("*", false);
59 
60  mt.track_branch<int>("nLepGood");
61  mt.track_branch_vec< int >("nLepGood", "LepGood_pdgId");
62  mt.track_branch_vec<float>("nLepGood", "LepGood_pt");
63  mt.track_branch_vec<float>("nLepGood", "LepGood_eta");
64  mt.track_branch_vec<float>("nLepGood", "LepGood_phi");
65  mt.track_branch_vec<float>("nLepGood", "LepGood_mass");
66  mt.track_branch_vec< int >("nLepGood", "LepGood_charge");
67  mt.track_branch_vec< int >("nLepGood", "LepGood_mcMatchId");
68  mt.track_branch_vec< int >("nLepGood", "LepGood_mcMatchPdgId");
69  mt.track_branch_vec< int >("nLepGood", "LepGood_mcMatchAny");
70  mt.track_branch_vec< int >("nLepGood", "LepGood_mcMatchTau");
71  mt.track_branch_vec< int >("nLepGood", "LepGood_mcPt");
72 
73  mt.track_branch<int>("nJet");
74  mt.track_branch_vec<float>("nJet", "Jet_pt");
75  mt.track_branch_vec<float>("nJet", "Jet_eta");
76  mt.track_branch_vec<float>("nJet", "Jet_phi");
77  mt.track_branch_vec<float>("nJet", "Jet_mass");
78  mt.track_branch_vec<float>("nJet", "Jet_btagCMVA");
79  mt.track_branch_vec< int >("nJet", "Jet_mcMatchFlav");
80  mt.track_branch_vec< int >("nJet", "Jet_mcMatchId");
81  mt.track_branch_vec< int >("nJet", "Jet_mcFlavour");
82 
83 
84  mt.track_branch<int>("nGenPart");
85  mt.track_branch_vec< int >("nGenPart", "GenPart_pdgId");
86  mt.track_branch_vec< int >("nGenPart", "GenPart_motherIndex");
87  mt.track_branch_vec< int >("nGenPart", "GenPart_motherId");
88  mt.track_branch_vec<float>("nGenPart", "GenPart_pt");
89  mt.track_branch_vec<float>("nGenPart", "GenPart_eta");
90  mt.track_branch_vec<float>("nGenPart", "GenPart_phi");
91  mt.track_branch_vec<float>("nGenPart", "GenPart_mass");
92  mt.track_branch_vec< int >("nGenPart", "GenPart_status");
93 
94  mt.track_branch<int>("nBJetLoose40");
95  mt.track_branch<int>("nBJetMedium40");
96  mt.track_branch<int>("nBJetTight40");
97 
98 
99  mt.track_branch<int>("nVert");
100 
101  mt.track_branch< int >("run" );
102  mt.track_branch< int >("lumi");
103  mt.track_branch< int >("evt" );
104  mt.track_branch<float>("xsec");
105 }
106 
107 struct Jet{
108  TLorentzVector v;
109  int idx;
110  int pdgid;
111  float b_cmva;
112  Jet() { }
113  Jet(const TLorentzVector& v, int idx, int pdgid, float b_cmva)
114  :v(v),idx(idx),pdgid(pdgid),b_cmva(b_cmva) { }
115 
116  static Jet reco(const TLorentzVector& v, int idx, float b_cmva){
117  return Jet(v, idx, 0, b_cmva);
118  }
119 
120  static Jet mc(const TLorentzVector& v, int idx, int pdgid){
121  return Jet(v, idx, pdgid, 0);
122  }
123 };
124 
125 void declare_values(MiniTreeDataSet& mt){
126 
127  // Define Lorentz Vector(TLorentzVector) object from fields in ntuple
128  lorentz_vectors("LepGood_pt", "LepGood_eta", "LepGood_phi", "LepGood_mass", "LepGood_4v");
129  lorentz_vectors("GenPart_pt", "GenPart_eta", "GenPart_phi", "GenPart_mass", "GenPart_4v");
130  lorentz_vectors("Jet_pt", "Jet_eta", "Jet_phi", "Jet_mass", "Jet_4v" );
131 
132  energies("GenPart_4v", "GenPart_energy");
133 
134 
135  // Define a couple selections to be used in the top-mass reconstruction.
136  auto& b_mva_filter = GenFunction::register_function<bool(Jet)>("b_mva_filter",
137  FUNC(([cut=0.0](const Jet& j){
138  return j.b_cmva > cut;
139  })));
140  auto& b_pdgid_filter = GenFunction::register_function<bool(Jet)>("b_pdgid_filter",
141  FUNC(([](const Jet& j){
142  return j.pdgid == 5 || j.pdgid==-5;
143  })));
144  auto& w_mass_filter = GenFunction::register_function<bool(Jet, Jet)>("w_mass_filter",
145  FUNC(([win_l=W_MASS-10, win_h=W_MASS+10](const Jet& j1, const Jet& j2){
146  float inv_mass = (j1.v + j2.v).M();
147  return inv_mass > win_l && inv_mass < win_h;
148  })));
149  auto& dup_filter = GenFunction::register_function<bool(std::tuple<Jet,Jet>,Jet)>("dup_filter",
150  FUNC(([](const std::tuple<Jet,Jet>& w, const Jet& b){
151  int j0 = b.idx;
152  int j1 = std::get<0>(w).idx;
153  int j2 = std::get<1>(w).idx;
154  return (j0 != j1) && (j0 != j2) && (j1 != j2);
155  })));
156  auto& qg_id_filter = GenFunction::register_function<bool(Jet, Jet)>("qg_id_filter",
157  FUNC(([](const Jet& j1, const Jet& j2){
158  // require both particles be either quarks(not Top) or gluons
159  int id1 = abs(j1.pdgid);
160  int id2 = abs(j2.pdgid);
161  return ((id1 >=1 && id1 <= 5) || id1 == 21) &&
162  ((id2 >=1 && id2 <= 5) || id2 == 21);
163  })));
164 
165  // Here is the calculation of the Top Reconstructed Mass from Jets
166  auto jets = apply(GenFunction::register_function<std::vector<Jet>(std::vector<TLorentzVector>,std::vector<float>)>("build_reco_jets",
167  FUNC(([](const std::vector<TLorentzVector>& vs, const std::vector<float>& b_cmvas){
168  std::vector<Jet> jets;
169  for(int i=0; i<vs.size(); i++){
170  jets.push_back(Jet::reco(vs[i],i, b_cmvas[i]));
171  }
172  return jets;
173  }))), fv::tuple(lookup<std::vector<TLorentzVector>>("Jet_4v"), lookup<std::vector<float>>("Jet_btagCMVA")), "reco_jets");
174 
175 
176  auto b_jets = filter(b_mva_filter, jets, "reco_b_jets");
177  auto w_dijets = tup_filter<Jet,Jet>(w_mass_filter, combinations<Jet,2>(jets, "reco_dijets"));
178 
179  auto top_cands = cart_product<std::tuple<Jet,Jet>, Jet>(w_dijets, b_jets);
180 
181  top_cands = tup_filter(dup_filter, top_cands);
182 
183  auto& t_mass = GenFunction::register_function<float(std::tuple<Jet,Jet>,Jet)>("t_mass",
184  FUNC(([](const std::tuple<Jet,Jet>& w, const Jet& b){
185  return (std::get<0>(w).v+std::get<1>(w).v+b.v).M();
186  })));
187 
188  fv::map(t_mass, top_cands, "reco_top_mass");
189 
190  // Here is the calculation of the Top Reconstructed Mass from Generator-Level objects
191  jets = apply(GenFunction::register_function<std::vector<Jet>(std::vector<TLorentzVector>,std::vector<int>)>("build_mcjets",
192  FUNC(([](const std::vector<TLorentzVector>& vs, const std::vector<int>& pdgid){
193  std::vector<Jet> jets;
194  for(int i=0; i<vs.size(); i++){
195  jets.push_back(Jet::mc(vs[i],i, pdgid[i]));
196  }
197  return jets;
198  }))), fv::tuple(lookup<std::vector<TLorentzVector>>("GenPart_4v"), lookup<std::vector<int>>("GenPart_pdgId")), "mcjets");
199 
200  b_jets = filter(b_pdgid_filter, jets);
201 
202  w_dijets = tup_filter(qg_id_filter, combinations<Jet,2>(jets));
203  w_dijets = tup_filter(w_mass_filter, w_dijets);
204 
205  top_cands = cart_product<std::tuple<Jet,Jet>, Jet>(w_dijets, b_jets);
206 
207  top_cands = tup_filter(dup_filter, top_cands);
208 
209 
210  fv::map(t_mass, top_cands, "mc_top_mass");
211 
212 
213  // calculation of di-jet inv-mass spectrum
214  auto& inv_mass2 = GenFunction::register_function<float(Jet, Jet)>("inv_mass2",
215  FUNC(([] (const Jet& j1, const Jet& j2){
216  TLorentzVector sum = j1.v + j2.v;
217  return (float)sum.M();
218  })));
219  fv::map(inv_mass2, lookup<std::vector<std::tuple<Jet,Jet>>>("reco_dijets"), "dijet_inv_mass");
220 
221 
222 
223  count<float>(GenFunction::register_function<bool(float)>("bJet_Selection",
224  FUNC(([](float x){
225  return x>0;
226  }))), "Jet_btagCMVA", "b_jet_count");
227 
228  auto &is_electron = GenFunction::register_function<bool(int)>("is_electron",
229  FUNC(([](int pdgId) {
230  return abs(pdgId) == 11;
231  })));
232  auto &is_muon = GenFunction::register_function<bool(int)>("is_muon",
233  FUNC(([](int pdgId) {
234  return abs(pdgId) == 13;
235  })));
236  auto &is_lepton = GenFunction::register_function<bool(int)>("is_lepton",
237  FUNC(([ie=&is_electron, im=&is_muon](int pdgId) {
238  return (*ie)(pdgId) || (*im)(pdgId);
239  })));
240 
241  count<int>(is_electron, "GenPart_pdgId", "genEle_count");
242  count<int>(is_muon, "GenPart_pdgId", "genMu_count");
243  count<int>(is_lepton, "GenPart_pdgId", "genLep_count");
244 
245  count<int>(is_electron, "LepGood_pdgId", "recEle_count");
246  count<int>(is_muon, "LepGood_pdgId", "recMu_count");
247  count<int>(is_lepton, "LepGood_pdgId", "recLep_count");
248 
249 
250  fv::pair<int, int>("genEle_count", "recEle_count", "genEle_count_v_recEle_count");
251  fv::pair<int, int>("genMu_count", "recMu_count", "genMu_count_v_recMu_count");
252  fv::pair<int, int>("genLep_count", "recLep_count", "genLep_count_v_recLep_count");
253 
254  obs_filter("trilepton", FUNC(([nLepGood=lookup<int>("nLepGood")]()
255  {
256  return dynamic_cast<Value<int>*>(nLepGood)->get_value() == 3;
257  })));
258 
259  obs_filter("os-dilepton", FUNC(([LepGood_charge=lookup<vector<int>>("LepGood_charge")]()
260  {
261  auto& charges = LepGood_charge->get_value();
262  return charges.size()==2 && (charges[0] != charges[1]);
263  })));
264 
265  obs_filter("ss-dilepton", FUNC(([LepGood_charge=lookup<vector<int>>("LepGood_charge")]()
266  {
267  auto& charges = LepGood_charge->get_value();
268  return charges.size()==2 && (charges[0] == charges[1]);
269  })));
270 
271 }
272 
273 void declare_containers(MiniTreeDataSet& mt){
274 
275  mt.register_container<ContainerTH1<int>>("lepton_count", "Lepton Multiplicity", lookup<int>("nLepGood"), 8, 0, 8);
276 
277  mt.register_container<ContainerTH1Many<float>>("Jet_pt_dist", "Jet P_T",
278  lookup<vector<float>>("Jet_pt"), 50, 0, 500,
279  "Jet P_T");
280  mt.register_container<ContainerTH1Many<float>>("Jet_eta_dist", "Jet Eta",
281  lookup<vector<float>>("Jet_eta"), 50, -3, 3,
282  "Jet Eta");
283  mt.register_container<ContainerTH1Many<float>>("Jet_phi_dist", "Jet Phi",
284  lookup<vector<float>>("Jet_phi"), 20, -PI, PI,
285  "Jet Phi");
286  mt.register_container<ContainerTH1Many<float>>("Jet_mass_dist", "Jet Mass",
287  lookup<vector<float>>("Jet_mass"), 50, 0, 200,
288  "Jet Mass");
289 
290  mt.register_container<ContainerTH1Many<float>>("dijet_inv_mass", "Di-Jet Inv. Mass - All",
291  lookup<vector<float>>("dijet_inv_mass"), 100, 0, 500,
292  "Di-Jet Mass");
293  mt.register_container<ContainerTH1Many<float>>("dijet_inv_mass_osdilepton", "Di-Jet Inv. Mass - OS Dilepton",
294  lookup<vector<float>>("dijet_inv_mass"), 100, 0, 500,
295  "Di-Jet Mass")->add_filter(lookup_obs_filter("os-dilepton"));
296  mt.register_container<ContainerTH1Many<float>>("dijet_inv_mass_ssdilepton", "Di-Jet Inv. Mass - SS Dilepton",
297  lookup<vector<float>>("dijet_inv_mass"), 100, 0, 500,
298  "Di-Jet Mass")->add_filter(lookup_obs_filter("ss-dilepton"));
299  mt.register_container<ContainerTH1Many<float>>("dijet_inv_mass_trilepton", "Di-Jet Inv. Mass - Trilepton",
300  lookup<vector<float>>("dijet_inv_mass"), 100, 0, 500,
301  "Di-Jet Mass")->add_filter(lookup_obs_filter("trilepton"));
302 
303  mt.register_container<ContainerTH1Many<float>>("reco_top_mass", "Reconstructed Top mass",
304  lookup<vector<float>>("reco_top_mass"), 100, 0, 500,
305  "Tri-jet invarient mass");
306  /* mt.register_container<ContainerTH1Many<float>>("reco_top_pt", "Reconstructed Top Pt", */
307  /* lookup<vector<float>>("reco_top_pt"), 100, 0, 500, */
308  /* "Tri-jet Pt"); */
309 
310  mt.register_container<ContainerTH1Many<float>>("mc_top_mass", "Reconstructed MC Top mass",
311  lookup<vector<float>>("mc_top_mass"), 100, 0, 500,
312  "Tri-jet invarient mass");
313 
314 
315  mt.register_container<ContainerTH2<int>>("nLepvsnJet", "Number of Leptons vs Number of Jets",
316  fv::pair<int, int>("nLepGood", "nJet"),
317  7, 0, 7, 20, 0, 20,
318  "Number of Leptons", "Number of Jets");
319 
320 
321  mt.register_container<ContainerTH2<int>>("genEle_count_v_recEle_count", "Number of Generated Electrons v. Number of Reconstructed Electrons",
322  lookup<std::pair<int,int>>("genEle_count_v_recEle_count"),
323  10, 0, 10, 10, 0, 10,
324  "Generated Electrons","Reconstructed Electrons");
325 
326  mt.register_container<ContainerTH2<int>>("genMu_count_v_recMu_count", "Number of Generated Muons v. Number of Reconstructed Muons",
327  lookup<std::pair<int,int>>("genMu_count_v_recMu_count"),
328  10, 0, 10, 10, 0, 10,
329  "Generated Muons","Reconstructed Muons");
330 
331  mt.register_container<ContainerTH2<int>>("genLep_count_v_recLep_count", "Number of Generated Leptons v. Number of Reconstructed Leptons (e & mu only)",
332  lookup<std::pair<int,int>>("genLep_count_v_recLep_count"),
333  10, 0, 10, 10, 0, 10,
334  "Generated Leptons","Reconstructed Leptons");
335 
336  mt.register_container<ContainerTH1<int>>("b_jet_count", "B-Jet Multiplicity", lookup<int>("b_jet_count"), 10, 0, 10);
337 
338 
339  mt.register_container<ContainerTH1<int>>("jet_count_os_dilepton", "Jet Multiplicity - OS Dilepton Events",
340  lookup<int>("nJet"), 14, 0, 14)->add_filter(lookup_obs_filter("os-dilepton"));
341  mt.register_container<ContainerTH1<int>>("jet_count_ss_dilepton", "Jet Multiplicity - SS Dilepton Events",
342  lookup<int>("nJet"), 14, 0, 14)->add_filter(lookup_obs_filter("ss-dilepton"));
343  mt.register_container<ContainerTH1<int>>("jet_count_trilepton", "Jet Multiplicity - Trilepton Events",
344  lookup<int>("nJet"), 14, 0, 14)->add_filter(lookup_obs_filter("trilepton"));
345 
346  mt.register_container<ContainerTH1<int>>("nVert", "Number of Primary Vertices", lookup<int>("nVert"), 50, 0, 50);
347 
348  mt.register_container<CounterMany<int>>("GenPart_pdgId_counter", lookup<vector<int>>("GenPart_pdgId"));
349 
350 
351  mt.register_container<Vector<std::vector< int >>>("GenPart_pdgId", lookup<std::vector< int >>("GenPart_pdgId"));
352  mt.register_container<Vector<std::vector< int >>>("GenPart_motherIndex", lookup<std::vector< int >>("GenPart_motherIndex"));
353  mt.register_container<Vector<std::vector< int >>>("GenPart_motherId", lookup<std::vector< int >>("GenPart_motherId"));
354  mt.register_container<Vector<std::vector<float>>>("GenPart_pt", lookup<std::vector<float>>("GenPart_pt"));
355  mt.register_container<Vector<std::vector<float>>>("GenPart_eta", lookup<std::vector<float>>("GenPart_eta"));
356  mt.register_container<Vector<std::vector<float>>>("GenPart_phi", lookup<std::vector<float>>("GenPart_phi"));
357  mt.register_container<Vector<std::vector<float>>>("GenPart_mass", lookup<std::vector<float>>("GenPart_mass"));
358  mt.register_container<Vector<std::vector<float>>>("GenPart_energy", lookup<std::vector<float>>("GenPart_energy"));
359  mt.register_container<Vector<std::vector< int >>>("GenPart_status", lookup<std::vector< int >>("GenPart_status"));
360 
361  mt.register_container<Vector<vector< int >>>("LepGood_mcMatchId", lookup<vector< int >>("LepGood_mcMatchId"));
362  mt.register_container<Vector<vector< int >>>("LepGood_mcMatchPdgId", lookup<vector< int >>("LepGood_mcMatchPdgId"));
363 
364  mt.register_container<Vector< int >>("run", lookup< int >("run") );
365  mt.register_container<Vector< int >>("lumi", lookup< int >("lumi"));
366  mt.register_container<Vector< int >>("evt", lookup< int >("evt") );
367  mt.register_container<Vector<float>>("xsec", lookup<float>("xsec"));
368 
369  mt.register_container<Vector<std::vector< int >>>("Jet_mcMatchFlav", lookup<std::vector< int >>("Jet_mcMatchFlav"));
370  mt.register_container<Vector<std::vector< int >>>("Jet_mcMatchId", lookup<std::vector< int >>("Jet_mcMatchId"));
371  mt.register_container<Vector<std::vector< int >>>("Jet_mcFlavour", lookup<std::vector< int >>("Jet_mcFlavour"));
372 
373  mt.register_container<Vector<std::vector<float>>>("Jet_pt", lookup<std::vector<float>>("Jet_pt"));
374  mt.register_container<Vector<std::vector<float>>>("Jet_eta", lookup<std::vector<float>>("Jet_eta"));
375  mt.register_container<Vector<std::vector<float>>>("Jet_phi", lookup<std::vector<float>>("Jet_phi"));
376  mt.register_container<Vector<std::vector<float>>>("Jet_mass", lookup<std::vector<float>>("Jet_mass"));
377 
378 }
379 
380 
381 void run_analysis(const std::string& input_filename, bool silent){
382  gSystem->Load("libfilval.so");
383  auto replace_suffix = [](const std::string& input, const std::string& new_suffix){
384  return input.substr(0, input.find_last_of(".")) + new_suffix;
385  };
386  string log_filename = replace_suffix(input_filename, "_result.log");
387  fv::util::Log::init_logger(log_filename, fv::util::LogPriority::kLogDebug);
388 
389  string output_filename = replace_suffix(input_filename, "_result.root");
390  MiniTreeDataSet mt(input_filename, output_filename);
391 
392  enable_branches(mt);
393  declare_values(mt);
394  declare_containers(mt);
395 
396  mt.process(silent);
397  mt.save_all();
398 }
399 
400 int main(int argc, char * argv[])
401 {
402  fv::util::ArgParser args(argc, argv);
403  if(!args.cmdOptionExists("-f")) {
404  cout << "Usage: ./main (-s) -f input_minitree.root" << endl;
405  return -1;
406  }
407  bool silent = args.cmdOptionExists("-s");
408  string input_filename = args.getCmdOption("-f");
409  run_analysis(input_filename, silent);
410 }
fv::root::CounterMany
Same as Counter but accepts multiple values per fill.
Definition: container.hpp:298
std
STL namespace.
fv
The namespace containing all filval classes and functions.
Definition: api.hpp:46
fv::root
Definition: api.hpp:8
Generated by   doxygen 1.8.13