LOGO
RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
example_ttbar_to_bWlnubWlnu.C
Go to the documentation of this file.
1 // RestFrames: particle physics event analysis library
3 // --------------------------------------------------------------------
4 // Copyright (c) 2014-2019, 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 #if (!defined(__CINT__) && !defined(__CLING__))
31 #define COMPILER
32 #endif
33 #if defined(__MAKECINT__) || defined(__ROOTCLING__) || defined(COMPILER)
34 #include "RestFrames/RestFrames.hh"
35 #else
36 RestFrames::RFKey ensure_autoload(1);
37 #endif
38 
39 using namespace RestFrames;
40 
41 void example_ttbar_to_bWlnubWlnu(const std::string output_name =
42  "output_ttbar_to_bWlnubWlnu.root"){
43 
44  double mT = 173.21; // GeV, PDG 2016
45  double mW = 80.385;
46  double mB = 4.18;
47  double mL = 0.106;
48  double mN = 0.;
49 
50  // number of events to generate
51  int Ngen = 100000;
52 
54  g_Log << LogInfo << "Initializing generator frames and tree..." << LogEnd;
56  ppLabGenFrame LAB_Gen("LAB_Gen","LAB");
57  DecayGenFrame TT_Gen("TT_Gen","t #bar{t}");
58  DecayGenFrame Ta_Gen("Ta_Gen","t_{a}");
59  DecayGenFrame Tb_Gen("Tb_Gen","t_{b}");
60  DecayGenFrame Wa_Gen("Wa_Gen","W_{a}");
61  DecayGenFrame Wb_Gen("Wb_Gen","W_{b}");
62  VisibleGenFrame Ba_Gen("Ba_Gen","b_{a}");
63  VisibleGenFrame La_Gen("La_Gen","#it{l}_{a}");
64  InvisibleGenFrame Na_Gen("Na_Gen","#nu_{a}");
65  VisibleGenFrame Bb_Gen("Bb_Gen","b_{b}");
66  VisibleGenFrame Lb_Gen("Lb_Gen","#it{l}_{b}");
67  InvisibleGenFrame Nb_Gen("Nb_Gen","#nu_{b}");
68 
69  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
70 
71  LAB_Gen.SetChildFrame(TT_Gen);
72  TT_Gen.AddChildFrame(Ta_Gen);
73  TT_Gen.AddChildFrame(Tb_Gen);
74  Ta_Gen.AddChildFrame(Ba_Gen);
75  Ta_Gen.AddChildFrame(Wa_Gen);
76  Tb_Gen.AddChildFrame(Bb_Gen);
77  Tb_Gen.AddChildFrame(Wb_Gen);
78  Wa_Gen.AddChildFrame(La_Gen);
79  Wa_Gen.AddChildFrame(Na_Gen);
80  Wb_Gen.AddChildFrame(Lb_Gen);
81  Wb_Gen.AddChildFrame(Nb_Gen);
82 
83  if(LAB_Gen.InitializeTree())
84  g_Log << LogInfo << "...Successfully initialized generator tree" << LogEnd;
85  else
86  g_Log << LogError << "...Failed initializing generator tree" << LogEnd;
87 
88  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
89 
90  // non-resonant ttbar production
91  TT_Gen.SetVariableMass();
92  // set top masses
93  Ta_Gen.SetMass(mT); Tb_Gen.SetMass(mT);
94  // set W masses
95  Wa_Gen.SetMass(mW); Wb_Gen.SetMass(mW);
96  // set B masses
97  Ba_Gen.SetMass(mB); Bb_Gen.SetMass(mB);
98  // set : masses
99  La_Gen.SetMass(mL); Lb_Gen.SetMass(mL);
100  // set neutrino masses
101  Na_Gen.SetMass(mN); Nb_Gen.SetMass(mN);
102 
103  // set b-jet/lepton pT/eta cuts
104  Ba_Gen.SetPtCut(20.); Bb_Gen.SetPtCut(20.);
105  Ba_Gen.SetEtaCut(2.5); Bb_Gen.SetEtaCut(2.5);
106  La_Gen.SetPtCut(15.); Lb_Gen.SetPtCut(15.);
107  La_Gen.SetEtaCut(2.5); Lb_Gen.SetEtaCut(2.5);
108 
109  if(LAB_Gen.InitializeAnalysis())
110  g_Log << LogInfo << "...Successfully initialized generator analysis" << LogEnd;
111  else
112  g_Log << LogError << "...Failed initializing generator analysis" << LogEnd;
116  g_Log << LogInfo << "Initializing reconstruction frames and trees..." << LogEnd;
118  LabRecoFrame LAB_R1("LAB_R1","LAB"); LabRecoFrame LAB_R2("LAB_R2","LAB");
119  DecayRecoFrame TT_R1("TT_R1","t #bar{t}"); DecayRecoFrame TT_R2("TT_R2","t #bar{t}");
120  DecayRecoFrame Ta_R1("Ta_R1","t_{a}"); DecayRecoFrame Ta_R2("Ta_R2","t_{a}");
121  DecayRecoFrame Tb_R1("Tb_R1","t_{b}"); DecayRecoFrame Tb_R2("Tb_R2","t_{b}");
122  DecayRecoFrame Wa_R1("Wa_R1","W_{a}"); DecayRecoFrame Wa_R2("Wa_R2","W_{a}");
123  DecayRecoFrame Wb_R1("Wb_R1","W_{b}"); DecayRecoFrame Wb_R2("Wb_R2","W_{b}");
124  VisibleRecoFrame Ba_R1("Ba_R1","b_{a}"); VisibleRecoFrame Ba_R2("Ba_R2","b_{a}");
125  VisibleRecoFrame La_R1("La_R1","#it{l}_{a}"); VisibleRecoFrame La_R2("La_R2","#it{l}_{a}");
126  InvisibleRecoFrame Na_R1("Na_R1","#nu_{a}"); InvisibleRecoFrame Na_R2("Na_R2","#nu_{a}");
127  VisibleRecoFrame Bb_R1("Bb_R1","b_{b}"); VisibleRecoFrame Bb_R2("Bb_R2","b_{b}");
128  VisibleRecoFrame Lb_R1("Lb_R1","#it{l}_{b}"); VisibleRecoFrame Lb_R2("Lb_R2","#it{l}_{b}");
129  InvisibleRecoFrame Nb_R1("Nb_R1","#nu_{b}"); InvisibleRecoFrame Nb_R2("Nb_R2","#nu_{b}");
130 
131  LabRecoFrame LAB_R3("LAB_R3","LAB"); LabRecoFrame LAB_R4("LAB_R4","LAB");
132  DecayRecoFrame TT_R3("TT_R3","t #bar{t}"); DecayRecoFrame TT_R4("TT_R4","t #bar{t}");
133  DecayRecoFrame Ta_R3("Ta_R3","t_{a}"); DecayRecoFrame Ta_R4("Ta_R4","t_{a}");
134  DecayRecoFrame Tb_R3("Tb_R3","t_{b}"); DecayRecoFrame Tb_R4("Tb_R4","t_{b}");
135  DecayRecoFrame Wa_R3("Wa_R3","W_{a}"); DecayRecoFrame Wa_R4("Wa_R4","W_{a}");
136  DecayRecoFrame Wb_R3("Wb_R3","W_{b}"); DecayRecoFrame Wb_R4("Wb_R4","W_{b}");
137  VisibleRecoFrame Ba_R3("Ba_R3","b_{a}"); VisibleRecoFrame Ba_R4("Ba_R4","b_{a}");
138  VisibleRecoFrame La_R3("La_R3","#it{l}_{a}"); VisibleRecoFrame La_R4("La_R4","#it{l}_{a}");
139  InvisibleRecoFrame Na_R3("Na_R3","#nu_{a}"); InvisibleRecoFrame Na_R4("Na_R4","#nu_{a}");
140  VisibleRecoFrame Bb_R3("Bb_R3","b_{b}"); VisibleRecoFrame Bb_R4("Bb_R4","b_{b}");
141  VisibleRecoFrame Lb_R3("Lb_R3","#it{l}_{b}"); VisibleRecoFrame Lb_R4("Lb_R4","#it{l}_{b}");
142  InvisibleRecoFrame Nb_R3("Nb_R3","#nu_{b}"); InvisibleRecoFrame Nb_R4("Nb_R4","#nu_{b}");
143 
144  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
145 
146  LAB_R1.SetChildFrame(TT_R1); LAB_R2.SetChildFrame(TT_R2);
147  TT_R1.AddChildFrame(Ta_R1); TT_R2.AddChildFrame(Ta_R2);
148  TT_R1.AddChildFrame(Tb_R1); TT_R2.AddChildFrame(Tb_R2);
149  Ta_R1.AddChildFrame(Ba_R1); Ta_R2.AddChildFrame(Ba_R2);
150  Ta_R1.AddChildFrame(Wa_R1); Ta_R2.AddChildFrame(Wa_R2);
151  Tb_R1.AddChildFrame(Bb_R1); Tb_R2.AddChildFrame(Bb_R2);
152  Tb_R1.AddChildFrame(Wb_R1); Tb_R2.AddChildFrame(Wb_R2);
153  Wa_R1.AddChildFrame(La_R1); Wa_R2.AddChildFrame(La_R2);
154  Wa_R1.AddChildFrame(Na_R1); Wa_R2.AddChildFrame(Na_R2);
155  Wb_R1.AddChildFrame(Lb_R1); Wb_R2.AddChildFrame(Lb_R2);
156  Wb_R1.AddChildFrame(Nb_R1); Wb_R2.AddChildFrame(Nb_R2);
157 
158  LAB_R3.SetChildFrame(TT_R3); LAB_R4.SetChildFrame(TT_R4);
159  TT_R3.AddChildFrame(Ta_R3); TT_R4.AddChildFrame(Ta_R4);
160  TT_R3.AddChildFrame(Tb_R3); TT_R4.AddChildFrame(Tb_R4);
161  Ta_R3.AddChildFrame(Ba_R3); Ta_R4.AddChildFrame(Ba_R4);
162  Ta_R3.AddChildFrame(Wa_R3); Ta_R4.AddChildFrame(Wa_R4);
163  Tb_R3.AddChildFrame(Bb_R3); Tb_R4.AddChildFrame(Bb_R4);
164  Tb_R3.AddChildFrame(Wb_R3); Tb_R4.AddChildFrame(Wb_R4);
165  Wa_R3.AddChildFrame(La_R3); Wa_R4.AddChildFrame(La_R4);
166  Wa_R3.AddChildFrame(Na_R3); Wa_R4.AddChildFrame(Na_R4);
167  Wb_R3.AddChildFrame(Lb_R3); Wb_R4.AddChildFrame(Lb_R4);
168  Wb_R3.AddChildFrame(Nb_R3); Wb_R4.AddChildFrame(Nb_R4);
169 
170  if(LAB_R1.InitializeTree() && LAB_R2.InitializeTree() &&
171  LAB_R3.InitializeTree() && LAB_R4.InitializeTree())
172  g_Log << LogInfo << "...Successfully initialized reconstruction trees" << LogEnd;
173  else
174  g_Log << LogError << "...Failed initializing reconstruction trees" << LogEnd;
175 
176  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
177 
179  std::string group_name;
180 
181  // Invisible Group
182  group_name = "#splitline{#nu #nu Jigsaws for}{min M_{top} , M_{top}^{ a} = M_{top}^{ b}}";
183  InvisibleGroup INV_R1("INV_R1", group_name);
184  INV_R1.AddFrame(Na_R1);
185  INV_R1.AddFrame(Nb_R1);
186  // Combinatoric Group for b's
187  CombinatoricGroup B_R1("VIS_R1","b-jet Jigsaws");
188  B_R1.AddFrame(Ba_R1);
189  B_R1.AddFrame(Bb_R1);
190  // b-jet frames must have at least one element
191  B_R1.SetNElementsForFrame(Ba_R1, 1);
192  B_R1.SetNElementsForFrame(Bb_R1, 1);
193 
194  group_name = "#splitline{#nu #nu Jigsaws for}{min M_{W}, M_{W}^{ a} = M_{W}^{ b}}";
195  InvisibleGroup INV_R2("INV_R2", group_name);
196  INV_R2.AddFrame(Na_R2);
197  INV_R2.AddFrame(Nb_R2);
198  CombinatoricGroup B_R2("VIS_R2","b-jet Jigsaws");
199  B_R2.AddFrame(Ba_R2);
200  B_R2.AddFrame(Bb_R2);
201  B_R2.SetNElementsForFrame(Ba_R2, 1);
202  B_R2.SetNElementsForFrame(Bb_R2, 1);
203 
204  group_name = "#splitline{#nu #nu Jigsaws for}{min M_{top a}^{2}+ M_{top b}^{2}}";
205  InvisibleGroup INV_R3("INV_R3", group_name);
206  INV_R3.AddFrame(Na_R3);
207  INV_R3.AddFrame(Nb_R3);
208  CombinatoricGroup B_R3("VIS_R3","b-jet Jigsaws");
209  B_R3.AddFrame(Ba_R3);
210  B_R3.AddFrame(Bb_R3);
211  B_R3.SetNElementsForFrame(Ba_R3, 1);
212  B_R3.SetNElementsForFrame(Bb_R3, 1);
213 
214  group_name = "#splitline{#nu #nu Jigsaws for}{min (M_{top a}- M_{top b})^{2}}";
215  InvisibleGroup INV_R4("INV_R4", group_name);
216  INV_R4.AddFrame(Na_R4);
217  INV_R4.AddFrame(Nb_R4);
218  CombinatoricGroup B_R4("VIS_R4","b-jet Jigsaws");
219  B_R4.AddFrame(Ba_R4);
220  B_R4.AddFrame(Bb_R4);
221  B_R4.SetNElementsForFrame(Ba_R4, 1);
222  B_R4.SetNElementsForFrame(Bb_R4, 1);
223 
225  std::string jigsaw_name;
226 
227  // Minimize equal top masses neutrino jigsaws
228  jigsaw_name = "M_{#nu#nu} = f(m_{b#it{l}b#it{l}} , m_{b#it{l}}^{ a} , m_{b#it{l}}^{ b})";
229  SetMassInvJigsaw NuNuM_R1("NuNuM_R1", jigsaw_name);
230  INV_R1.AddJigsaw(NuNuM_R1);
231 
232  jigsaw_name = "#eta_{#nu#nu} = #eta_{b #it{l} b #it{l}}";
233  SetRapidityInvJigsaw NuNuR_R1("NuNuR_R1", jigsaw_name);
234  INV_R1.AddJigsaw(NuNuR_R1);
235  NuNuR_R1.AddVisibleFrames(La_R1+Ba_R1+Lb_R1+Bb_R1);
236 
237  jigsaw_name = "min M_{top}, M_{top}^{ a} = M_{top}^{ b}";
238  ContraBoostInvJigsaw MinMt_R1("MinMt_R1", jigsaw_name);
239  INV_R1.AddJigsaw(MinMt_R1);
240  MinMt_R1.AddVisibleFrames(La_R1+Ba_R1, 0);
241  MinMt_R1.AddVisibleFrames(Lb_R1+Bb_R1, 1);
242  MinMt_R1.AddInvisibleFrame(Na_R1, 0);
243  MinMt_R1.AddInvisibleFrame(Nb_R1, 1);
244 
245  // Minimize equal W masses neutrino jigsaws
246  jigsaw_name = "M_{#nu#nu} = f(m_{#it{l}#it{l}} , m_{#it{l}}^{ a} , m_{#it{l}}^{ b})";
247  SetMassInvJigsaw NuNuM_R2("NuNuM_R2", jigsaw_name);
248  INV_R2.AddJigsaw(NuNuM_R2);
249 
250  jigsaw_name = "#eta_{#nu#nu} = #eta_{b #it{l} b #it{l}}";
251  SetRapidityInvJigsaw NuNuR_R2("NuNuR_R2", jigsaw_name);
252  INV_R2.AddJigsaw(NuNuR_R2);
253  NuNuR_R2.AddVisibleFrames(La_R2+Ba_R2+Lb_R2+Bb_R2);
254 
255  jigsaw_name = "min M_{W}, M_{W}^{ a} = M_{W}^{ b}";
256  ContraBoostInvJigsaw MinMW_R2("MinMW_R2", jigsaw_name);
257  INV_R2.AddJigsaw(MinMW_R2);
258  MinMW_R2.AddVisibleFrame(La_R2, 0);
259  MinMW_R2.AddVisibleFrame(Lb_R2, 1);
260  MinMW_R2.AddInvisibleFrame(Na_R2, 0);
261  MinMW_R2.AddInvisibleFrame(Nb_R2, 1);
262 
263  // Minimize sum Mt^2 jigsaws
264  jigsaw_name = "M_{#nu#nu} = f(m_{#it{l}#it{l}} , m_{#it{l}}^{ a} , m_{#it{l}}^{ b})";
265  SetMassInvJigsaw NuNuM_R3("NuNuM_R3", jigsaw_name);
266  INV_R3.AddJigsaw(NuNuM_R3);
267 
268  jigsaw_name = "#eta_{#nu#nu} = #eta_{b #it{l} b #it{l}}";
269  SetRapidityInvJigsaw NuNuR_R3("NuNuR_R3", jigsaw_name);
270  INV_R3.AddJigsaw(NuNuR_R3);
271  NuNuR_R3.AddVisibleFrames(LAB_R3.GetListVisibleFrames());
272 
273  jigsaw_name = "min #Sigma M_{top}^{2}";
274  MinMassesSqInvJigsaw MinMt_R3("MinMt_R3", jigsaw_name, 2);
275  INV_R3.AddJigsaw(MinMt_R3);
276  MinMt_R3.AddInvisibleFrame(Na_R3, 0);
277  MinMt_R3.AddInvisibleFrame(Nb_R3, 1);
278  MinMt_R3.AddVisibleFrames(La_R3+Ba_R3, 0);
279  MinMt_R3.AddVisibleFrames(Lb_R3+Bb_R3, 1);
280  MinMt_R3.AddMassFrame(La_R3, 0);
281  MinMt_R3.AddMassFrame(Lb_R3, 1);
282 
283  // Minimize difference Mt jigsaws
284  jigsaw_name = "M_{#nu#nu} = f(m_{#it{l}#it{l}} , m_{#it{l}}^{ a} , m_{#it{l}}^{ b})";
285  SetMassInvJigsaw NuNuM_R4("NuNuM_R4", jigsaw_name);
286  INV_R4.AddJigsaw(NuNuM_R4);
287 
288  jigsaw_name = "#eta_{#nu#nu} = #eta_{b #it{l} b #it{l}}";
289  SetRapidityInvJigsaw NuNuR_R4("NuNuR_R4", jigsaw_name);
290  INV_R4.AddJigsaw(NuNuR_R4);
291  NuNuR_R4.AddVisibleFrames(LAB_R4.GetListVisibleFrames());
292 
293  jigsaw_name = "min ( M_{top a}- M_{top b} )^{2}";
294  MinMassDiffInvJigsaw MinDeltaMt_R4("MinDeltaMt_R4", jigsaw_name, 2);
295  INV_R4.AddJigsaw(MinDeltaMt_R4);
296  MinDeltaMt_R4.AddInvisibleFrame(Na_R4, 0);
297  MinDeltaMt_R4.AddInvisibleFrame(Nb_R4, 1);
298  MinDeltaMt_R4.AddVisibleFrames(La_R4+Ba_R4, 0);
299  MinDeltaMt_R4.AddVisibleFrames(Lb_R4+Bb_R4, 1);
300  MinDeltaMt_R4.AddMassFrame(La_R4, 0);
301  MinDeltaMt_R4.AddMassFrame(Lb_R4, 1);
302 
303  // b-jet combinatoric jigsaws for all trees
304  jigsaw_name = "Minimize M(b #it{l} )_{a} , M(b #it{l} )_{b}";
305 
306  MinMassesCombJigsaw MinBL_R1("MinBL_R1", jigsaw_name);
307  B_R1.AddJigsaw(MinBL_R1);
308  MinBL_R1.AddFrames(La_R1+Ba_R1,0);
309  MinBL_R1.AddFrames(Lb_R1+Bb_R1,1);
310 
311  MinMassesCombJigsaw MinBL_R2("MinBL_R2", jigsaw_name);
312  B_R2.AddJigsaw(MinBL_R2);
313  MinBL_R2.AddFrames(La_R2+Ba_R2,0);
314  MinBL_R2.AddFrames(Lb_R2+Bb_R2,1);
315 
316  MinMassesCombJigsaw MinBL_R3("MinBL_R3", jigsaw_name);
317  B_R3.AddJigsaw(MinBL_R3);
318  MinBL_R3.AddFrames(La_R3+Ba_R3,0);
319  MinBL_R3.AddFrames(Lb_R3+Bb_R3,1);
320 
321  MinMassesCombJigsaw MinBL_R4("MinBL_R4", jigsaw_name);
322  B_R4.AddJigsaw(MinBL_R4);
323  MinBL_R4.AddFrames(La_R4+Ba_R4,0);
324  MinBL_R4.AddFrames(Lb_R4+Bb_R4,1);
325 
326  if(LAB_R1.InitializeAnalysis() && LAB_R2.InitializeAnalysis() &&
327  LAB_R3.InitializeAnalysis() && LAB_R4.InitializeAnalysis())
328  g_Log << LogInfo << "...Successfully initialized analysis" << LogEnd;
329  else
330  g_Log << LogError << "...Failed initializing analysis" << LogEnd;
331 
334 
335  TreePlot* treePlot = new TreePlot("TreePlot","TreePlot");
336 
337  treePlot->SetTree(LAB_Gen);
338  treePlot->Draw("GenTree", "Generator Tree", true);
339 
340  treePlot->SetTree(LAB_R1);
341  treePlot->Draw("RecoTree", "Reconstruction Tree");
342 
343  treePlot->SetTree(B_R1);
344  treePlot->Draw("VisTree", "b-jet Jigsaws", true);
345 
346  treePlot->SetTree(INV_R1);
347  treePlot->Draw("InvTree_R1", "Inivisibl Jigsaws", true);
348 
349  treePlot->SetTree(INV_R2);
350  treePlot->Draw("InvTree_R2", "Inivisibl Jigsaws", true);
351 
352  treePlot->SetTree(INV_R3);
353  treePlot->Draw("InvTree_R3", "Inivisibl Jigsaws", true);
354 
355  treePlot->SetTree(INV_R4);
356  treePlot->Draw("InvTree_R4", "Inivisibl Jigsaws", true);
357 
358  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
359 
360  HistPlot* histPlot = new HistPlot("HistPlot", "t #bar{t} #rightarrow b W(#it{l} #nu) b W(#it{l} #nu)");
361 
362  const HistPlotCategory& cat_Gen = histPlot->GetNewCategory("Gen", "Generator");
363  const HistPlotCategory& cat_R1 = histPlot->GetNewCategory("Reco1", "M_{top}^{ a} = M_{top}^{ b} Reco");
364  const HistPlotCategory& cat_R2 = histPlot->GetNewCategory("Reco2", "M_{W}^{ a} = M_{W}^{ b} Reco");
365  const HistPlotCategory& cat_R3 = histPlot->GetNewCategory("Reco3", "min #Sigma M_{top}^{ 2} Reco");
366  const HistPlotCategory& cat_R4 = histPlot->GetNewCategory("Reco4", "min #Delta M_{top} Reco");
367 
368  const HistPlotVar& Mtt = histPlot->GetNewVar("Mtt", "M_{t #bar{t}} / m_{t #bar{t}}", 0., 2.);
369  const HistPlotVar& Eb_ta = histPlot->GetNewVar("Eb_ta", "E_{b a}^{top a} / E_{b a}^{top a gen}", 0., 2.);
370  const HistPlotVar& Eb_tb = histPlot->GetNewVar("Eb_tb", "E_{b b}^{top b} / E_{b b}^{top b gen}", 0., 2.);
371  const HistPlotVar& El_Wa = histPlot->GetNewVar("El_Wa", "E_{#it{l} a}^{W a} / E_{#it{l} a}^{W a gen}", 0., 2.);
372  const HistPlotVar& El_Wb = histPlot->GetNewVar("El_Wb", "E_{#it{l} b}^{W b} / E_{#it{l} b}^{W b gen}", 0., 2.);
373  const HistPlotVar& costt = histPlot->GetNewVar("costt","cos #theta_{t #bar{t}}", -1., 1.);
374  const HistPlotVar& costa = histPlot->GetNewVar("costa","cos #theta_{top a}", -1., 1.);
375  const HistPlotVar& costb = histPlot->GetNewVar("costb","cos #theta_{top b}", -1., 1.);
376  const HistPlotVar& cosWa = histPlot->GetNewVar("cosWa","cos #theta_{W a}", -1., 1.);
377  const HistPlotVar& cosWb = histPlot->GetNewVar("cosWb","cos #theta_{W b}", -1., 1.);
378  const HistPlotVar& Dcostt = histPlot->GetNewVar("Dcostt","#theta_{t #bar{t}} - #theta_{t #bar{t}}^{gen}",
379  -acos(-1.)/2., acos(-1.)/2.);
380  const HistPlotVar& Dcosta = histPlot->GetNewVar("Dcosta","#theta_{top a} - #theta_{top a}^{gen}",
381  -acos(-1.)/2., acos(-1.)/2.);
382  const HistPlotVar& Dcostb = histPlot->GetNewVar("Dcostb","#theta_{top b} - #theta_{top b}^{gen}",
383  -acos(-1.)/2., acos(-1.)/2.);
384  const HistPlotVar& DcosWa = histPlot->GetNewVar("DcosWa","#theta_{W a} - #theta_{W a}^{gen}",
385  -acos(-1.)/2., acos(-1.)/2.);
386  const HistPlotVar& DcosWb = histPlot->GetNewVar("DcosWb","#theta_{W b} - #theta_{W b}^{gen}",
387  -acos(-1.)/2., acos(-1.)/2.);
388 
389  histPlot->AddPlot(Mtt, cat_R1+cat_R2+cat_R3+cat_R4);
390  histPlot->AddPlot(Eb_ta, cat_R1+cat_R2+cat_R3+cat_R4);
391  histPlot->AddPlot(El_Wa, cat_R1+cat_R2+cat_R3+cat_R4);
392  histPlot->AddPlot(Dcostt, cat_R1+cat_R2+cat_R3+cat_R4);
393  histPlot->AddPlot(Dcosta, cat_R1+cat_R2+cat_R3+cat_R4);
394  histPlot->AddPlot(DcosWa, cat_R1+cat_R2+cat_R3+cat_R4);
395 
396  histPlot->AddPlot(Mtt, Eb_ta, cat_R4);
397  histPlot->AddPlot(Mtt, El_Wa, cat_R4);
398  histPlot->AddPlot(Eb_ta, Eb_tb, cat_R4);
399  histPlot->AddPlot(El_Wa, El_Wb, cat_R4);
400  histPlot->AddPlot(Eb_ta, El_Wa, cat_R4);
401  histPlot->AddPlot(Eb_ta, El_Wb, cat_R4);
402  histPlot->AddPlot(Dcostt, Mtt, cat_R4);
403  histPlot->AddPlot(Dcosta, Eb_ta, cat_R4);
404  histPlot->AddPlot(DcosWa, El_Wa, cat_R4);
405  histPlot->AddPlot(Dcostt, Dcosta, cat_R4);
406  histPlot->AddPlot(Dcosta, Dcostb, cat_R4);
407  histPlot->AddPlot(DcosWa, DcosWb, cat_R4);
408  histPlot->AddPlot(Dcosta, DcosWa, cat_R4);
409 
412 
413  for(int igen = 0; igen < Ngen; igen++){
414  if(igen%((std::max(Ngen,10))/10) == 0)
415  g_Log << LogInfo << "Generating event " << igen << " of " << Ngen << LogEnd;
416 
417  // generate event
418  LAB_Gen.ClearEvent(); // clear the gen tree
419 
420  LAB_Gen.AnalyzeEvent(); // generate a new event
421 
422  // analyze event three different ways
423  TVector3 MET = LAB_Gen.GetInvisibleMomentum(); // Get the MET from gen tree
424  MET.SetZ(0.);
425 
426  LAB_R1.ClearEvent();
427  LAB_R2.ClearEvent();
428  LAB_R3.ClearEvent();
429  LAB_R4.ClearEvent();
430 
431  INV_R1.SetLabFrameThreeVector(MET);
432  INV_R2.SetLabFrameThreeVector(MET);
433  INV_R3.SetLabFrameThreeVector(MET);
434  INV_R4.SetLabFrameThreeVector(MET);
435 
436  La_R1.SetLabFrameFourVector(La_Gen.GetFourVector());
437  Lb_R1.SetLabFrameFourVector(Lb_Gen.GetFourVector());
438  La_R2.SetLabFrameFourVector(La_Gen.GetFourVector());
439  Lb_R2.SetLabFrameFourVector(Lb_Gen.GetFourVector());
440  La_R3.SetLabFrameFourVector(La_Gen.GetFourVector());
441  Lb_R3.SetLabFrameFourVector(Lb_Gen.GetFourVector());
442  La_R4.SetLabFrameFourVector(La_Gen.GetFourVector());
443  Lb_R4.SetLabFrameFourVector(Lb_Gen.GetFourVector());
444 
445  std::vector<RFKey> B_R1_ID; // ID for tracking jets in tree
446  B_R1_ID.push_back(B_R1.AddLabFrameFourVector(Ba_Gen.GetFourVector()));
447  B_R1_ID.push_back(B_R1.AddLabFrameFourVector(Bb_Gen.GetFourVector()));
448  B_R2.AddLabFrameFourVector(Ba_Gen.GetFourVector());
449  B_R2.AddLabFrameFourVector(Bb_Gen.GetFourVector());
450  B_R3.AddLabFrameFourVector(Ba_Gen.GetFourVector());
451  B_R3.AddLabFrameFourVector(Bb_Gen.GetFourVector());
452  B_R4.AddLabFrameFourVector(Ba_Gen.GetFourVector());
453  B_R4.AddLabFrameFourVector(Bb_Gen.GetFourVector());
454 
455  LAB_R1.AnalyzeEvent(); // analyze the event
456  LAB_R2.AnalyzeEvent();
457  LAB_R3.AnalyzeEvent();
458  LAB_R4.AnalyzeEvent();
459 
461  // Observable Calculations
463 
464  double Mttgen = TT_Gen.GetMass();
465  double Eb_tagen = Ba_Gen.GetFourVector(Ta_Gen).E();
466  double Eb_tbgen = Bb_Gen.GetFourVector(Tb_Gen).E();
467  double El_Wagen = La_Gen.GetFourVector(Wa_Gen).E();
468  double El_Wbgen = Lb_Gen.GetFourVector(Wb_Gen).E();
469  double costtgen = TT_Gen.GetCosDecayAngle();
470  double costagen = Ta_Gen.GetCosDecayAngle();
471  double costbgen = Tb_Gen.GetCosDecayAngle();
472  double cosWagen = Wa_Gen.GetCosDecayAngle();
473  double cosWbgen = Wb_Gen.GetCosDecayAngle();
474 
475  Mtt = TT_R1.GetMass()/Mttgen;
476  // Mta = Ta_R1.GetMass();
477  // Mtb = Tb_R1.GetMass();
478  // MWa = Wa_R1.GetMass();
479  // MWb = Wb_R1.GetMass();
480  Eb_ta = Ba_R1.GetFourVector(Ta_R1).E()/Eb_tagen;
481  Eb_tb = Bb_R1.GetFourVector(Tb_R1).E()/Eb_tbgen;
482  El_Wa = La_R1.GetFourVector(Wa_R1).E()/El_Wagen;
483  El_Wb = Lb_R1.GetFourVector(Wb_R1).E()/El_Wbgen;
484  costt = TT_R1.GetCosDecayAngle();
485  costa = Ta_R1.GetCosDecayAngle();
486  costb = Tb_R1.GetCosDecayAngle();
487  cosWa = Wa_R1.GetCosDecayAngle();
488  cosWb = Wb_R1.GetCosDecayAngle();
489  Dcostt = asin(sqrt(1.-costt*costt)*costtgen-sqrt(1.-costtgen*costtgen)*costt);
490  Dcosta = asin(sqrt(1.-costa*costa)*costagen-sqrt(1.-costagen*costagen)*costa);
491  Dcostb = asin(sqrt(1.-costb*costb)*costbgen-sqrt(1.-costbgen*costbgen)*costb);
492  DcosWa = asin(sqrt(1.-cosWa*cosWa)*cosWagen-sqrt(1.-cosWagen*cosWagen)*cosWa);
493  DcosWb = asin(sqrt(1.-cosWb*cosWb)*cosWbgen-sqrt(1.-cosWbgen*cosWbgen)*cosWb);
494 
495  histPlot->Fill(cat_R1);
496 
497  Mtt = TT_R2.GetMass()/Mttgen;
498  // Mta = Ta_R2.GetMass();
499  // Mtb = Tb_R2.GetMass();
500  // MWa = Wa_R2.GetMass();
501  // MWb = Wb_R2.GetMass();
502  Eb_ta = Ba_R2.GetFourVector(Ta_R2).E()/Eb_tagen;
503  Eb_tb = Bb_R2.GetFourVector(Tb_R2).E()/Eb_tbgen;
504  El_Wa = La_R2.GetFourVector(Wa_R2).E()/El_Wagen;
505  El_Wb = Lb_R2.GetFourVector(Wb_R2).E()/El_Wbgen;
506  costt = TT_R2.GetCosDecayAngle();
507  costa = Ta_R2.GetCosDecayAngle();
508  costb = Tb_R2.GetCosDecayAngle();
509  cosWa = Wa_R2.GetCosDecayAngle();
510  cosWb = Wb_R2.GetCosDecayAngle();
511  Dcostt = asin(sqrt(1.-costt*costt)*costtgen-sqrt(1.-costtgen*costtgen)*costt);
512  Dcosta = asin(sqrt(1.-costa*costa)*costagen-sqrt(1.-costagen*costagen)*costa);
513  Dcostb = asin(sqrt(1.-costb*costb)*costbgen-sqrt(1.-costbgen*costbgen)*costb);
514  DcosWa = asin(sqrt(1.-cosWa*cosWa)*cosWagen-sqrt(1.-cosWagen*cosWagen)*cosWa);
515  DcosWb = asin(sqrt(1.-cosWb*cosWb)*cosWbgen-sqrt(1.-cosWbgen*cosWbgen)*cosWb);
516 
517  histPlot->Fill(cat_R2);
518 
519  Mtt = TT_R3.GetMass()/Mttgen;
520  // Mta = Ta_R3.GetMass();
521  // Mtb = Tb_R3.GetMass();
522  // MWa = Wa_R3.GetMass();
523  // MWb = Wb_R3.GetMass();
524  Eb_ta = Ba_R3.GetFourVector(Ta_R3).E()/Eb_tagen;
525  Eb_tb = Bb_R3.GetFourVector(Tb_R3).E()/Eb_tbgen;
526  El_Wa = La_R3.GetFourVector(Wa_R3).E()/El_Wagen;
527  El_Wb = Lb_R3.GetFourVector(Wb_R3).E()/El_Wbgen;
528  costt = TT_R3.GetCosDecayAngle();
529  costa = Ta_R3.GetCosDecayAngle();
530  costb = Tb_R3.GetCosDecayAngle();
531  cosWa = Wa_R3.GetCosDecayAngle();
532  cosWb = Wb_R3.GetCosDecayAngle();
533  Dcostt = asin(sqrt(1.-costt*costt)*costtgen-sqrt(1.-costtgen*costtgen)*costt);
534  Dcosta = asin(sqrt(1.-costa*costa)*costagen-sqrt(1.-costagen*costagen)*costa);
535  Dcostb = asin(sqrt(1.-costb*costb)*costbgen-sqrt(1.-costbgen*costbgen)*costb);
536  DcosWa = asin(sqrt(1.-cosWa*cosWa)*cosWagen-sqrt(1.-cosWagen*cosWagen)*cosWa);
537  DcosWb = asin(sqrt(1.-cosWb*cosWb)*cosWbgen-sqrt(1.-cosWbgen*cosWbgen)*cosWb);
538 
539  histPlot->Fill(cat_R3);
540 
541  Mtt = TT_R4.GetMass()/Mttgen;
542  // Mta = Ta_R4.GetMass();
543  // Mtb = Tb_R4.GetMass();
544  // MWa = Wa_R4.GetMass();
545  // MWb = Wb_R4.GetMass();
546  Eb_ta = Ba_R4.GetFourVector(Ta_R4).E()/Eb_tagen;
547  Eb_tb = Bb_R4.GetFourVector(Tb_R4).E()/Eb_tbgen;
548  El_Wa = La_R4.GetFourVector(Wa_R4).E()/El_Wagen;
549  El_Wb = Lb_R4.GetFourVector(Wb_R4).E()/El_Wbgen;
550  costt = TT_R4.GetCosDecayAngle();
551  costa = Ta_R4.GetCosDecayAngle();
552  costb = Tb_R4.GetCosDecayAngle();
553  cosWa = Wa_R4.GetCosDecayAngle();
554  cosWb = Wb_R4.GetCosDecayAngle();
555  Dcostt = asin(sqrt(1.-costt*costt)*costtgen-sqrt(1.-costtgen*costtgen)*costt);
556  Dcosta = asin(sqrt(1.-costa*costa)*costagen-sqrt(1.-costagen*costagen)*costa);
557  Dcostb = asin(sqrt(1.-costb*costb)*costbgen-sqrt(1.-costbgen*costbgen)*costb);
558  DcosWa = asin(sqrt(1.-cosWa*cosWa)*cosWagen-sqrt(1.-cosWagen*cosWagen)*cosWa);
559  DcosWb = asin(sqrt(1.-cosWb*cosWb)*cosWbgen-sqrt(1.-cosWbgen*cosWbgen)*cosWb);
560 
561  histPlot->Fill(cat_R4);
562 
563  }
564 
565  histPlot->Draw();
566 
567  LAB_Gen.PrintGeneratorEfficiency();
568 
569  TFile fout(output_name.c_str(),"RECREATE");
570  fout.Close();
571  histPlot->WriteOutput(output_name);
572  histPlot->WriteHist(output_name);
573  treePlot->WriteOutput(output_name);
574 
575  g_Log << LogInfo << "Finished" << LogEnd;
576 }
577 
578 # ifndef __CINT__ // main function for stand-alone compilation
579 int main(){
580  example_ttbar_to_bWlnubWlnu();
581  return 0;
582 }
583 #endif
RestFrames::TreePlot::Draw
void Draw(const std::string &name="", const std::string &title="", bool invert_colors=false, bool invert_node_colors=false)
Draws tree.
Definition: TreePlot.cc:127
RestFrames::MinMassDiffInvJigsaw
Definition: MinMassDiffInvJigsaw.hh:40
RestFrames::HistPlot::WriteHist
void WriteHist(const std::string &filename)
Stores all histograms in root file.
Definition: HistPlot.cc:581
RestFrames::InvisibleGenFrame
Definition: InvisibleGenFrame.hh:41
RestFrames::HistPlotCategory
Definition: HistPlotCategory.hh:40
RestFrames::ppLabGenFrame
Definition: ppLabGenFrame.hh:40
RestFrames::DecayGenFrame
Definition: DecayGenFrame.hh:43
RestFrames::HistPlot
Definition: HistPlot.hh:44
RestFrames::RFPlot::WriteOutput
void WriteOutput(const std::string &filename)
Stores all plots in root file.
Definition: RFPlot.cc:59
RestFrames.hh
RestFrames::ContraBoostInvJigsaw
Definition: ContraBoostInvJigsaw.hh:37
RestFrames::VisibleRecoFrame
Definition: VisibleRecoFrame.hh:41
RestFrames::MinMassesSqInvJigsaw
Definition: MinMassesSqInvJigsaw.hh:37
RestFrames::SetRapidityInvJigsaw
Definition: SetRapidityInvJigsaw.hh:37
RestFrames::LabRecoFrame
Definition: LabRecoFrame.hh:44
RestFrames::MinMassesCombJigsaw
Definition: MinMassesCombJigsaw.hh:40
RestFrames::HistPlot::GetNewVar
HistPlotVar const & GetNewVar(const std::string &name, const std::string &title, double minval, double maxval, const std::string &unit="")
Defines new variable for plotting.
Definition: HistPlot.cc:86
RestFrames::VisibleGenFrame
Definition: VisibleGenFrame.hh:41
RestFrames::HistPlot::Fill
void Fill(double weight=1.)
Fills an entry for every variable according to each's current value.
Definition: HistPlot.cc:242
RestFrames::InvisibleRecoFrame
Definition: InvisibleRecoFrame.hh:41
RestFrames::TreePlot::SetTree
void SetTree(const RestFrame &frame)
Sets a tree that shows how a frame is connected to its children.
Definition: TreePlot.cc:160
RestFrames::TreePlot
Definition: TreePlot.hh:50
RestFrames::HistPlot::GetNewCategory
HistPlotCategory const & GetNewCategory(const std::string &name, const std::string &title)
Defines new category for plotting.
Definition: HistPlot.cc:95
RestFrames::SetMassInvJigsaw
Definition: SetMassInvJigsaw.hh:37
RestFrames::CombinatoricGroup
Definition: CombinatoricGroup.hh:39
RestFrames::RFKey
Definition: RFKey.hh:38
RestFrames::HistPlotVar
Definition: HistPlotVar.hh:40
RestFrames::HistPlot::AddPlot
void AddPlot(const HistPlotVar &var, RestFrames::RFList< const HistPlotCategory > cats=RestFrames::RFList< const RestFrames::HistPlotCategory >(), bool invert_colors=false)
Defines new 1D plot.
Definition: HistPlot.cc:102
RestFrames::DecayRecoFrame
Definition: DecayRecoFrame.hh:43
RestFrames::InvisibleGroup
Definition: InvisibleGroup.hh:40
RestFrames::HistPlot::Draw
void Draw(bool invert_colors=false)
Draws every histogram.
Definition: HistPlot.cc:272
Generated on Thu Dec 17 2020 14:04:44 for RestFrames by   doxygen 1.8.17