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