LOGO
RestFrames  v1.0.0
RestFrames HEP Event Analysis Software Library
example_H_to_WlnuWlnu.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_WlnuWlnu(const std::string& output_name =
40  "output_H_to_WlnuWlnu.root"){
41 
42  // set particle masses and widths
43  double mW = 80.385; // GeV, PDG 2016
44  double wW = 2.085;
45  double mL = 0.106; // muons
46  double mN = 0.;
47 
48  std::vector<double> mH; // vary neutral Higgs mass
49  mH.push_back(125.);
50  mH.push_back(400.);
51  mH.push_back(750.);
52  mH.push_back(1000.);
53  mH.push_back(1500.);
54 
55  // Number of events to generate (per H mass)
56  int Ngen = 10000;
57 
59  g_Log << LogInfo << "Initializing generator frames and tree..." << LogEnd;
61  ppLabGenFrame LAB_Gen("LAB_Gen","LAB");
62  DecayGenFrame H_Gen("H_Gen","H^{0}");
63  ResonanceGenFrame Wa_Gen("Wa_Gen","W_{a}");
64  ResonanceGenFrame Wb_Gen("Wb_Gen","W_{b}");
65  VisibleGenFrame La_Gen("La_Gen","#it{l}_{a}");
66  InvisibleGenFrame Na_Gen("Na_Gen","#nu_{a}");
67  VisibleGenFrame Lb_Gen("Lb_Gen","#it{l}_{b}");
68  InvisibleGenFrame Nb_Gen("Nb_Gen","#nu_{b}");
69 
70  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
71 
72  LAB_Gen.SetChildFrame(H_Gen);
73  H_Gen.AddChildFrame(Wa_Gen);
74  H_Gen.AddChildFrame(Wb_Gen);
75  Wa_Gen.AddChildFrame(La_Gen);
76  Wa_Gen.AddChildFrame(Na_Gen);
77  Wb_Gen.AddChildFrame(Lb_Gen);
78  Wb_Gen.AddChildFrame(Nb_Gen);
79 
80  if(LAB_Gen.InitializeTree())
81  g_Log << LogInfo << "...Successfully initialized generator tree" << LogEnd;
82  else
83  g_Log << LogError << "...Failed initializing generator tree" << LogEnd;
84  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
85 
86  if(mH.size() < 1) return;
87 
88  // set Higgs masses
89  H_Gen.SetMass(mH[0]);
90  // set W masses and widths
91  Wa_Gen.SetMass(mW); Wa_Gen.SetWidth(wW);
92  Wb_Gen.SetMass(mW); Wb_Gen.SetWidth(wW);
93  // set lepton and neutrino masses
94  La_Gen.SetMass(mL); Lb_Gen.SetMass(mL);
95 
96  // set lepton pT and eta cuts
97  La_Gen.SetPtCut(10.); Lb_Gen.SetPtCut(10.);
98  La_Gen.SetEtaCut(2.5); Lb_Gen.SetEtaCut(2.5);
99 
100  if(LAB_Gen.InitializeAnalysis())
101  g_Log << LogInfo << "...Successfully initialized generator analysis" << LogEnd;
102  else
103  g_Log << LogError << "...Failed initializing generator analysis" << LogEnd;
106 
108  g_Log << LogInfo << "Initializing reconstruction frames and trees..." << LogEnd;
110  LabRecoFrame LAB("LAB","LAB");
111  DecayRecoFrame H("H","H^{ 0}");
112  DecayRecoFrame Wa("Wa","W_{a}");
113  DecayRecoFrame Wb("Wb","W_{b}");
114  VisibleRecoFrame La("La","#it{l}_{a}");
115  InvisibleRecoFrame Na("Na","#nu_{a}");
116  VisibleRecoFrame Lb("Lb","#it{l}_{b}");
117  InvisibleRecoFrame Nb("Nb","#nu_{b}");
118 
119  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
120 
121  LAB.SetChildFrame(H);
122  H.AddChildFrame(Wa);
123  H.AddChildFrame(Wb);
124  Wa.AddChildFrame(La);
125  Wa.AddChildFrame(Na);
126  Wb.AddChildFrame(Lb);
127  Wb.AddChildFrame(Nb);
128 
129  if(LAB.InitializeTree())
130  g_Log << LogInfo << "...Successfully initialized reconstruction tree" << LogEnd;
131  else
132  g_Log << LogError << "...Failed initializing reconstruction tree" << LogEnd;
133 
134  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
135 
136  // Invisible Group
137  InvisibleGroup INV("INV","#nu #nu Jigsaws");
138  INV.AddFrame(Na);
139  INV.AddFrame(Nb);
140 
141  // Set nu nu mass equal to l l mass
142  SetMassInvJigsaw NuNuM("NuNuM", "M_{#nu#nu} = m_{#it{l}#it{l}}");
143  INV.AddJigsaw(NuNuM);
144 
145  SetRapidityInvJigsaw NuNuR("NuNuR", "#eta_{#nu#nu} = #eta_{#it{l}#it{l}}");
146  INV.AddJigsaw(NuNuR);
147  NuNuR.AddVisibleFrames(LAB.GetListVisibleFrames());
148 
149  //MinMassesSqInvJigsaw MinMW("MinMW","min M_{W}, M_{Wa}= M_{Wb}",2);
150  ContraBoostInvJigsaw MinMW("MinMW","min M_{W}, M_{Wa}= M_{Wb}");
151  INV.AddJigsaw(MinMW);
152  MinMW.AddVisibleFrame(La, 0);
153  MinMW.AddVisibleFrame(Lb, 1);
154  MinMW.AddInvisibleFrame(Na, 0);
155  MinMW.AddInvisibleFrame(Nb, 1);
156 
157  if(LAB.InitializeAnalysis())
158  g_Log << LogInfo << "...Successfully initialized analysis" << LogEnd;
159  else
160  g_Log << LogError << "...Failed initializing analysis" << LogEnd;
161 
164 
165  TreePlot* treePlot = new TreePlot("TreePlot","TreePlot");
166 
167  treePlot->SetTree(LAB_Gen);
168  treePlot->Draw("GenTree", "Generator Tree", true);
169 
170  treePlot->SetTree(LAB);
171  treePlot->Draw("RecoTree", "Reconstruction Tree");
172 
173  treePlot->SetTree(INV);
174  treePlot->Draw("InvTree", "Invisible Jigsaws", true);
175 
176  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
177 
178  std::string plot_title = "H^{ 0} #rightarrow W(#it{l} #nu) W(#it{l} #nu)";
179  HistPlot* histPlot = new HistPlot("HistPlot", plot_title);
180 
181  RFList<const HistPlotCategory> cat_list;
182  int Nmass = mH.size();
183  for(int m = 0; m < Nmass; m++){
184  char smass[50], scat[50];
185  sprintf(scat, "MH%.0f", mH[m]);
186  sprintf(smass, "m_{H^{ 0}} = %.0f", mH[m]);
187  cat_list += histPlot->GetNewCategory(scat, smass);
188  }
189 
190  const HistPlotCategory& cat_Gen = histPlot->GetNewCategory("Gen", "Generator");
191 
192  const HistPlotVar& MH = histPlot->GetNewVar("MH", "M_{H^{ 0}}", 0., 3000., "[GeV]");
193  const HistPlotVar& MHN = histPlot->GetNewVar("MHN", "M_{H^{ 0}} / m_{H^{ 0}}^{ true}", 0., 2.);
194  const HistPlotVar& MWaN = histPlot->GetNewVar("MWaN", "M_{Wa} / m_{Wa}^{ true}", 0., 3.);
195  const HistPlotVar& cosH = histPlot->GetNewVar("cosH","cos #theta_{H^{ 0}}", -1., 1.);
196  const HistPlotVar& cosWa = histPlot->GetNewVar("cosWa","cos #theta_{W_{a}}", -1., 1.);
197  const HistPlotVar& dphiH = histPlot->GetNewVar("dphiH", "#Delta #phi_{H^{ 0}}", 0., 2.*acos(-1.));
198  const HistPlotVar& DcosH = histPlot->GetNewVar("DcosH","#theta_{H^{ 0}} - #theta_{H^{ 0}}^{true}",
199  -acos(-1.)/2., acos(-1.)/2.);
200  const HistPlotVar& DcosWa = histPlot->GetNewVar("DcosWa","#theta_{W_{a}} - #theta_{W_{a}}^{true}",
201  -acos(-1.)/2., acos(-1.)/2.);
202 
203  histPlot->AddPlot(MH, cat_list);
204  histPlot->AddPlot(MHN, cat_list);
205  histPlot->AddPlot(MWaN, cat_list);
206  histPlot->AddPlot(DcosH, cat_list);
207  histPlot->AddPlot(DcosWa, cat_list);
208  histPlot->AddPlot(MHN, MWaN, cat_list[1]);
209  histPlot->AddPlot(MHN, DcosH, cat_list[0]);
210  histPlot->AddPlot(MWaN, DcosWa, cat_list[1]);
211 
214 
215  for(int m = 0; m < Nmass; m++){
216  g_Log << LogInfo << "Generating events for H^{0} mass = " << mH[m] << LogEnd;
217 
218  H_Gen.SetMass(mH[m]);
219  LAB_Gen.InitializeAnalysis();
220 
221  for(int igen = 0; igen < Ngen; igen++){
222  if(igen%((std::max(Ngen,10))/10) == 0)
223  g_Log << LogInfo << "Generating event " << igen << " of " << Ngen << LogEnd;
224 
225  // generate event
226  LAB_Gen.ClearEvent(); // clear the gen tree
227 
228  LAB_Gen.AnalyzeEvent(); // generate a new event
229 
230  // analyze event
231  LAB.ClearEvent(); // clear the reco tree
232 
233  TVector3 MET = LAB_Gen.GetInvisibleMomentum(); // Get the MET from gen tree
234  MET.SetZ(0.);
235  INV.SetLabFrameThreeVector(MET); // Set the MET in reco tree
236 
237  La.SetLabFrameFourVector(La_Gen.GetFourVector()); // set lepton 4-vectors
238  Lb.SetLabFrameFourVector(Lb_Gen.GetFourVector());
239 
240  LAB.AnalyzeEvent(); // analyze the event
241 
242  // Generator-level observables
243  double cosHgen = H_Gen.GetCosDecayAngle();
244  cosH = cosHgen;
245  double cosWagen = Wa_Gen.GetCosDecayAngle();
246  cosWa = cosWagen;
247 
248  if(m == 0)
249  histPlot->Fill(cat_Gen);
250 
251  // Reconstruction-level observables
252  MH = H.GetMass();
253  MHN = H.GetMass()/H_Gen.GetMass();
254  MWaN = Wa.GetMass()/Wa_Gen.GetMass();
255  cosH = H.GetCosDecayAngle();
256  cosWa = Wa.GetCosDecayAngle();
257  DcosH = asin(sqrt(1.-cosH*cosH)*cosHgen-sqrt(1.-cosHgen*cosHgen)*cosH);
258  DcosWa = asin(sqrt(1.-cosWa*cosWa)*cosWagen-sqrt(1.-cosWagen*cosWagen)*cosWa);
259 
260  histPlot->Fill(cat_list[m]);
261  }
262 
263  LAB_Gen.PrintGeneratorEfficiency();
264  }
265 
266  histPlot->Draw();
267 
268  TFile fout(output_name.c_str(),"RECREATE");
269  fout.Close();
270  histPlot->WriteOutput(output_name);
271  histPlot->WriteHist(output_name);
272  treePlot->WriteOutput(output_name);
273 
274  g_Log << LogInfo << "Finished" << LogEnd;
275 }
276 
277 # ifndef __CINT__ // main function for stand-alone compilation
278 int main(){
279  example_H_to_WlnuWlnu();
280  return 0;
281 }
282 #endif
Generated on Fri Aug 12 2016 14:55:13 for RestFrames by   doxygen 1.8.9.1