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