LOGO
RestFrames  v1.0.1
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-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_H_to_ttbar_to_bWlnubWlnu(std::string output_name =
42  "output_H_to_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  std::vector<double> mH; // vary neutral Higgs mass
51  mH.push_back(400.);
52  mH.push_back(750.);
53  mH.push_back(1000.);
54  mH.push_back(1500.);
55 
56  // number of events to generate (per Higgs mass)
57  int Ngen = 10000;
58 
60  g_Log << LogInfo << "Initializing generator frames and tree..." << LogEnd;
62  ppLabGenFrame LAB_Gen("LAB_Gen","LAB");
63  DecayGenFrame H_Gen("H_Gen","H^{ 0}");
64  DecayGenFrame Ta_Gen("Ta_Gen","t_{a}");
65  DecayGenFrame Tb_Gen("Tb_Gen","t_{b}");
66  DecayGenFrame Wa_Gen("Wa_Gen","W_{a}");
67  DecayGenFrame Wb_Gen("Wb_Gen","W_{b}");
68  VisibleGenFrame Ba_Gen("Ba_Gen","b_{a}");
69  VisibleGenFrame La_Gen("La_Gen","#it{l}_{a}");
70  InvisibleGenFrame Na_Gen("Na_Gen","#nu_{a}");
71  VisibleGenFrame Bb_Gen("Bb_Gen","b_{b}");
72  VisibleGenFrame Lb_Gen("Lb_Gen","#it{l}_{b}");
73  InvisibleGenFrame Nb_Gen("Nb_Gen","#nu_{b}");
74 
75  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
76 
77  LAB_Gen.SetChildFrame(H_Gen);
78  H_Gen.AddChildFrame(Ta_Gen);
79  H_Gen.AddChildFrame(Tb_Gen);
80  Ta_Gen.AddChildFrame(Ba_Gen);
81  Ta_Gen.AddChildFrame(Wa_Gen);
82  Tb_Gen.AddChildFrame(Bb_Gen);
83  Tb_Gen.AddChildFrame(Wb_Gen);
84  Wa_Gen.AddChildFrame(La_Gen);
85  Wa_Gen.AddChildFrame(Na_Gen);
86  Wb_Gen.AddChildFrame(Lb_Gen);
87  Wb_Gen.AddChildFrame(Nb_Gen);
88 
89  if(LAB_Gen.InitializeTree())
90  g_Log << LogInfo << "...Successfully initialized generator tree" << LogEnd;
91  else
92  g_Log << LogError << "...Failed initializing generator tree" << LogEnd;
93 
94  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
95 
96  if(mH.size() < 1) return;
97 
98  // Higgs mass
99  H_Gen.SetMass(mH[0]);
100  // set top masses
101  Ta_Gen.SetMass(mT); Tb_Gen.SetMass(mT);
102  // set W masses
103  Wa_Gen.SetMass(mW); Wb_Gen.SetMass(mW);
104  // set B masses
105  Ba_Gen.SetMass(mB); Bb_Gen.SetMass(mB);
106  // set : masses
107  La_Gen.SetMass(mL); Lb_Gen.SetMass(mL);
108  // set neutrino masses
109  Na_Gen.SetMass(mN); Nb_Gen.SetMass(mN);
110 
111  // set b-jet/lepton pT/eta cuts
112  Ba_Gen.SetPtCut(20.); Bb_Gen.SetPtCut(20.);
113  Ba_Gen.SetEtaCut(2.5); Bb_Gen.SetEtaCut(2.5);
114  La_Gen.SetPtCut(15.); Lb_Gen.SetPtCut(15.);
115  La_Gen.SetEtaCut(2.5); Lb_Gen.SetEtaCut(2.5);
116 
117  if(LAB_Gen.InitializeAnalysis())
118  g_Log << LogInfo << "...Successfully initialized generator analysis" << LogEnd;
119  else
120  g_Log << LogError << "...Failed initializing generator analysis" << LogEnd;
124  g_Log << LogInfo << "Initializing reconstruction frames and tree..." << LogEnd;
126 
127  LabRecoFrame LAB("LAB","LAB");
128  DecayRecoFrame H("H","H^{ 0}");
129  DecayRecoFrame Ta("Ta","t_{a}");
130  DecayRecoFrame Tb("Tb","t_{b}");
131  DecayRecoFrame Wa("Wa","W_{a}");
132  DecayRecoFrame Wb("Wb","W_{b}");
133  VisibleRecoFrame Ba("Ba","b_{a}");
134  VisibleRecoFrame La("La","#it{l}_{a}");
135  InvisibleRecoFrame Na("Na","#nu_{a}");
136  VisibleRecoFrame Bb("Bb","b_{b}");
137  VisibleRecoFrame Lb("Lb","#it{l}_{b}");
138  InvisibleRecoFrame Nb("Nb","#nu_{b}");
139 
140  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
141 
142  LAB.SetChildFrame(H);
143  H.AddChildFrame(Ta);
144  H.AddChildFrame(Tb);
145  Ta.AddChildFrame(Ba);
146  Ta.AddChildFrame(Wa);
147  Tb.AddChildFrame(Bb);
148  Tb.AddChildFrame(Wb);
149  Wa.AddChildFrame(La);
150  Wa.AddChildFrame(Na);
151  Wb.AddChildFrame(Lb);
152  Wb.AddChildFrame(Nb);
153 
154  if(LAB.InitializeTree())
155  g_Log << LogInfo << "...Successfully initialized reconstruction tree" << LogEnd;
156  else
157  g_Log << LogError << "...Failed initializing reconstruction tree" << LogEnd;
158 
159  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
160 
162  std::string group_name;
163 
164  // Invisible Group
165  group_name = "#nu #nu Jigsaws";
166  InvisibleGroup INV("INV", group_name);
167  INV.AddFrame(Na);
168  INV.AddFrame(Nb);
169  CombinatoricGroup B("VIS","b-jet Jigsaws");
170  B.AddFrame(Ba);
171  B.AddFrame(Bb);
172  B.SetNElementsForFrame(Ba, 1);
173  B.SetNElementsForFrame(Bb, 1);
174 
176  std::string jigsaw_name;
177 
178  // Minimize difference Mt jigsaws
179  jigsaw_name = "M_{#nu #nu} ~ m_{#it{l} #it{l}}";
180  SetMassInvJigsaw NuNuM("NuNuM", jigsaw_name);
181  INV.AddJigsaw(NuNuM);
182 
183  jigsaw_name = "#eta_{#nu#nu} = #eta_{b #it{l} b #it{l}}";
184  SetRapidityInvJigsaw NuNuR("NuNuR", jigsaw_name);
185  INV.AddJigsaw(NuNuR);
186  NuNuR.AddVisibleFrames(H.GetListVisibleFrames());
187 
188  jigsaw_name = "min ( M_{top a}- M_{top b} )^{2}";
189  MinMassDiffInvJigsaw MinDeltaMt("MinDeltaMt", jigsaw_name, 2);
190  INV.AddJigsaw(MinDeltaMt);
191  MinDeltaMt.AddInvisibleFrame(Na, 0);
192  MinDeltaMt.AddInvisibleFrame(Nb, 1);
193  MinDeltaMt.AddVisibleFrames(La+Ba, 0);
194  MinDeltaMt.AddVisibleFrames(Lb+Bb, 1);
195  MinDeltaMt.AddMassFrame(La, 0);
196  MinDeltaMt.AddMassFrame(Lb, 1);
197 
198  // b-jet combinatoric jigsaws for all trees
199  jigsaw_name = "Minimize M(b #it{l} )_{a} , M(b #it{l} )_{b}";
200 
201  MinMassesCombJigsaw MinBL("MinBL", jigsaw_name);
202  B.AddJigsaw(MinBL);
203  MinBL.AddFrames(La+Ba,0);
204  MinBL.AddFrames(Lb+Bb,1);
205 
206  if(LAB.InitializeAnalysis())
207  g_Log << LogInfo << "...Successfully initialized analysis" << LogEnd;
208  else
209  g_Log << LogError << "...Failed initializing analysis" << LogEnd;
210 
213 
214  TreePlot* treePlot = new TreePlot("TreePlot","TreePlot");
215 
216  treePlot->SetTree(LAB_Gen);
217  treePlot->Draw("GenTree", "Generator Tree", true);
218 
219  treePlot->SetTree(LAB);
220  treePlot->Draw("RecoTree", "Reconstruction Tree");
221 
222  treePlot->SetTree(B);
223  treePlot->Draw("VisTree", "b-jet Jigsaws", true);
224 
225  treePlot->SetTree(INV);
226  treePlot->Draw("InvTree", "Inivisible Jigsaws", true);
227 
228  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
229 
230  std::string plot_title = "H^{ 0} #rightarrow t #bar{t} #rightarrow b W(#it{l} #nu) b W(#it{l} #nu)";
231  HistPlot* histPlot = new HistPlot("HistPlot", plot_title);
232 
233  RFList<const HistPlotCategory> cat_list;
234  int Nmass = mH.size();
235  for(int m = 0; m < Nmass; m++){
236  char smass[50], scat[50];
237  sprintf(scat, "MH%.0f", mH[m]);
238  sprintf(smass, "m_{H^{ 0}} = %.0f", mH[m]);
239  cat_list += histPlot->GetNewCategory(scat, smass);
240  }
241 
242  const HistPlotVar& MH = histPlot->GetNewVar("MH", "M_{H^{ 0}} / m_{H^{ 0}}^{ true}", 0., 2.);
243  const HistPlotVar& Eb_ta = histPlot->GetNewVar("Eb_ta", "E_{b a}^{top a} / E_{b a}^{top a gen}", 0., 2.);
244  const HistPlotVar& Eb_tb = histPlot->GetNewVar("Eb_tb", "E_{b b}^{top b} / E_{b b}^{top b gen}", 0., 2.);
245  const HistPlotVar& El_Wa = histPlot->GetNewVar("El_Wa", "E_{#it{l} a}^{W a} / E_{#it{l} a}^{W a true}", 0., 2.);
246  const HistPlotVar& El_Wb = histPlot->GetNewVar("El_Wb", "E_{#it{l} b}^{W b} / E_{#it{l} b}^{W b true}", 0., 2.);
247  const HistPlotVar& cosH = histPlot->GetNewVar("cosH","cos #theta_{H^{ 0}}", -1., 1.);
248  const HistPlotVar& costa = histPlot->GetNewVar("costa","cos #theta_{top a}", -1., 1.);
249  const HistPlotVar& costb = histPlot->GetNewVar("costb","cos #theta_{top b}", -1., 1.);
250  const HistPlotVar& cosWa = histPlot->GetNewVar("cosWa","cos #theta_{W a}", -1., 1.);
251  const HistPlotVar& cosWb = histPlot->GetNewVar("cosWb","cos #theta_{W b}", -1., 1.);
252  const HistPlotVar& DcosH = histPlot->GetNewVar("Dcostt","#theta_{H^{ 0}} - #theta_{H^{ 0}}^{true}",
253  -acos(-1.)/2., acos(-1.)/2.);
254  const HistPlotVar& Dcosta = histPlot->GetNewVar("Dcosta","#theta_{top a} - #theta_{top a}^{true}",
255  -acos(-1.)/2., acos(-1.)/2.);
256  const HistPlotVar& Dcostb = histPlot->GetNewVar("Dcostb","#theta_{top b} - #theta_{top b}^{true}",
257  -acos(-1.)/2., acos(-1.)/2.);
258  const HistPlotVar& DcosWa = histPlot->GetNewVar("DcosWa","#theta_{W a} - #theta_{W a}^{true}",
259  -acos(-1.)/2., acos(-1.)/2.);
260  const HistPlotVar& DcosWb = histPlot->GetNewVar("DcosWb","#theta_{W b} - #theta_{W b}^{true}",
261  -acos(-1.)/2., acos(-1.)/2.);
262 
263  histPlot->AddPlot(MH, cat_list);
264  histPlot->AddPlot(Eb_ta, cat_list);
265  histPlot->AddPlot(El_Wa, cat_list);
266  histPlot->AddPlot(DcosH, cat_list);
267  histPlot->AddPlot(Dcosta, cat_list);
268  histPlot->AddPlot(DcosWa, cat_list);
269 
270  histPlot->AddPlot(MH, Eb_ta, cat_list[1]);
271  histPlot->AddPlot(MH, El_Wa, cat_list[1]);
272  histPlot->AddPlot(Eb_ta, Eb_tb, cat_list[1]);
273  histPlot->AddPlot(El_Wa, El_Wb, cat_list[1]);
274  histPlot->AddPlot(Eb_ta, El_Wa, cat_list[1]);
275  histPlot->AddPlot(DcosH, MH, cat_list[1]);
276  histPlot->AddPlot(Dcosta, Eb_ta, cat_list[1]);
277  histPlot->AddPlot(DcosWa, El_Wa, cat_list[1]);
278  histPlot->AddPlot(DcosH, Dcosta, cat_list[1]);
279  histPlot->AddPlot(Dcosta, Dcostb, cat_list[1]);
280  histPlot->AddPlot(DcosWa, DcosWb, cat_list[1]);
281  histPlot->AddPlot(Dcosta, DcosWa, cat_list[1]);
282 
285 
286  for(int m = 0; m < Nmass; m++){
287  g_Log << LogInfo << "Generating events for H^{0} mass = " << mH[m] << LogEnd;
288 
289  H_Gen.SetMass(mH[m]);
290  LAB_Gen.InitializeAnalysis();
291 
292  for(int igen = 0; igen < Ngen; igen++){
293  if(igen%((std::max(Ngen,10))/10) == 0)
294  g_Log << LogInfo << "Generating event " << igen << " of " << Ngen << LogEnd;
295 
296  // generate event
297  LAB_Gen.ClearEvent(); // clear the gen tree
298 
299  LAB_Gen.AnalyzeEvent(); // generate a new event
300 
301  TVector3 MET = LAB_Gen.GetInvisibleMomentum(); // Get the MET from gen tree
302  MET.SetZ(0.);
303 
304  // analyze event
305  LAB.ClearEvent(); // clear the event
306 
307  INV.SetLabFrameThreeVector(MET); // Set MET in reco tree
308 
309  La.SetLabFrameFourVector(La_Gen.GetFourVector()); // Set lepton 4-vectors
310  Lb.SetLabFrameFourVector(Lb_Gen.GetFourVector());
311 
312  B.AddLabFrameFourVector(Ba_Gen.GetFourVector()); // Set b-jet 4-vectors
313  B.AddLabFrameFourVector(Bb_Gen.GetFourVector());
314 
315  LAB.AnalyzeEvent(); // analyze the event
316 
318  // Observable Calculations
320 
321  double MHgen = H_Gen.GetMass();
322  double Eb_tagen = Ba_Gen.GetFourVector(Ta_Gen).E();
323  double Eb_tbgen = Bb_Gen.GetFourVector(Tb_Gen).E();
324  double El_Wagen = La_Gen.GetFourVector(Wa_Gen).E();
325  double El_Wbgen = Lb_Gen.GetFourVector(Wb_Gen).E();
326  double cosHgen = H_Gen.GetCosDecayAngle();
327  double costagen = Ta_Gen.GetCosDecayAngle();
328  double costbgen = Tb_Gen.GetCosDecayAngle();
329  double cosWagen = Wa_Gen.GetCosDecayAngle();
330  double cosWbgen = Wb_Gen.GetCosDecayAngle();
331 
332  MH = H.GetMass()/MHgen;
333  Eb_ta = Ba.GetFourVector(Ta).E()/Eb_tagen;
334  Eb_tb = Bb.GetFourVector(Tb).E()/Eb_tbgen;
335  El_Wa = La.GetFourVector(Wa).E()/El_Wagen;
336  El_Wb = Lb.GetFourVector(Wb).E()/El_Wbgen;
337  cosH = H.GetCosDecayAngle();
338  costa = Ta.GetCosDecayAngle();
339  costb = Tb.GetCosDecayAngle();
340  cosWa = Wa.GetCosDecayAngle();
341  cosWb = Wb.GetCosDecayAngle();
342  DcosH = asin(sqrt(1.-cosH*cosH)*cosHgen-sqrt(1.-cosHgen*cosHgen)*cosH);
343  Dcosta = asin(sqrt(1.-costa*costa)*costagen-sqrt(1.-costagen*costagen)*costa);
344  Dcostb = asin(sqrt(1.-costb*costb)*costbgen-sqrt(1.-costbgen*costbgen)*costb);
345  DcosWa = asin(sqrt(1.-cosWa*cosWa)*cosWagen-sqrt(1.-cosWagen*cosWagen)*cosWa);
346  DcosWb = asin(sqrt(1.-cosWb*cosWb)*cosWbgen-sqrt(1.-cosWbgen*cosWbgen)*cosWb);
347 
348  histPlot->Fill(cat_list[m]);
349  }
350 
351  LAB_Gen.PrintGeneratorEfficiency();
352  }
353 
354  histPlot->Draw();
355 
356  TFile fout(output_name.c_str(),"RECREATE");
357  fout.Close();
358  histPlot->WriteOutput(output_name);
359  histPlot->WriteHist(output_name);
360  treePlot->WriteOutput(output_name);
361 
362  g_Log << LogInfo << "Finished" << LogEnd;
363 }
364 
365 # ifndef __CINT__ // main function for stand-alone compilation
366 int main(){
367  example_H_to_ttbar_to_bWlnubWlnu();
368  return 0;
369 }
370 #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::ppLabGenFrame
Definition: ppLabGenFrame.hh:40
RestFrames::DecayGenFrame
Definition: DecayGenFrame.hh:43
RestFrames::RFList
Definition: RFList.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::VisibleRecoFrame
Definition: VisibleRecoFrame.hh:41
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