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