LOGO
RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
example_top_to_bWlnu.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_top_to_bWlnu(const std::string& output_name =
42  "output_top_to_bWlnu.root"){
43 
44  // set particle masses and widths
45  double mtop = 173.21; // GeV, PDG 2016
46  double wtop = 1.41;
47  double mW = 80.385;
48  double wW = 2.085;
49 
50  // Number of events to generate
51  int Ngen = 100000;
52 
54  g_Log << LogInfo << "Initializing generator frames and tree..." << LogEnd;
56  ppLabGenFrame LAB_Gen("LAB_Gen","LAB");
57  ResonanceGenFrame T_Gen("T_Gen","t");
58  ResonanceGenFrame W_Gen("W_Gen","W");
59  VisibleGenFrame B_Gen("B_Gen","b");
60  VisibleGenFrame L_Gen("L_Gen","#it{l}");
61  InvisibleGenFrame NU_Gen("NU_Gen","#nu");
62 
63  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
64 
65  LAB_Gen.SetChildFrame(T_Gen);
66  T_Gen.AddChildFrame(B_Gen);
67  T_Gen.AddChildFrame(W_Gen);
68  W_Gen.AddChildFrame(L_Gen);
69  W_Gen.AddChildFrame(NU_Gen);
70 
71  if(LAB_Gen.InitializeTree())
72  g_Log << LogInfo << "...Successfully initialized generator tree" << LogEnd;
73  else
74  g_Log << LogError << "...Failed initializing generator tree" << LogEnd;
75 
76  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
77 
78  // set top mass and width
79  T_Gen.SetMass(mtop); T_Gen.SetWidth(wtop);
80  // set W mass and width
81  W_Gen.SetMass(mW); W_Gen.SetWidth(wW);
82 
83  // set b-jet/lepton pT and eta cuts
84  L_Gen.SetPtCut(20.); L_Gen.SetEtaCut(2.5);
85  B_Gen.SetPtCut(20.); B_Gen.SetEtaCut(2.5);
86 
87  if(LAB_Gen.InitializeAnalysis())
88  g_Log << LogInfo << "...Successfully initialized generator analysis" << std::endl << LogEnd;
89  else
90  g_Log << LogError << "...Failed initializing generator analysis" << LogEnd;
93 
95  g_Log << LogInfo << "Initializing reconstruction frames and trees..." << LogEnd;
97  LabRecoFrame LAB_Mt("LAB_Mt","LAB"); LabRecoFrame LAB_MW("LAB_MW","LAB");
98  DecayRecoFrame T_Mt("T_Mt","t"); DecayRecoFrame T_MW("T_MW","t");
99  DecayRecoFrame W_Mt("W_Mt","W"); DecayRecoFrame W_MW("W_MW","W");
100  VisibleRecoFrame B_Mt("B_Mt","b"); VisibleRecoFrame B_MW("B_MW","b");
101  VisibleRecoFrame L_Mt("L_Mt","#it{l}"); VisibleRecoFrame L_MW("L_MW","#it{l}");
102  InvisibleRecoFrame NU_Mt("NU_Mt","#nu"); InvisibleRecoFrame NU_MW("NU_MW","#nu");
103 
104  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
105 
106  LAB_Mt.SetChildFrame(T_Mt); LAB_MW.SetChildFrame(T_MW);
107  T_Mt.AddChildFrame(B_Mt); T_MW.AddChildFrame(B_MW);
108  T_Mt.AddChildFrame(W_Mt); T_MW.AddChildFrame(W_MW);
109  W_Mt.AddChildFrame(L_Mt); W_MW.AddChildFrame(L_MW);
110  W_Mt.AddChildFrame(NU_Mt); W_MW.AddChildFrame(NU_MW);
111 
112  if(LAB_Mt.InitializeTree() && LAB_MW.InitializeTree())
113  g_Log << LogInfo << "...Successfully initialized reconstruction trees" << LogEnd;
114  else
115  g_Log << LogError << "...Failed initializing reconstruction trees" << LogEnd;
116 
117  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
118 
119  // Invisible Groups
120  InvisibleGroup INV_Mt("INV_Mt","#nu Jigsaws");
121  INV_Mt.AddFrame(NU_Mt);
122  InvisibleGroup INV_MW("INV_MW","#nu Jigsaws");
123  INV_MW.AddFrame(NU_MW);
124 
125  // Set neutrino masses to zero
126  SetMassInvJigsaw NuM_Mt("NuM_Mt","M_{#nu} = 0");
127  INV_Mt.AddJigsaw(NuM_Mt);
128  SetMassInvJigsaw NuM_MW("NuM_MW","M_{#nu} = 0");
129  INV_MW.AddJigsaw(NuM_MW);
130 
131  // Set neutrino rapidity to that of visible particles
132  SetRapidityInvJigsaw NuR_Mt("NuR_Mt","#eta_{#nu} = #eta_{b+#it{l}}");
133  INV_Mt.AddJigsaw(NuR_Mt);
134  NuR_Mt.AddVisibleFrame(L_Mt);
135  NuR_Mt.AddVisibleFrame(B_Mt);
136  SetRapidityInvJigsaw NuR_MW("NuR_MW","#eta_{#nu} = #eta_{#it{l}}");
137  INV_MW.AddJigsaw(NuR_MW);
138  NuR_MW.AddVisibleFrame(L_MW);
139 
140  if(LAB_Mt.InitializeAnalysis() && LAB_MW.InitializeAnalysis())
141  g_Log << LogInfo << "...Successfully initialized analyses" << LogEnd;
142  else
143  g_Log << LogError << "...Failed initializing analyses" << LogEnd;
144 
147 
148  TreePlot* treePlot = new TreePlot("TreePlot","TreePlot");
149 
150  treePlot->SetTree(LAB_Gen);
151  treePlot->Draw("GenTree", "Generator Tree", true);
152 
153  treePlot->SetTree(LAB_Mt);
154  treePlot->Draw("RecoTree", "Reconstruction Tree");
155 
156  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
157 
158  // Declare observables for histogram booking
159  HistPlot* histPlot = new HistPlot("HistPlot","pp #rightarrow t #rightarrow W(#it{l} #nu) b");
160 
161  const HistPlotCategory& cat_Gen = histPlot->GetNewCategory("Gen", "Generator");
162  const HistPlotCategory& cat_minMt = histPlot->GetNewCategory("minMt", "min M_{t} Reco");
163  const HistPlotCategory& cat_minMW = histPlot->GetNewCategory("minMW", "min M_{W} Reco");
164 
165  const HistPlotVar& Mt = histPlot->GetNewVar("Mt", "M_{t}", 0., 280., "[GeV]");
166  const HistPlotVar& MW = histPlot->GetNewVar("MW", "M_{W}", 0., 150., "[GeV]");
167  const HistPlotVar& pTtoMt = histPlot->GetNewVar("pTtoMt","p_{T}^{top} / m_{top}", 0., 1.);
168  const HistPlotVar& cosT = histPlot->GetNewVar("cosT","cos #theta_{t}", -1., 1.);
169  const HistPlotVar& cosW = histPlot->GetNewVar("cosW","cos #theta_{W}", -1., 1.);
170  const HistPlotVar& dphiT = histPlot->GetNewVar("dphiT", "#Delta #phi_{t}", 0., 2.*acos(-1.));
171  const HistPlotVar& dphiW = histPlot->GetNewVar("dphiW", "#Delta #phi_{W}", 0., 2.*acos(-1.));
172  const HistPlotVar& DcosT = histPlot->GetNewVar("DcosT","#theta_{t} - #theta_{t}^{gen}", -1., 1.);
173  const HistPlotVar& DcosW = histPlot->GetNewVar("DcosW","#theta_{W} - #theta_{W}^{gen}", -1., 1.);
174  const HistPlotVar& DdphiT = histPlot->GetNewVar("DdphiT","#Delta #phi_{t} - #Delta #phi_{t}^{gen}", -1., 1.);
175  const HistPlotVar& DdphiW = histPlot->GetNewVar("DdphiW","#Delta #phi_{W} - #Delta #phi_{W}^{gen}", -1., 1.);
176 
177  histPlot->AddPlot(DcosT, cat_minMt+cat_minMW);
178  histPlot->AddPlot(DcosW, cat_minMt+cat_minMW);
179  histPlot->AddPlot(DdphiT, cat_minMt+cat_minMW);
180  histPlot->AddPlot(DdphiW, cat_minMt+cat_minMW);
181  histPlot->AddPlot(MW, cat_minMt+cat_minMW+cat_Gen);
182  histPlot->AddPlot(Mt, cat_minMt+cat_minMW+cat_Gen);
183  histPlot->AddPlot(Mt, MW, cat_minMt+cat_minMW+cat_Gen);
184  histPlot->AddPlot(Mt, pTtoMt, cat_minMt+cat_minMW);
185  histPlot->AddPlot(MW, pTtoMt, cat_minMt+cat_minMW);
186  histPlot->AddPlot(DcosW, MW, cat_minMt+cat_minMW);
187  histPlot->AddPlot(DcosT, Mt, cat_minMt+cat_minMW);
188  histPlot->AddPlot(Mt, pTtoMt, cat_minMt+cat_minMW);
189  histPlot->AddPlot(MW, pTtoMt, cat_minMt+cat_minMW);
190  histPlot->AddPlot(DcosT, DcosW, cat_minMt+cat_minMW);
191 
192 
195 
196  for(int igen = 0; igen < Ngen; igen++){
197  if(igen%((std::max(Ngen,10))/10) == 0)
198  g_Log << LogInfo << "Generating event " << igen << " of " << Ngen << LogEnd;
199 
200  // generate event
201  LAB_Gen.ClearEvent(); // clear the gen tree
202 
203  pTtoMt = gRandom->Rndm(); // give the Top some Pt
204  LAB_Gen.SetPToverM(pTtoMt);
205 
206  LAB_Gen.AnalyzeEvent(); // generate a new event
207 
208  TVector3 MET = LAB_Gen.GetInvisibleMomentum(); // Get the MET from gen tree
209  MET.SetZ(0.);
210 
211  // analyze event one way
212  LAB_Mt.ClearEvent(); // clear the reco tree
213 
214  L_Mt.SetLabFrameFourVector(L_Gen.GetFourVector()); // Set lepton 4-vec
215  B_Mt.SetLabFrameFourVector(B_Gen.GetFourVector()); // Set b-jet 4-vec
216  INV_Mt.SetLabFrameThreeVector(MET); // Set the MET in reco tree
217 
218  LAB_Mt.AnalyzeEvent(); //analyze the event
219 
220  // analyze event another way
221  LAB_MW.ClearEvent(); // clear the reco tree
222 
223  L_MW.SetLabFrameFourVector(L_Gen.GetFourVector()); // Set lepton 4-vec
224  B_MW.SetLabFrameFourVector(B_Gen.GetFourVector()); // Set b-jet 4-vec
225  INV_MW.SetLabFrameThreeVector(MET); // Set the MET in reco tree
226 
227  LAB_MW.AnalyzeEvent(); //analyze the event
228 
229  // Generator-level observables
230  double MTgen = T_Gen.GetMass();
231  double cosTgen = T_Gen.GetCosDecayAngle();
232  double dphiTgen = LAB_Gen.GetDeltaPhiDecayPlanes(T_Gen);
233  double MWgen = W_Gen.GetMass();
234  double cosWgen = W_Gen.GetCosDecayAngle();
235  double dphiWgen = T_Gen.GetDeltaPhiDecayPlanes(W_Gen);
236 
237  Mt = MTgen;
238  MW = MWgen;
239  histPlot->Fill(cat_Gen);
240 
241  // minMt observables
242  Mt = T_Mt.GetMass();
243  cosT = T_Mt.GetCosDecayAngle();
244  dphiT = LAB_Mt.GetDeltaPhiDecayPlanes(T_Mt);
245  MW = W_Mt.GetMass();
246  cosW = W_Mt.GetCosDecayAngle();
247  dphiW = T_Mt.GetDeltaPhiDecayPlanes(W_Mt);
248  DcosT = asin(sqrt(1.-cosT*cosT)*cosTgen-sqrt(1.-cosTgen*cosTgen)*cosT);
249  DdphiT = asin(sin(dphiT-dphiTgen));
250  DcosW = asin(sqrt(1.-cosW*cosW)*cosWgen-sqrt(1.-cosWgen*cosWgen)*cosW);
251  DdphiW = asin(sin(dphiW-dphiWgen));
252 
253  histPlot->Fill(cat_minMt);
254 
255  // minMW observables
256  Mt = T_MW.GetMass();
257  cosT = T_MW.GetCosDecayAngle();
258  dphiT = LAB_MW.GetDeltaPhiDecayPlanes(T_MW);
259  MW = W_MW.GetMass();
260  cosW = W_MW.GetCosDecayAngle();
261  dphiW = T_MW.GetDeltaPhiDecayPlanes(W_MW);
262  DcosT = asin(sqrt(1.-cosT*cosT)*cosTgen-sqrt(1.-cosTgen*cosTgen)*cosT);
263  DdphiT = asin(sin(dphiT-dphiTgen));
264  DcosW = asin(sqrt(1.-cosW*cosW)*cosWgen-sqrt(1.-cosWgen*cosWgen)*cosW);
265  DdphiW = asin(sin(dphiW-dphiWgen));
266 
267  histPlot->Fill(cat_minMW);
268  }
269 
270  histPlot->Draw();
271 
272  LAB_Gen.PrintGeneratorEfficiency();
273 
274  TFile fout(output_name.c_str(),"RECREATE");
275  fout.Close();
276  histPlot->WriteOutput(output_name);
277  histPlot->WriteHist(output_name);
278  treePlot->WriteOutput(output_name);
279 
280  g_Log << LogInfo << "Finished" << LogEnd;
281 
282 }
283 
284 # ifndef __CINT__ // main function for stand-alone compilation
285 int main(){
286  example_top_to_bWlnu();
287  return 0;
288 }
289 #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::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::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