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