TTTT_Analysis.cpp 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382
  1. /**
  2. * @file
  3. * @author Caleb Fangmeier <caleb@fangmeier.tech>
  4. * @version 0.1
  5. *
  6. * @section LICENSE
  7. *
  8. *
  9. * MIT License
  10. *
  11. * Copyright (c) 2017 Caleb Fangmeier
  12. *
  13. * Permission is hereby granted, free of charge, to any person obtaining a copy
  14. * of this software and associated documentation files (the "Software"), to deal
  15. * in the Software without restriction, including without limitation the rights
  16. * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  17. * copies of the Software, and to permit persons to whom the Software is
  18. * furnished to do so, subject to the following conditions:
  19. *
  20. * The above copyright notice and this permission notice shall be included in all
  21. * copies or substantial portions of the Software.
  22. *
  23. * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  24. * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  25. * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  26. * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  27. * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  28. * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  29. * SOFTWARE.
  30. *
  31. * @section DESCRIPTION
  32. * Main analysis routine file. This file declares the Histogram/Graph objects
  33. * that will end up in the final root file. It also declares the values that
  34. * are used to populate the histogram, as well as how these values are
  35. * calculated. See the Fil-Val documentation for how the system works.
  36. */
  37. #include <iostream>
  38. #include <vector>
  39. #include <utility>
  40. #include <numeric>
  41. #include <limits>
  42. #include "filval/filval.hpp"
  43. #include "filval/root/root_filval.hpp"
  44. #include "analysis/common/obj_types.hpp"
  45. #include "analysis/common/constants.hpp"
  46. #include "analysis/selection.hpp"
  47. #include "MiniTreeDataSet.hpp"
  48. #include <TSystem.h>
  49. using namespace std;
  50. using namespace fv;
  51. using namespace fv::root;
  52. void enable_extra_branches(MiniTreeDataSet& mt){
  53. mt.track_branch<int>("nBJetLoose40");
  54. mt.track_branch<int>("nBJetMedium40");
  55. mt.track_branch<int>("nBJetTight40");
  56. mt.track_branch< int >("run" );
  57. mt.track_branch< int >("lumi");
  58. mt.track_branch< int >("evt" );
  59. mt.track_branch<float>("xsec");
  60. }
  61. void declare_values(MiniTreeDataSet& mt){
  62. // Define a couple selections to be used in the top-mass reconstruction.
  63. auto& b_pdgid_filter = GenFunction::reg_func<bool(Particle)>("b_pdgid_filter",
  64. FUNC(([](const Particle& j){
  65. return j.genpart.pdgId == 5 || j.genpart.pdgId==-5;
  66. })));
  67. auto& w_mass_filter = GenFunction::reg_func<bool(Particle, Particle)>("w_mass_filter",
  68. FUNC(([win_l=W_MASS-10, win_h=W_MASS+10](const Particle& j1, const Particle& j2){
  69. float inv_mass = (j1.v + j2.v).M();
  70. return inv_mass > win_l && inv_mass < win_h;
  71. })));
  72. auto& dup_filter = GenFunction::reg_func<bool(std::tuple<Particle,Particle>,Particle)>("dup_filter",
  73. FUNC(([](const std::tuple<Particle,Particle>& w, const Particle& b){
  74. int j0 = b.idx;
  75. int j1 = std::get<0>(w).idx;
  76. int j2 = std::get<1>(w).idx;
  77. return (j0 != j1) && (j0 != j2) && (j1 != j2);
  78. })));
  79. auto& qg_id_filter = GenFunction::reg_func<bool(Particle, Particle)>("qg_id_filter",
  80. FUNC(([](const Particle& j1, const Particle& j2){
  81. // require both particles be either quarks(not Top) or gluons
  82. int id1 = abs(j1.genpart.pdgId);
  83. int id2 = abs(j2.genpart.pdgId);
  84. return ((id1 >=1 && id1 <= 5) || id1 == 21) &&
  85. ((id2 >=1 && id2 <= 5) || id2 == 21);
  86. })));
  87. auto jets = lookup<std::vector<Particle>>("jets");
  88. auto leptons = lookup<std::vector<Particle>>("leptons");
  89. auto leading_jet = fv::reduce(leading_particle, jets, "leading_jet");
  90. fv::apply(particle_pt, fv::tuple(leading_jet), "leading_jet_pt");
  91. auto leading_lepton = fv::reduce(leading_particle, leptons);
  92. fv::apply(particle_pt, fv::tuple(leading_lepton), "leading_lepton_pt");
  93. fv::apply(lepton_relIso, fv::tuple(leading_lepton), "leading_lepton_relIso");
  94. auto b_jets = lookup<std::vector<Particle>>("b_jets");
  95. // Here is the calculation of the Top Reconstructed Mass from Particle
  96. auto w_dijets = tup_filter<Particle,Particle>(w_mass_filter, combinations<Particle,2>(jets, "reco_dijets"));
  97. auto top_cands = cart_product<std::tuple<Particle,Particle>, Particle>(w_dijets, b_jets);
  98. top_cands = tup_filter(dup_filter, top_cands);
  99. auto& t_mass = GenFunction::reg_func<float(std::tuple<Particle,Particle>,Particle)>("t_mass",
  100. FUNC(([](const std::tuple<Particle,Particle>& w, const Particle& b){
  101. return (std::get<0>(w).v+std::get<1>(w).v+b.v).M();
  102. })));
  103. fv::map(t_mass, top_cands, "reco_top_mass");
  104. // Here is the calculation of the Top Reconstructed Mass from Generator-Level objects
  105. jets = lookup<std::vector<Particle>>("mc_jets");
  106. b_jets = filter(b_pdgid_filter, jets);
  107. w_dijets = tup_filter(qg_id_filter, combinations<Particle,2>(jets));
  108. w_dijets = tup_filter(w_mass_filter, w_dijets);
  109. top_cands = cart_product<std::tuple<Particle,Particle>, Particle>(w_dijets, b_jets);
  110. top_cands = tup_filter(dup_filter, top_cands);
  111. fv::map(t_mass, top_cands, "mc_top_mass");
  112. // calculation of di-jet inv-mass spectrum
  113. auto& inv_mass2 = GenFunction::reg_func<float(Particle, Particle)>("inv_mass2",
  114. FUNC(([] (const Particle& j1, const Particle& j2){
  115. TLorentzVector sum = j1.v + j2.v;
  116. return (float)sum.M();
  117. })));
  118. fv::map(inv_mass2, lookup<std::vector<std::tuple<Particle,Particle>>>("reco_dijets"), "dijet_inv_mass");
  119. count<float>(GenFunction::reg_func<bool(float)>("bJet_Selection",
  120. FUNC(([](float x){
  121. return x>0;
  122. }))), "Jet_btagCMVA", "b_jet_count");
  123. auto &is_electron = GenFunction::reg_func<bool(int)>("is_electron",
  124. FUNC(([](int pdgId) {
  125. return abs(pdgId) == 11;
  126. })));
  127. auto &is_muon = GenFunction::reg_func<bool(int)>("is_muon",
  128. FUNC(([](int pdgId) {
  129. return abs(pdgId) == 13;
  130. })));
  131. auto &is_lepton = GenFunction::reg_func<bool(int)>("is_lepton",
  132. FUNC(([ie=&is_electron, im=&is_muon](int pdgId) {
  133. return (*ie)(pdgId) || (*im)(pdgId);
  134. })));
  135. count<int>(is_electron, "GenPart_pdgId", "genEle_count");
  136. count<int>(is_muon, "GenPart_pdgId", "genMu_count");
  137. count<int>(is_lepton, "GenPart_pdgId", "genLep_count");
  138. count<int>(is_electron, "LepGood_pdgId", "recEle_count");
  139. count<int>(is_muon, "LepGood_pdgId", "recMu_count");
  140. count<int>(is_lepton, "LepGood_pdgId", "recLep_count");
  141. fv::pair<int, int>("genEle_count", "recEle_count", "genEle_count_v_recEle_count");
  142. fv::pair<int, int>("genMu_count", "recMu_count", "genMu_count_v_recMu_count");
  143. fv::pair<int, int>("genLep_count", "recLep_count", "genLep_count_v_recLep_count");
  144. /* auto& sum = GenFunction::reg_func<float(std::vector<float>)>("sum", */
  145. /* FUNC(([](const std::vector<float>& v){ */
  146. /* return std::accumulate(v.begin(), v.end(), 0); */
  147. /* }))); */
  148. /* auto sum_jet_pt = fv::apply(sum, lookup<std::vector<float>>("Jet_pt")); */
  149. obs_filter("trilepton", FUNC(([nLepGood=lookup<int>("nLepGood")]()
  150. {
  151. return dynamic_cast<Value<int>*>(nLepGood)->get_value() == 3;
  152. })));
  153. obs_filter("os-dilepton", FUNC(([LepGood_charge=lookup<vector<int>>("LepGood_charge")]()
  154. {
  155. auto& charges = LepGood_charge->get_value();
  156. return charges.size()==2 && (charges[0] != charges[1]);
  157. })));
  158. obs_filter("ss-dilepton", FUNC(([LepGood_charge=lookup<vector<int>>("LepGood_charge")]()
  159. {
  160. auto& charges = LepGood_charge->get_value();
  161. return charges.size()==2 && (charges[0] == charges[1]);
  162. })));
  163. }
  164. void declare_containers(MiniTreeDataSet& mt){
  165. mt.register_container<ContainerTH1<int>>("lepton_count", "Lepton Multiplicity", lookup<int>("nLepGood"), 8, 0, 8);
  166. mt.register_container<ContainerTH1Many<float>>("Jet_pt_dist", "Jet P_T",
  167. lookup<vector<float>>("Jet_pt"), 50, 0, 500,
  168. "Jet P_T");
  169. mt.register_container<ContainerTH1Many<float>>("Jet_eta_dist", "Jet Eta",
  170. lookup<vector<float>>("Jet_eta"), 50, -3, 3,
  171. "Jet Eta");
  172. mt.register_container<ContainerTH1Many<float>>("Jet_phi_dist", "Jet Phi",
  173. lookup<vector<float>>("Jet_phi"), 20, -PI, PI,
  174. "Jet Phi");
  175. mt.register_container<ContainerTH1Many<float>>("Jet_mass_dist", "Jet Mass",
  176. lookup<vector<float>>("Jet_mass"), 50, 0, 200,
  177. "Jet Mass");
  178. mt.register_container<ContainerTH1Many<float>>("dijet_inv_mass", "Di-Jet Inv. Mass - All",
  179. lookup<vector<float>>("dijet_inv_mass"), 100, 0, 500,
  180. "Di-Jet Mass");
  181. mt.register_container<ContainerTH1Many<float>>("dijet_inv_mass_osdilepton", "Di-Jet Inv. Mass - OS Dilepton",
  182. lookup<vector<float>>("dijet_inv_mass"), 100, 0, 500,
  183. "Di-Jet Mass")->add_filter(lookup_obs_filter("os-dilepton"));
  184. mt.register_container<ContainerTH1Many<float>>("dijet_inv_mass_ssdilepton", "Di-Jet Inv. Mass - SS Dilepton",
  185. lookup<vector<float>>("dijet_inv_mass"), 100, 0, 500,
  186. "Di-Jet Mass")->add_filter(lookup_obs_filter("ss-dilepton"));
  187. mt.register_container<ContainerTH1Many<float>>("dijet_inv_mass_trilepton", "Di-Jet Inv. Mass - Trilepton",
  188. lookup<vector<float>>("dijet_inv_mass"), 100, 0, 500,
  189. "Di-Jet Mass")->add_filter(lookup_obs_filter("trilepton"));
  190. mt.register_container<ContainerTH1Many<float>>("reco_top_mass", "Reconstructed Top mass",
  191. lookup<vector<float>>("reco_top_mass"), 100, 0, 500,
  192. "Tri-jet invarient mass");
  193. /* mt.register_container<ContainerTH1Many<float>>("reco_top_pt", "Reconstructed Top Pt", */
  194. /* lookup<vector<float>>("reco_top_pt"), 100, 0, 500, */
  195. /* "Tri-jet Pt"); */
  196. mt.register_container<ContainerTH1Many<float>>("mc_top_mass", "Reconstructed MC Top mass",
  197. lookup<vector<float>>("mc_top_mass"), 100, 0, 500,
  198. "Tri-jet invarient mass");
  199. mt.register_container<ContainerTH2<int>>("nLepvsnJet", "Number of Leptons vs Number of Jets",
  200. fv::pair<int, int>("nLepGood", "nJet"),
  201. 7, 0, 7, 20, 0, 20,
  202. "Number of Leptons", "Number of Jets");
  203. mt.register_container<ContainerTH2<int>>("genEle_count_v_recEle_count", "Number of Generated Electrons v. Number of Reconstructed Electrons",
  204. lookup<std::pair<int,int>>("genEle_count_v_recEle_count"),
  205. 10, 0, 10, 10, 0, 10,
  206. "Generated Electrons","Reconstructed Electrons");
  207. mt.register_container<ContainerTH2<int>>("genMu_count_v_recMu_count", "Number of Generated Muons v. Number of Reconstructed Muons",
  208. lookup<std::pair<int,int>>("genMu_count_v_recMu_count"),
  209. 10, 0, 10, 10, 0, 10,
  210. "Generated Muons","Reconstructed Muons");
  211. mt.register_container<ContainerTH2<int>>("genLep_count_v_recLep_count", "Number of Generated Leptons v. Number of Reconstructed Leptons (e & mu only)",
  212. lookup<std::pair<int,int>>("genLep_count_v_recLep_count"),
  213. 10, 0, 10, 10, 0, 10,
  214. "Generated Leptons","Reconstructed Leptons");
  215. auto jet_count = mt.register_container<ContainerTH1<int>>("jet_count", "B-Jet Multiplicity", lookup<int>("nJet"), 15, 0, 15);
  216. mt.cut_set(jet_count,{
  217. {event_selection.trilepton, "jet_count_trilepton"},
  218. {event_selection.eem_trilepton, "jet_count_eem_trilepton"},
  219. {event_selection.mme_trilepton, "jet_count_mme_trilepton"},
  220. {event_selection.b_jet3, "jet_count_b_jet3"},
  221. {event_selection.z_mass_veto, "jet_count_z_mass_veto"},
  222. {event_selection.base_sel, "jet_count_base_selection"},
  223. {event_selection.eem_base_sel, "jet_count_eem_base_selection"},
  224. {event_selection.mme_base_sel, "jet_count_mme_base_selection"},
  225. });
  226. auto b_jet_count = mt.register_container<ContainerTH1<int>>("b_jet_count", "B-Jet Multiplicity", lookup<int>("b_jet_count"), 10, 0, 10);
  227. mt.cut_set(b_jet_count,{
  228. {event_selection.trilepton, "b_jet_count_trilepton"},
  229. {event_selection.eem_trilepton, "b_jet_count_eem_trilepton"},
  230. {event_selection.mme_trilepton, "b_jet_count_mme_trilepton"},
  231. {event_selection.b_jet3, "b_jet_count_b_jet3"},
  232. {event_selection.z_mass_veto, "b_jet_count_z_mass_veto"},
  233. {event_selection.base_sel, "b_jet_count_base_selection"},
  234. {event_selection.eem_base_sel, "b_jet_count_eem_base_selection"},
  235. {event_selection.mme_base_sel, "b_jet_count_mme_base_selection"},
  236. });
  237. mt.register_container<ContainerTH1<int>>("jet_count_os_dilepton", "Jet Multiplicity - OS Dilepton Events",
  238. lookup<int>("nJet"), 14, 0, 14)->add_filter(lookup_obs_filter("os-dilepton"));
  239. mt.register_container<ContainerTH1<int>>("jet_count_ss_dilepton", "Jet Multiplicity - SS Dilepton Events",
  240. lookup<int>("nJet"), 14, 0, 14)->add_filter(lookup_obs_filter("ss-dilepton"));
  241. mt.register_container<CounterMany<int>>("GenPart_pdgId_counter", lookup<vector<int>>("GenPart_pdgId"));
  242. mt.register_container<Vector<std::vector< int >>>("GenPart_pdgId", lookup<std::vector< int >>("GenPart_pdgId"));
  243. mt.register_container<Vector<std::vector< int >>>("GenPart_motherIndex", lookup<std::vector< int >>("GenPart_motherIndex"));
  244. mt.register_container<Vector<std::vector< int >>>("GenPart_motherId", lookup<std::vector< int >>("GenPart_motherId"));
  245. mt.register_container<Vector<std::vector<float>>>("GenPart_pt", lookup<std::vector<float>>("GenPart_pt"));
  246. mt.register_container<Vector<std::vector<float>>>("GenPart_eta", lookup<std::vector<float>>("GenPart_eta"));
  247. mt.register_container<Vector<std::vector<float>>>("GenPart_phi", lookup<std::vector<float>>("GenPart_phi"));
  248. mt.register_container<Vector<std::vector<float>>>("GenPart_mass", lookup<std::vector<float>>("GenPart_mass"));
  249. mt.register_container<Vector<std::vector<float>>>("GenPart_energy", lookup<std::vector<float>>("GenPart_energy"));
  250. mt.register_container<Vector<std::vector< int >>>("GenPart_status", lookup<std::vector< int >>("GenPart_status"));
  251. mt.register_container<Vector<vector< int >>>("LepGood_mcMatchPdgId", lookup<vector< int >>("LepGood_mcMatchPdgId"));
  252. mt.register_container<Vector<vector< int >>>("LepGood_pdgId", lookup<vector< int >>("LepGood_pdgId"));
  253. mt.register_container<Vector< int >>("run", lookup< int >("run") );
  254. mt.register_container<Vector< int >>("lumi", lookup< int >("lumi"));
  255. mt.register_container<Vector< int >>("evt", lookup< int >("evt") );
  256. mt.register_container<Vector<float>>("xsec", lookup<float>("xsec"));
  257. mt.register_container<Vector<std::vector< int >>>("Jet_mcMatchFlav", lookup<std::vector< int >>("Jet_mcMatchFlav"));
  258. mt.register_container<Vector<std::vector<float>>>("Jet_pt", lookup<std::vector<float>>("Jet_pt"));
  259. mt.register_container<Vector<std::vector<float>>>("Jet_eta", lookup<std::vector<float>>("Jet_eta"));
  260. mt.register_container<Vector<std::vector<float>>>("Jet_phi", lookup<std::vector<float>>("Jet_phi"));
  261. mt.register_container<Vector<std::vector<float>>>("Jet_mass", lookup<std::vector<float>>("Jet_mass"));
  262. mt.register_container<ContainerTH1<float>>("leading_jet_pt", "Leading Jet Pt", lookup<float>("leading_jet_pt"), 100, 0, 500);
  263. mt.register_container<ContainerTH1<float>>("leading_lepton_pt", "Leading Lepton Pt", lookup<float>("leading_lepton_pt"), 100, 0, 500);
  264. mt.register_container<ContainerTH1<float>>("leading_lepton_relIso", "Leading Lepton Relative Isolation", lookup<float>("leading_lepton_relIso"), 100, 0, 0.05);
  265. mt.register_container<Vector<float>>("leading_lepton_relIso_all", lookup<float>("leading_lepton_relIso"));
  266. mt.register_container<PassCount>("trilepton_count", event_selection.trilepton);
  267. mt.register_container<PassCount>("b_jet3_count", event_selection.b_jet3);
  268. mt.register_container<PassCount>("z_mass_veto_count", event_selection.z_mass_veto);
  269. mt.register_container<PassCount>("J4_count", event_selection.J4);
  270. mt.register_container<PassCount>("J5_count", event_selection.J5);
  271. mt.register_container<PassCount>("J6_count", event_selection.J6);
  272. mt.register_container<PassCount>("SR4j_count", event_selection.SR4j);
  273. mt.register_container<PassCount>("SR5j_count", event_selection.SR5j);
  274. mt.register_container<PassCount>("SR6j_count", event_selection.SR6j);
  275. }
  276. void run_analysis(const std::string& input_filename, const std::string& data_label, bool silent){
  277. gSystem->Load("libfilval.so");
  278. auto replace_suffix = [](const std::string& input, const std::string& new_suffix){
  279. return input.substr(0, input.find_last_of(".")) + new_suffix;
  280. };
  281. string log_filename = replace_suffix(input_filename, "_result.log");
  282. fv::util::Log::init_logger(log_filename, fv::util::LogPriority::kLogDebug);
  283. string output_filename = replace_suffix(input_filename, "_result.root");
  284. MiniTreeDataSet mt(output_filename, input_filename, data_label);
  285. create_all_common_values(mt);
  286. enable_extra_branches(mt);
  287. declare_values(mt);
  288. init_selection();
  289. declare_containers(mt);
  290. mt.process(silent);
  291. mt.save_all();
  292. }
  293. int main(int argc, char * argv[])
  294. {
  295. fv::util::ArgParser args(argc, argv);
  296. if(!args.cmdOptionExists("-l") || !args.cmdOptionExists("-f")) {
  297. cout << "Usage: ./main (-s) -l DATA_LABEL -f input_minitree.root" << endl;
  298. return -1;
  299. }
  300. bool silent = args.cmdOptionExists("-s");
  301. string input_filename = args.getCmdOption("-f");
  302. string data_label = args.getCmdOption("-l");
  303. run_analysis(input_filename, data_label, silent);
  304. }