LOGO
RestFrames  v1.0.0
RestFrames HEP Event Analysis Software Library
example_H_to_ttbar_to_bWlnubWlnu.C
Go to the documentation of this file.
1 // RestFrames: particle physics event analysis library
3 // --------------------------------------------------------------------
4 // Copyright (c) 2014-2016, Christopher Rogan
13 //
14 // This file is part of RestFrames.
15 //
16 // RestFrames is free software; you can redistribute it and/or modify
17 // it under the terms of the GNU General Public License as published by
18 // the Free Software Foundation; either version 2 of the License, or
19 // (at your option) any later version.
20 //
21 // RestFrames is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with RestFrames. If not, see <http://www.gnu.org/licenses/>.
29 
30 #define COMPILER (!defined(__CINT__) && !defined(__CLING__))
31 #if defined(__MAKECINT__) || defined(__ROOTCLING__) || COMPILER
32 #include "RestFrames/RestFrames.hh"
33 #else
34 RestFrames::RFKey ensure_autoload(1);
35 #endif
36 
37 using namespace RestFrames;
38 
39 void example_H_to_ttbar_to_bWlnubWlnu(std::string output_name =
40  "output_H_to_ttbar_to_bWlnubWlnu.root"){
41 
42  double mT = 173.21; // GeV, PDG 2016
43  double mW = 80.385;
44  double mB = 4.18;
45  double mL = 0.106;
46  double mN = 0.;
47 
48  std::vector<double> mH; // vary neutral Higgs mass
49  mH.push_back(400.);
50  mH.push_back(750.);
51  mH.push_back(1000.);
52  mH.push_back(1500.);
53 
54  // number of events to generate (per Higgs mass)
55  int Ngen = 10000;
56 
58  g_Log << LogInfo << "Initializing generator frames and tree..." << LogEnd;
60  ppLabGenFrame LAB_Gen("LAB_Gen","LAB");
61  DecayGenFrame H_Gen("H_Gen","H^{ 0}");
62  DecayGenFrame Ta_Gen("Ta_Gen","t_{a}");
63  DecayGenFrame Tb_Gen("Tb_Gen","t_{b}");
64  DecayGenFrame Wa_Gen("Wa_Gen","W_{a}");
65  DecayGenFrame Wb_Gen("Wb_Gen","W_{b}");
66  VisibleGenFrame Ba_Gen("Ba_Gen","b_{a}");
67  VisibleGenFrame La_Gen("La_Gen","#it{l}_{a}");
68  InvisibleGenFrame Na_Gen("Na_Gen","#nu_{a}");
69  VisibleGenFrame Bb_Gen("Bb_Gen","b_{b}");
70  VisibleGenFrame Lb_Gen("Lb_Gen","#it{l}_{b}");
71  InvisibleGenFrame Nb_Gen("Nb_Gen","#nu_{b}");
72 
73  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
74 
75  LAB_Gen.SetChildFrame(H_Gen);
76  H_Gen.AddChildFrame(Ta_Gen);
77  H_Gen.AddChildFrame(Tb_Gen);
78  Ta_Gen.AddChildFrame(Ba_Gen);
79  Ta_Gen.AddChildFrame(Wa_Gen);
80  Tb_Gen.AddChildFrame(Bb_Gen);
81  Tb_Gen.AddChildFrame(Wb_Gen);
82  Wa_Gen.AddChildFrame(La_Gen);
83  Wa_Gen.AddChildFrame(Na_Gen);
84  Wb_Gen.AddChildFrame(Lb_Gen);
85  Wb_Gen.AddChildFrame(Nb_Gen);
86 
87  if(LAB_Gen.InitializeTree())
88  g_Log << LogInfo << "...Successfully initialized generator tree" << LogEnd;
89  else
90  g_Log << LogError << "...Failed initializing generator tree" << LogEnd;
91 
92  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
93 
94  if(mH.size() < 1) return;
95 
96  // Higgs mass
97  H_Gen.SetMass(mH[0]);
98  // set top masses
99  Ta_Gen.SetMass(mT); Tb_Gen.SetMass(mT);
100  // set W masses
101  Wa_Gen.SetMass(mW); Wb_Gen.SetMass(mW);
102  // set B masses
103  Ba_Gen.SetMass(mB); Bb_Gen.SetMass(mB);
104  // set : masses
105  La_Gen.SetMass(mL); Lb_Gen.SetMass(mL);
106  // set neutrino masses
107  Na_Gen.SetMass(mN); Nb_Gen.SetMass(mN);
108 
109  // set b-jet/lepton pT/eta cuts
110  Ba_Gen.SetPtCut(20.); Bb_Gen.SetPtCut(20.);
111  Ba_Gen.SetEtaCut(2.5); Bb_Gen.SetEtaCut(2.5);
112  La_Gen.SetPtCut(15.); Lb_Gen.SetPtCut(15.);
113  La_Gen.SetEtaCut(2.5); Lb_Gen.SetEtaCut(2.5);
114 
115  if(LAB_Gen.InitializeAnalysis())
116  g_Log << LogInfo << "...Successfully initialized generator analysis" << LogEnd;
117  else
118  g_Log << LogError << "...Failed initializing generator analysis" << LogEnd;
122  g_Log << LogInfo << "Initializing reconstruction frames and tree..." << LogEnd;
124 
125  LabRecoFrame LAB("LAB","LAB");
126  DecayRecoFrame H("H","H^{ 0}");
127  DecayRecoFrame Ta("Ta","t_{a}");
128  DecayRecoFrame Tb("Tb","t_{b}");
129  DecayRecoFrame Wa("Wa","W_{a}");
130  DecayRecoFrame Wb("Wb","W_{b}");
131  VisibleRecoFrame Ba("Ba","b_{a}");
132  VisibleRecoFrame La("La","#it{l}_{a}");
133  InvisibleRecoFrame Na("Na","#nu_{a}");
134  VisibleRecoFrame Bb("Bb","b_{b}");
135  VisibleRecoFrame Lb("Lb","#it{l}_{b}");
136  InvisibleRecoFrame Nb("Nb","#nu_{b}");
137 
138  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
139 
140  LAB.SetChildFrame(H);
141  H.AddChildFrame(Ta);
142  H.AddChildFrame(Tb);
143  Ta.AddChildFrame(Ba);
144  Ta.AddChildFrame(Wa);
145  Tb.AddChildFrame(Bb);
146  Tb.AddChildFrame(Wb);
147  Wa.AddChildFrame(La);
148  Wa.AddChildFrame(Na);
149  Wb.AddChildFrame(Lb);
150  Wb.AddChildFrame(Nb);
151 
152  if(LAB.InitializeTree())
153  g_Log << LogInfo << "...Successfully initialized reconstruction tree" << LogEnd;
154  else
155  g_Log << LogError << "...Failed initializing reconstruction tree" << LogEnd;
156 
157  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
158 
160  std::string group_name;
161 
162  // Invisible Group
163  group_name = "#nu #nu Jigsaws";
164  InvisibleGroup INV("INV", group_name);
165  INV.AddFrame(Na);
166  INV.AddFrame(Nb);
167  CombinatoricGroup B("VIS","b-jet Jigsaws");
168  B.AddFrame(Ba);
169  B.AddFrame(Bb);
170  B.SetNElementsForFrame(Ba, 1);
171  B.SetNElementsForFrame(Bb, 1);
172 
174  std::string jigsaw_name;
175 
176  // Minimize difference Mt jigsaws
177  jigsaw_name = "M_{#nu #nu} ~ m_{#it{l} #it{l}}";
178  SetMassInvJigsaw NuNuM("NuNuM", jigsaw_name);
179  INV.AddJigsaw(NuNuM);
180 
181  jigsaw_name = "#eta_{#nu#nu} = #eta_{b #it{l} b #it{l}}";
182  SetRapidityInvJigsaw NuNuR("NuNuR", jigsaw_name);
183  INV.AddJigsaw(NuNuR);
184  NuNuR.AddVisibleFrames(H.GetListVisibleFrames());
185 
186  jigsaw_name = "min ( M_{top a}- M_{top b} )^{2}";
187  MinMassDiffInvJigsaw MinDeltaMt("MinDeltaMt", jigsaw_name, 2);
188  INV.AddJigsaw(MinDeltaMt);
189  MinDeltaMt.AddInvisibleFrame(Na, 0);
190  MinDeltaMt.AddInvisibleFrame(Nb, 1);
191  MinDeltaMt.AddVisibleFrames(La+Ba, 0);
192  MinDeltaMt.AddVisibleFrames(Lb+Bb, 1);
193  MinDeltaMt.AddMassFrame(La, 0);
194  MinDeltaMt.AddMassFrame(Lb, 1);
195 
196  // b-jet combinatoric jigsaws for all trees
197  jigsaw_name = "Minimize M(b #it{l} )_{a} , M(b #it{l} )_{b}";
198 
199  MinMassesCombJigsaw MinBL("MinBL", jigsaw_name);
200  B.AddJigsaw(MinBL);
201  MinBL.AddFrames(La+Ba,0);
202  MinBL.AddFrames(Lb+Bb,1);
203 
204  if(LAB.InitializeAnalysis())
205  g_Log << LogInfo << "...Successfully initialized analysis" << LogEnd;
206  else
207  g_Log << LogError << "...Failed initializing analysis" << LogEnd;
208 
211 
212  TreePlot* treePlot = new TreePlot("TreePlot","TreePlot");
213 
214  treePlot->SetTree(LAB_Gen);
215  treePlot->Draw("GenTree", "Generator Tree", true);
216 
217  treePlot->SetTree(LAB);
218  treePlot->Draw("RecoTree", "Reconstruction Tree");
219 
220  treePlot->SetTree(B);
221  treePlot->Draw("VisTree", "b-jet Jigsaws", true);
222 
223  treePlot->SetTree(INV);
224  treePlot->Draw("InvTree", "Inivisible Jigsaws", true);
225 
226  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
227 
228  std::string plot_title = "H^{ 0} #rightarrow t #bar{t} #rightarrow b W(#it{l} #nu) b W(#it{l} #nu)";
229  HistPlot* histPlot = new HistPlot("HistPlot", plot_title);
230 
231  RFList<const HistPlotCategory> cat_list;
232  int Nmass = mH.size();
233  for(int m = 0; m < Nmass; m++){
234  char smass[50], scat[50];
235  sprintf(scat, "MH%.0f", mH[m]);
236  sprintf(smass, "m_{H^{ 0}} = %.0f", mH[m]);
237  cat_list += histPlot->GetNewCategory(scat, smass);
238  }
239 
240  const HistPlotVar& MH = histPlot->GetNewVar("MH", "M_{H^{ 0}} / m_{H^{ 0}}^{ true}", 0., 2.);
241  const HistPlotVar& Eb_ta = histPlot->GetNewVar("Eb_ta", "E_{b a}^{top a} / E_{b a}^{top a gen}", 0., 2.);
242  const HistPlotVar& Eb_tb = histPlot->GetNewVar("Eb_tb", "E_{b b}^{top b} / E_{b b}^{top b gen}", 0., 2.);
243  const HistPlotVar& El_Wa = histPlot->GetNewVar("El_Wa", "E_{#it{l} a}^{W a} / E_{#it{l} a}^{W a true}", 0., 2.);
244  const HistPlotVar& El_Wb = histPlot->GetNewVar("El_Wb", "E_{#it{l} b}^{W b} / E_{#it{l} b}^{W b true}", 0., 2.);
245  const HistPlotVar& cosH = histPlot->GetNewVar("cosH","cos #theta_{H^{ 0}}", -1., 1.);
246  const HistPlotVar& costa = histPlot->GetNewVar("costa","cos #theta_{top a}", -1., 1.);
247  const HistPlotVar& costb = histPlot->GetNewVar("costb","cos #theta_{top b}", -1., 1.);
248  const HistPlotVar& cosWa = histPlot->GetNewVar("cosWa","cos #theta_{W a}", -1., 1.);
249  const HistPlotVar& cosWb = histPlot->GetNewVar("cosWb","cos #theta_{W b}", -1., 1.);
250  const HistPlotVar& DcosH = histPlot->GetNewVar("Dcostt","#theta_{H^{ 0}} - #theta_{H^{ 0}}^{true}",
251  -acos(-1.)/2., acos(-1.)/2.);
252  const HistPlotVar& Dcosta = histPlot->GetNewVar("Dcosta","#theta_{top a} - #theta_{top a}^{true}",
253  -acos(-1.)/2., acos(-1.)/2.);
254  const HistPlotVar& Dcostb = histPlot->GetNewVar("Dcostb","#theta_{top b} - #theta_{top b}^{true}",
255  -acos(-1.)/2., acos(-1.)/2.);
256  const HistPlotVar& DcosWa = histPlot->GetNewVar("DcosWa","#theta_{W a} - #theta_{W a}^{true}",
257  -acos(-1.)/2., acos(-1.)/2.);
258  const HistPlotVar& DcosWb = histPlot->GetNewVar("DcosWb","#theta_{W b} - #theta_{W b}^{true}",
259  -acos(-1.)/2., acos(-1.)/2.);
260 
261  histPlot->AddPlot(MH, cat_list);
262  histPlot->AddPlot(Eb_ta, cat_list);
263  histPlot->AddPlot(El_Wa, cat_list);
264  histPlot->AddPlot(DcosH, cat_list);
265  histPlot->AddPlot(Dcosta, cat_list);
266  histPlot->AddPlot(DcosWa, cat_list);
267 
268  histPlot->AddPlot(MH, Eb_ta, cat_list[1]);
269  histPlot->AddPlot(MH, El_Wa, cat_list[1]);
270  histPlot->AddPlot(Eb_ta, Eb_tb, cat_list[1]);
271  histPlot->AddPlot(El_Wa, El_Wb, cat_list[1]);
272  histPlot->AddPlot(Eb_ta, El_Wa, cat_list[1]);
273  histPlot->AddPlot(DcosH, MH, cat_list[1]);
274  histPlot->AddPlot(Dcosta, Eb_ta, cat_list[1]);
275  histPlot->AddPlot(DcosWa, El_Wa, cat_list[1]);
276  histPlot->AddPlot(DcosH, Dcosta, cat_list[1]);
277  histPlot->AddPlot(Dcosta, Dcostb, cat_list[1]);
278  histPlot->AddPlot(DcosWa, DcosWb, cat_list[1]);
279  histPlot->AddPlot(Dcosta, DcosWa, cat_list[1]);
280 
283 
284  for(int m = 0; m < Nmass; m++){
285  g_Log << LogInfo << "Generating events for H^{0} mass = " << mH[m] << LogEnd;
286 
287  H_Gen.SetMass(mH[m]);
288  LAB_Gen.InitializeAnalysis();
289 
290  for(int igen = 0; igen < Ngen; igen++){
291  if(igen%((std::max(Ngen,10))/10) == 0)
292  g_Log << LogInfo << "Generating event " << igen << " of " << Ngen << LogEnd;
293 
294  // generate event
295  LAB_Gen.ClearEvent(); // clear the gen tree
296 
297  LAB_Gen.AnalyzeEvent(); // generate a new event
298 
299  TVector3 MET = LAB_Gen.GetInvisibleMomentum(); // Get the MET from gen tree
300  MET.SetZ(0.);
301 
302  // analyze event
303  LAB.ClearEvent(); // clear the event
304 
305  INV.SetLabFrameThreeVector(MET); // Set MET in reco tree
306 
307  La.SetLabFrameFourVector(La_Gen.GetFourVector()); // Set lepton 4-vectors
308  Lb.SetLabFrameFourVector(Lb_Gen.GetFourVector());
309 
310  B.AddLabFrameFourVector(Ba_Gen.GetFourVector()); // Set b-jet 4-vectors
311  B.AddLabFrameFourVector(Bb_Gen.GetFourVector());
312 
313  LAB.AnalyzeEvent(); // analyze the event
314 
316  // Observable Calculations
318 
319  double MHgen = H_Gen.GetMass();
320  double Eb_tagen = Ba_Gen.GetFourVector(Ta_Gen).E();
321  double Eb_tbgen = Bb_Gen.GetFourVector(Tb_Gen).E();
322  double El_Wagen = La_Gen.GetFourVector(Wa_Gen).E();
323  double El_Wbgen = Lb_Gen.GetFourVector(Wb_Gen).E();
324  double cosHgen = H_Gen.GetCosDecayAngle();
325  double costagen = Ta_Gen.GetCosDecayAngle();
326  double costbgen = Tb_Gen.GetCosDecayAngle();
327  double cosWagen = Wa_Gen.GetCosDecayAngle();
328  double cosWbgen = Wb_Gen.GetCosDecayAngle();
329 
330  MH = H.GetMass()/MHgen;
331  Eb_ta = Ba.GetFourVector(Ta).E()/Eb_tagen;
332  Eb_tb = Bb.GetFourVector(Tb).E()/Eb_tbgen;
333  El_Wa = La.GetFourVector(Wa).E()/El_Wagen;
334  El_Wb = Lb.GetFourVector(Wb).E()/El_Wbgen;
335  cosH = H.GetCosDecayAngle();
336  costa = Ta.GetCosDecayAngle();
337  costb = Tb.GetCosDecayAngle();
338  cosWa = Wa.GetCosDecayAngle();
339  cosWb = Wb.GetCosDecayAngle();
340  DcosH = asin(sqrt(1.-cosH*cosH)*cosHgen-sqrt(1.-cosHgen*cosHgen)*cosH);
341  Dcosta = asin(sqrt(1.-costa*costa)*costagen-sqrt(1.-costagen*costagen)*costa);
342  Dcostb = asin(sqrt(1.-costb*costb)*costbgen-sqrt(1.-costbgen*costbgen)*costb);
343  DcosWa = asin(sqrt(1.-cosWa*cosWa)*cosWagen-sqrt(1.-cosWagen*cosWagen)*cosWa);
344  DcosWb = asin(sqrt(1.-cosWb*cosWb)*cosWbgen-sqrt(1.-cosWbgen*cosWbgen)*cosWb);
345 
346  histPlot->Fill(cat_list[m]);
347  }
348 
349  LAB_Gen.PrintGeneratorEfficiency();
350  }
351 
352  histPlot->Draw();
353 
354  TFile fout(output_name.c_str(),"RECREATE");
355  fout.Close();
356  histPlot->WriteOutput(output_name);
357  histPlot->WriteHist(output_name);
358  treePlot->WriteOutput(output_name);
359 
360  g_Log << LogInfo << "Finished" << LogEnd;
361 }
362 
363 # ifndef __CINT__ // main function for stand-alone compilation
364 int main(){
365  example_H_to_ttbar_to_bWlnubWlnu();
366  return 0;
367 }
368 #endif
Generated on Fri Aug 12 2016 14:55:13 for RestFrames by   doxygen 1.8.9.1