MiniTree.cpp 10.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289
  1. #ifndef minitree_cxx
  2. #define minitree_cxx
  3. #include <algorithm>
  4. #include <functional>
  5. #include <cstdio>
  6. #include <string>
  7. #include <vector>
  8. #include <map>
  9. #include <TStyle.h>
  10. #include <TLorentzVector.h>
  11. #include "MiniTree.h"
  12. #include "HistCollection.h"
  13. using namespace std;
  14. HistCollection* build_histograms(const char* sample_name, const char* filename){
  15. HistCollection *hc = new HistCollection(sample_name);
  16. TFile *f = TFile::Open(filename);
  17. TTree *tree = (TTree*) f->Get("tree");
  18. MiniTree minitree(tree);
  19. minitree.Loop(hc);
  20. delete tree;
  21. f->Close();
  22. delete f;
  23. printf("Finished filling histograms for file %s with %d events.\n",
  24. filename, hc->get_event_count());
  25. return hc;
  26. }
  27. /* HistCollection* load_histograms(const char* sample_name){ */
  28. /* string filename = string(sample_name) + "_histcollection.root"; */
  29. /* TFile *f = TFile::Open(filename.c_str()); */
  30. /* HistCollection *hc = new HistCollection(sample_name); */
  31. /* TFile *f = TFile::Open(filename); */
  32. /* TTree *tree = (TTree*) f->Get("tree"); */
  33. /* MiniTree minitree(tree); */
  34. /* minitree.Loop(hc); */
  35. /* delete tree; */
  36. /* f->Close(); */
  37. /* delete f; */
  38. /* printf("Finished filling histograms for file %s with %d events.\n", */
  39. /* filename, hc->get_event_count()); */
  40. /* return hc; */
  41. /* } */
  42. void set_branch_statuses(MiniTree *minitree){
  43. auto on = [minitree](const char* bname){
  44. minitree->fChain->SetBranchStatus(bname, true);
  45. };
  46. auto off = [minitree](const char* bname){
  47. minitree->fChain->SetBranchStatus(bname, false);
  48. };
  49. off("*");
  50. on("nLepGood");
  51. on("LepGood_pt");
  52. on("LepGood_eta");
  53. on("LepGood_phi");
  54. on("LepGood_mcPt");
  55. on("LepGood_charge");
  56. on("nJet");
  57. on("Jet_btagCMVA");
  58. on("Jet_pt");
  59. on("Jet_eta");
  60. on("Jet_phi");
  61. on("nGenTop");
  62. on("ngenLep");
  63. on("genLep_sourceId");
  64. on("genLep_pt");
  65. }
  66. double delta_r2(double eta1, double phi1, double eta2, double phi2){
  67. return pow(eta1 - eta2, 2) + pow(phi1 - phi2, 2);
  68. }
  69. double isolation(int i, int n, function<double(int, int)> delta_r2_calculator){
  70. double min_iso = -1;
  71. for (int j=0; j<n; j++){
  72. if (i == j) continue;
  73. double curr_iso = delta_r2_calculator(i, j);
  74. if ((curr_iso < min_iso && curr_iso >= 0) || min_iso == -1){
  75. min_iso = curr_iso;
  76. }
  77. }
  78. if (min_iso == -1)
  79. return min_iso;
  80. else
  81. return sqrt(min_iso);
  82. }
  83. TLorentzVector g_vec1;
  84. TLorentzVector g_vec2;
  85. double diobject_mass(double pt1, double eta1, double phi1, double mass1,
  86. double pt2, double eta2, double phi2, double mass2){
  87. g_vec1.SetPtEtaPhiM(pt1, eta1, phi1, mass1);
  88. g_vec2.SetPtEtaPhiM(pt2, eta2, phi2, mass2);
  89. g_vec1 += g_vec2;
  90. return g_vec1.M();
  91. }
  92. void MiniTree::Loop(HistCollection* hc){
  93. /* Define first a few little subroutines that will be used later in the
  94. * main loop.
  95. */
  96. // Functions to calculate the \delta R^2 between two objects
  97. auto lepton_lepton_delta_r2 = [this](int i, int j){
  98. return delta_r2(this->LepGood_eta[i], this->LepGood_phi[i],
  99. this->LepGood_eta[j], this->LepGood_phi[j]);
  100. };
  101. auto lepton_jet_delta_r2 = [this](int i, int j){
  102. return delta_r2(this->LepGood_eta[i], this->LepGood_phi[i],
  103. this->Jet_eta[j], this->Jet_phi[j]);
  104. };
  105. auto jet_jet_delta_r2 = [this](int i, int j){
  106. return delta_r2(this->Jet_eta[i], this->Jet_phi[i],
  107. this->Jet_eta[j], this->Jet_phi[j]);
  108. };
  109. // Funtions to calculate isolation of objects
  110. auto lepton_lepton_isolation = [this, lepton_lepton_delta_r2](int i){
  111. if (i >= this->nLepGood) return -1.0;
  112. return isolation(i, this->nLepGood, function<double(int, int)>(lepton_lepton_delta_r2));
  113. };
  114. auto lepton_jet_isolation = [this, lepton_jet_delta_r2](int i){
  115. if (i >= this->nLepGood) return -1.0;
  116. return isolation(i, this->nJet, function<double(int, int)>(lepton_jet_delta_r2));
  117. };
  118. auto mini_isolation_cut = [this, lepton_jet_isolation](int i){
  119. /* Implement the mini-isolation cut which utilizes a dynamic isolation cone size
  120. * based on lepton pt
  121. */
  122. double r0 = 0.5;
  123. double pt0 = 7.5; // GeV
  124. double iso_cut = min(r0, pt0/this->LepGood_pt[i]);
  125. double isolation = lepton_jet_isolation(i);
  126. return isolation > iso_cut;
  127. };
  128. auto jet_jet_isolation = [this, jet_jet_delta_r2](int i){
  129. if (i >= this->nJet) return -1.0;
  130. return isolation(i, this->nJet, function<double(int, int)>(jet_jet_delta_r2));
  131. };
  132. auto dijet_mass = [this](int i, int j){
  133. if (i >= this->nJet || j >= this->nJet) return -1.0;
  134. return diobject_mass(this->Jet_pt[i], this->Jet_eta[i],
  135. this->Jet_phi[i], this->Jet_mass[i],
  136. this->Jet_pt[j], this->Jet_eta[j],
  137. this->Jet_phi[j], this->Jet_mass[j]);
  138. };
  139. auto dilepton_mass = [this](int i, int j){
  140. if (i >= this->nLepGood || j >= this->nLepGood) return -1.0;
  141. return diobject_mass(this->LepGood_pt[i], this->LepGood_eta[i],
  142. this->LepGood_phi[i], this->LepGood_mass[i],
  143. this->LepGood_pt[j], this->LepGood_eta[j],
  144. this->LepGood_phi[j], this->LepGood_mass[j]);
  145. };
  146. auto dijet_mass_set = [this, dijet_mass](){
  147. auto *mass_set = new vector<double>();
  148. for(int i=0; i<(this->nJet-1); i++){
  149. for(int j=i+1; j<this->nJet; j++){
  150. mass_set->push_back(dijet_mass(i, j));
  151. }
  152. }
  153. return mass_set;
  154. };
  155. auto nBJets = [this](float mva_cut = 0){
  156. int n = 0;
  157. for(int i=0; i<this->nJet; ++i)
  158. if(this->Jet_btagCMVA[i] > mva_cut) ++n;
  159. return n;
  160. };
  161. auto lepton_mass_window_pass = [this, dilepton_mass](double low, double high){
  162. for (int i=0; i<this->nLepGood, ++i){
  163. for (int j=0; j<this->nLepGood, ++j){
  164. double mass = dilepton_mass(i, j);
  165. if (low < mass && mass < high){
  166. return false;
  167. }
  168. }
  169. }
  170. return true;
  171. }
  172. auto ss_dilepton_base_selection[this, nBJets, lepton_mass_window_pass](){
  173. if (this->nLepGood != 2) return false;
  174. if (this->LepGood_charge[0] != this->LepGood_charge[1]) return false;
  175. if (this->nJet < 4) return false;
  176. if (nBJets() < 3) return false;
  177. if (!lepton_mass_window_pass(70, 105)) return false;
  178. return true;
  179. }
  180. auto trilepton_base_selection[this](){
  181. if (nLepGood != 3) return false;
  182. if (nJet < 4) return false;
  183. if (nBJets() < 3) return false;
  184. if (!lepton_mass_window_pass(70, 105)) return false;
  185. return true;
  186. }
  187. if (fChain == 0) return;
  188. /* Set branch statuses so only used branches are loaded. This
  189. * can greatly reduce execution time.
  190. */
  191. set_branch_statuses(this);
  192. Long64_t nentries = fChain->GetEntriesFast();
  193. Long64_t nbytes = 0, nb = 0;
  194. for (Long64_t jentry=0; jentry<nentries;jentry++) {
  195. Long64_t ientry = LoadTree(jentry);
  196. if (ientry < 0) break;
  197. nb = fChain->GetEntry(jentry); nbytes += nb;
  198. /* Main analysis code begins here
  199. *
  200. */
  201. ++(*hc); // Increments event count
  202. hc->lepton_count->Fill(nLepGood);
  203. hc->top_quark_count->Fill(nGenTop);
  204. for(int i=0; i<nLepGood; i++){
  205. double delta_pt_val = LepGood_pt[i] - LepGood_mcPt[i];
  206. hc->delta_pt->Fill(delta_pt_val);
  207. hc->lepton_count_vs_delta_pt->Fill(nLepGood, delta_pt_val);
  208. }
  209. hc->jet_count->Fill(nJet);
  210. int b_jet_count = 0;
  211. for(int i=0; i<nJet; i++){
  212. hc->b_jet_discriminator->Fill(Jet_btagCMVA[i]);
  213. b_jet_count += (Jet_btagCMVA[i] > 0);
  214. }
  215. hc->b_jet_count->Fill(b_jet_count);
  216. for(int i=0; i<nJet; i++){
  217. hc->b_jet_pt_all->Fill(Jet_pt[i]);
  218. if (b_jet_count >= 3){
  219. hc->b_jet_pt_3_or_more->Fill(Jet_pt[i]);
  220. }
  221. }
  222. int lepton_count_pass_miniiso= 0;
  223. for(int i=0; i<nLepGood; i++){
  224. lepton_count_pass_miniiso += mini_isolation_cut(i);
  225. }
  226. hc->lepton_count_pass_miniiso->Fill(lepton_count_pass_miniiso);
  227. if (lepton_count_pass_miniiso == nLepGood){
  228. hc->lepton_count_pass_miniiso_event->Fill(nLepGood);
  229. }
  230. if (nLepGood == 3){
  231. hc->jet_count_trilepton->Fill(nJet);
  232. }
  233. if (nLepGood == 2 && (LepGood_charge[0] * LepGood_charge[1]) == 1){
  234. hc->jet_count_ss_dilepton->Fill(nJet);
  235. }
  236. if (nLepGood == 2 && (LepGood_charge[0] * LepGood_charge[1]) == -1){
  237. hc->jet_count_os_dilepton->Fill(nJet);
  238. }
  239. if (nLepGood == 3){
  240. vector<double> isos = {lepton_jet_isolation(0),
  241. lepton_jet_isolation(1),
  242. lepton_jet_isolation(2)};
  243. sort(isos.begin(), isos.end()); //Sorts ascending by default
  244. hc->trilepton_iso_1->Fill(isos[2]);
  245. hc->trilepton_iso_2->Fill(isos[1]);
  246. hc->trilepton_iso_3->Fill(isos[0]);
  247. vector<double> jet_isos = {jet_jet_isolation(0),
  248. jet_jet_isolation(1),
  249. jet_jet_isolation(2)};
  250. sort(jet_isos.begin(), jet_isos.end());
  251. hc->trilepton_jet_iso_1->Fill(jet_isos[2]);
  252. hc->trilepton_jet_iso_2->Fill(jet_isos[1]);
  253. hc->trilepton_jet_iso_3->Fill(jet_isos[0]);
  254. }
  255. vector<double> *mass_set = dijet_mass_set();
  256. for( double &mass : *mass_set ){
  257. hc->dijet_mass->Fill(mass);
  258. }
  259. }
  260. hc->lepton_count_pass_miniiso_ratio->Divide(hc->lepton_count_pass_miniiso_event,
  261. hc->lepton_count);
  262. }
  263. #endif // tree_cxx