LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
example_Wlnu.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_Wlnu(const std::string& output_name = "output_Wlnu.root"){
40 
41  double mW = 80.385; // GeV, PDG 2016
42  double wW = 2.085;
43 
44  // Number of events to generate
45  int Ngen = 100000;
46 
48  g_Log << LogInfo << "Initializing generator frames and tree..." << LogEnd;
50  LabGenFrame LAB_Gen("LAB_Gen","LAB");
51  ResonanceGenFrame W_Gen("W_Gen","W");
52  VisibleGenFrame L_Gen("L_Gen","#it{l}");
53  InvisibleGenFrame NU_Gen("NU_Gen","#nu");
54 
55  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
56 
57  LAB_Gen.SetChildFrame(W_Gen);
58  W_Gen.AddChildFrame(L_Gen);
59  W_Gen.AddChildFrame(NU_Gen);
60 
61  if(LAB_Gen.InitializeTree())
62  g_Log << LogInfo << "...Successfully initialized generator tree" << LogEnd;
63  else
64  g_Log << LogError << "...Failed initializing generator tree" << LogEnd;
65 
66  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
67 
68  // set W pole mass and width
69  W_Gen.SetMass(mW); W_Gen.SetWidth(wW);
70 
71  // set lepton pT and eta cuts
72  L_Gen.SetPtCut(20.); L_Gen.SetEtaCut(2.5);
73 
74  if(LAB_Gen.InitializeAnalysis())
75  g_Log << LogInfo << "...Successfully initialized generator analysis" << std::endl << LogEnd;
76  else
77  g_Log << LogError << "...Failed initializing generator analysis" << LogEnd;
80 
82  g_Log << LogInfo << "Initializing reconstruction frames and trees..." << LogEnd;
84  LabRecoFrame LAB("LAB","LAB");
85  DecayRecoFrame W("W","W");
86  VisibleRecoFrame L("L","#it{l}");
87  InvisibleRecoFrame NU("NU","#nu");
88 
89  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
90 
91  LAB.SetChildFrame(W);
92  W.AddChildFrame(L);
93  W.AddChildFrame(NU);
94 
95  if(LAB.InitializeTree())
96  g_Log << LogInfo << "...Successfully initialized reconstruction trees" << LogEnd;
97  else
98  g_Log << LogError << "...Failed initializing reconstruction trees" << LogEnd;
99 
100  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
101 
102  // Now we add invisible jigsaws
103  InvisibleGroup INV("INV","Neutrino Jigsaws");
104  INV.AddFrame(NU);
105 
106  // Set the neutrino mass
107  SetMassInvJigsaw MassJigsaw("MassJigsaw","m_{#nu} = 0");
108  INV.AddJigsaw(MassJigsaw);
109 
110  // Set the neutrino rapidity
111  SetRapidityInvJigsaw RapidityJigsaw("RapidityJigsaw","#eta_{#nu} = #eta_{#it{l}}");
112  INV.AddJigsaw(RapidityJigsaw);
113  RapidityJigsaw.AddVisibleFrame(L);
114 
115  if(LAB.InitializeAnalysis())
116  g_Log << LogInfo << "...Successfully initialized analyses" << LogEnd;
117  else
118  g_Log << LogError << "...Failed initializing analyses" << LogEnd;
119 
122 
123  TreePlot* tree_plot = new TreePlot("TreePlot","TreePlot");
124 
125  // generator tree
126  tree_plot->SetTree(LAB_Gen);
127  tree_plot->Draw("GenTree", "Generator Tree", true);
128 
129  // reconstruction tree
130  tree_plot->SetTree(LAB);
131  tree_plot->Draw("RecoTree", "Reconstruction Tree");
132 
133  // Invisible Jigsaws tree
134  tree_plot->SetTree(INV);
135  tree_plot->Draw("InvTree", "InvisibleJigsaws", true);
136 
137  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
138 
139  // Declare observables for histogram booking
140  HistPlot* hist_plot = new HistPlot("HistPlot","W #rightarrow #it{l} #nu");
141 
142  const HistPlotVar& MW = hist_plot->GetNewVar("MW", "M_{W}", 30., 120., "[GeV]");
143  const HistPlotVar& cosW = hist_plot->GetNewVar("cosW","cos #phi_{W}", -1., 1.);
144  const HistPlotVar& dphiW = hist_plot->GetNewVar("dphiW", "#Delta #phi_{W}", 0., 2.*acos(-1.));
145  const HistPlotVar& DcosW = hist_plot->GetNewVar("DcosW","#phi_{W} - #phi_{W}^{true}", -0.5, 0.5);
146  const HistPlotVar& DdphiW = hist_plot->GetNewVar("DdphiW","#Delta #phi_{W} - #Delta #phi_{W}^{true}", -0.5, 0.5);
147  const HistPlotVar& pTWoMW = hist_plot->GetNewVar("pTW","p_{T}^{W} / m_{W}",0.,1.);
148 
149  const HistPlotCategory& cat_Gen = hist_plot->GetNewCategory("Reco", "Generator");
150  const HistPlotCategory& cat_Reco = hist_plot->GetNewCategory("Reco", "Reconstruction");
151 
152  hist_plot->AddPlot(DcosW, cat_Reco);
153  hist_plot->AddPlot(DdphiW, cat_Reco);
154  hist_plot->AddPlot(MW, cat_Gen+cat_Reco);
155  hist_plot->AddPlot(DcosW, MW, cat_Reco);
156  hist_plot->AddPlot(DdphiW, MW, cat_Reco);
157  hist_plot->AddPlot(MW, pTWoMW, cat_Reco);
158  hist_plot->AddPlot(DcosW, pTWoMW, cat_Reco);
159  hist_plot->AddPlot(DdphiW, pTWoMW, cat_Reco);
160 
161  for(int igen = 0; igen < Ngen; igen++){
162  if(igen%((std::max(Ngen,10))/10) == 0)
163  g_Log << LogInfo << "Generating event " << igen << " of " << Ngen << LogEnd;
164 
165  // generate event
166  LAB_Gen.ClearEvent(); // clear the gen tree
167 
168  pTWoMW = LAB_Gen.GetRandom();
169  LAB_Gen.SetPToverM(pTWoMW); // give the W some Pt
170  double PzW = mW*(2.*LAB_Gen.GetRandom()-1.);
171  LAB_Gen.SetLongitudinalMomentum(PzW); // give the W some Pz
172 
173  LAB_Gen.AnalyzeEvent(); // generate a new event
174 
175  // analyze event
176  LAB.ClearEvent(); // clear the reco tree
177 
178  L.SetLabFrameFourVector(L_Gen.GetFourVector()); // Set lepton 4-vec
179 
180  TVector3 MET = LAB_Gen.GetInvisibleMomentum(); // Get the MET from gen tree
181  MET.SetZ(0.);
182  INV.SetLabFrameThreeVector(MET); // Set the MET in reco tree
183 
184  LAB.AnalyzeEvent(); // analyze the event
185 
186  // generator-level observables
187  MW = W_Gen.GetMass();
188  cosW = cos(W_Gen.GetDeltaPhiDecayAngle());
189  dphiW = LAB_Gen.GetDeltaPhiDecayPlanes(W_Gen);
190 
191  hist_plot->Fill(cat_Gen);
192 
193  // calculate observables
194  MW = W.GetMass();
195  cosW = cos(W.GetDeltaPhiDecayAngle());
196  dphiW = LAB.GetDeltaPhiDecayPlanes(W);
197 
198  double cosWgen = cos(W_Gen.GetDeltaPhiDecayAngle());
199  double dphiWgen = LAB_Gen.GetDeltaPhiDecayPlanes(W_Gen);
200  DcosW = asin(sqrt(1.-cosW*cosW)*cosWgen-sqrt(1.-cosWgen*cosWgen)*cosW);
201  DdphiW = asin(sin(dphiW-dphiWgen));
202 
203  hist_plot->Fill(cat_Reco);
204  }
205 
206  hist_plot->Draw();
207 
208  TFile fout(output_name.c_str(),"RECREATE");
209  fout.Close();
210  hist_plot->WriteOutput(output_name);
211  hist_plot->WriteHist(output_name);
212  tree_plot->WriteOutput(output_name);
213 
214 }
215 
216 # ifndef __CINT__ // main function for stand-alone compilation
217 int main(){
218  example_Wlnu();
219  return 0;
220 }
221 #endif