LOGO
RestFrames  v1.0.0
RestFrames HEP Event Analysis Software Library
example_Hp_to_HggWlnu.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_Hp_to_HggWlnu(const std::string& output_name =
40  "output_Hp_to_HggWlnu.root"){
41 
42  // set particle masses and widths
43  double mH = 125.;
44  double mW = 80.385; // GeV, PDG 2016
45  double wW = 2.085;
46 
47  std::vector<double> mHp; // vary charged Higgs mass
48  mHp.push_back(300.);
49  mHp.push_back(500.);
50  mHp.push_back(750.);
51  mHp.push_back(1000.);
52  mHp.push_back(1500.);
53 
54  // Number of events to generate (per H+ mass)
55  int Ngen = 10000;
56 
58  g_Log << LogInfo << "Initializing generator frames and tree..." << LogEnd;
60  ppLabGenFrame LAB_Gen("LAB_Gen","LAB");
61  DecayGenFrame Hp_Gen("Hp_Gen","H^{ +}");
62  DecayGenFrame H_Gen("H_Gen","h^{ 0}");
63  ResonanceGenFrame W_Gen("W_Gen","W^{ +}");
64  VisibleGenFrame G1_Gen("G1_Gen","#gamma_{1}");
65  VisibleGenFrame G2_Gen("G2_Gen","#gamma_{2}");
66  VisibleGenFrame L_Gen("L_Gen","#it{l}^{ +}");
67  InvisibleGenFrame NU_Gen("NU_Gen","#nu");
68 
69  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
70 
71  LAB_Gen.SetChildFrame(Hp_Gen);
72  Hp_Gen.AddChildFrame(H_Gen);
73  Hp_Gen.AddChildFrame(W_Gen);
74  H_Gen.AddChildFrame(G1_Gen);
75  H_Gen.AddChildFrame(G2_Gen);
76  W_Gen.AddChildFrame(L_Gen);
77  W_Gen.AddChildFrame(NU_Gen);
78 
79  if(LAB_Gen.InitializeTree())
80  g_Log << LogInfo << "...Successfully initialized generator tree" << LogEnd;
81  else
82  g_Log << LogError << "...Failed initializing generator tree" << LogEnd;
83 
84  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
85 
86  if(mHp.size() < 1) return;
87 
88  // set neutral and charged Higgs mass
89  Hp_Gen.SetMass(mHp[0]); H_Gen.SetMass(mH);
90  // set W mass and width
91  W_Gen.SetMass(mW); W_Gen.SetWidth(wW);
92 
93  // set photon/lepton pT and eta cuts
94  L_Gen.SetPtCut(15.); L_Gen.SetEtaCut(2.5);
95  G1_Gen.SetPtCut(20.); G1_Gen.SetEtaCut(3.);
96  G2_Gen.SetPtCut(20.); G2_Gen.SetEtaCut(3.);
97 
98  if(LAB_Gen.InitializeAnalysis())
99  g_Log << LogInfo << "...Successfully initialized generator analysis" << std::endl << LogEnd;
100  else
101  g_Log << LogError << "...Failed initializing generator analysis" << LogEnd;
104 
106  g_Log << LogInfo << "Initializing reconstruction frames and trees..." << LogEnd;
108  LabRecoFrame LAB("LAB","LAB");
109  DecayRecoFrame Hp("Hp","H^{ +}");
110  DecayRecoFrame H("H","h^{ 0}");
111  DecayRecoFrame W("W","W^{ +}");
112  VisibleRecoFrame G1("G1","#gamma_{1}");
113  VisibleRecoFrame G2("G2","#gamma_{2}");
114  VisibleRecoFrame L("L","#it{l}^{ +}");
115  InvisibleRecoFrame NU("NU","#nu");
116 
117  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
118 
119  LAB.SetChildFrame(Hp);
120  Hp.AddChildFrame(H);
121  Hp.AddChildFrame(W);
122  H.AddChildFrame(G1);
123  H.AddChildFrame(G2);
124  W.AddChildFrame(L);
125  W.AddChildFrame(NU);
126 
127  if(LAB.InitializeTree())
128  g_Log << LogInfo << "...Successfully initialized reconstruction trees" << LogEnd;
129  else
130  g_Log << LogError << "...Failed initializing reconstruction trees" << LogEnd;
131 
132  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
133 
134  // Invisible Group
135  InvisibleGroup INV("INV","#nu Jigsaws");
136  INV.AddFrame(NU);
137 
138  // Set neutrino masses to zero
139  SetMassInvJigsaw NuM("NuM","M_{#nu} = 0");
140  INV.AddJigsaw(NuM);
141 
142  // Set neutrino rapidity to that of visible particles
143  SetRapidityInvJigsaw NuR("NuR","#eta_{#nu} = #eta_{#gamma#gamma#it{l}}");
144  INV.AddJigsaw(NuR);
145  NuR.AddVisibleFrames(Hp.GetListVisibleFrames());
146 
147  if(LAB.InitializeAnalysis())
148  g_Log << LogInfo << "...Successfully initialized analyses" << LogEnd;
149  else
150  g_Log << LogError << "...Failed initializing analyses" << LogEnd;
151 
154 
155  TreePlot* treePlot = new TreePlot("TreePlot","TreePlot");
156 
157  treePlot->SetTree(LAB_Gen);
158  treePlot->Draw("GenTree", "Generator Tree", true);
159 
160  treePlot->SetTree(LAB);
161  treePlot->Draw("RecoTree", "Reconstruction Tree");
162 
163  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
164 
165  // Declare observables for histogram booking
166  std::string plot_title = "pp #rightarrow H^{ +} #rightarrow h^{ 0}(#gamma #gamma ) W(#it{l} #nu)";
167  HistPlot* histPlot = new HistPlot("HistPlot", plot_title);
168 
169  RFList<const HistPlotCategory> cat_list;
170  int Nmass = mHp.size();
171  for(int m = 0; m < Nmass; m++){
172  char smass[50], scat[50];
173  sprintf(scat, "MHp%.0f", mHp[m]);
174  sprintf(smass, "m_{H^{ +}} = %.0f", mHp[m]);
175  cat_list += histPlot->GetNewCategory(scat, smass);
176  }
177 
178  const HistPlotCategory& cat_Hp = histPlot->GetNewCategory("HprodHp", "h^{ 0} prod. frame = H^{ +}");
179  const HistPlotCategory& cat_LAB = histPlot->GetNewCategory("HprodLAB", "h^{ 0} prod. frame = LAB");
180 
181  const HistPlotVar& MHp = histPlot->GetNewVar("MHp", "M_{H^{ +}}", 0., 2400., "[GeV]");
182  const HistPlotVar& MHpN = histPlot->GetNewVar("MHpN", "M_{H^{ +}} / m_{H^{ +}}^{true}", 0.7, 1.05);
183  const HistPlotVar& MWN = histPlot->GetNewVar("MWN", "M_{W} / m_{W}^{true}", 0., 2.);
184  const HistPlotVar& cosHp = histPlot->GetNewVar("cosHp","cos #theta_{H^{ +}}", -1., 1.);
185  const HistPlotVar& cosW = histPlot->GetNewVar("cosW","cos #theta_{W}", -1., 1.);
186  const HistPlotVar& cosH = histPlot->GetNewVar("cosH","cos #theta_{h^{ 0}}", -1., 1.);
187  const HistPlotVar& DcosHp = histPlot->GetNewVar("DcosHp","#theta_{H^{ +}} - #theta_{H^{ +}}^{true}", -1., 1.);
188  const HistPlotVar& DcosW = histPlot->GetNewVar("DcosW","#theta_{W} - #theta_{W}^{true}", -1., 1.);
189  const HistPlotVar& DcosH = histPlot->GetNewVar("DcosW","#theta_{h^{ 0}} - #theta_{h^{ 0}}^{true}", -1., 1.);
190 
191  histPlot->AddPlot(DcosH, cat_list);
192  histPlot->AddPlot(DcosW, cat_list);
193  histPlot->AddPlot(DcosHp, cat_list);
194  histPlot->AddPlot(MWN, cat_list);
195  histPlot->AddPlot(MHpN, cat_list);
196  histPlot->AddPlot(MHp, cat_list);
197  histPlot->AddPlot(MWN, MHpN, cat_list[2]);
198  histPlot->AddPlot(DcosW, MWN, cat_list[2]);
199  histPlot->AddPlot(DcosHp, MHpN, cat_list[2]);
200  histPlot->AddPlot(DcosHp, DcosW, cat_list[2]);
201  histPlot->AddPlot(DcosH, MHpN, cat_list[2]);
202 
203  histPlot->AddPlot(DcosH, cat_Hp+cat_LAB);
204 
207 
208  for(int m = 0; m < Nmass; m++){
209  g_Log << LogInfo << "Generating events for H^{+} mass = " << mHp[m] << LogEnd;
210 
211  Hp_Gen.SetMass(mHp[m]);
212  LAB_Gen.InitializeAnalysis();
213 
214  for(int igen = 0; igen < Ngen; igen++){
215  if(igen%((std::max(Ngen,10))/10) == 0)
216  g_Log << LogInfo << "Generating event " << igen << " of " << Ngen << LogEnd;
217 
218  // generate event
219  LAB_Gen.ClearEvent(); // clear the gen tree
220 
221  LAB_Gen.SetPToverM(gRandom->Rndm()); // give charged Higgs some Pt
222 
223  LAB_Gen.AnalyzeEvent(); // generate a new event
224 
225  TVector3 MET = LAB_Gen.GetInvisibleMomentum(); // Get the MET from gen tree
226  MET.SetZ(0.);
227 
228  // analyze event
229  LAB.ClearEvent(); // clear the reco tree
230 
231  L.SetLabFrameFourVector(L_Gen.GetFourVector()); // Set lepton 4-vec
232  G1.SetLabFrameFourVector(G1_Gen.GetFourVector()); // Set photons' 4-vec
233  G2.SetLabFrameFourVector(G2_Gen.GetFourVector());
234  INV.SetLabFrameThreeVector(MET); // Set the MET in reco tree
235 
236  LAB.AnalyzeEvent(); //analyze the event
237 
238  // Generator-level observables
239  double MHpgen = Hp_Gen.GetMass();
240  double MWgen = W_Gen.GetMass();
241  double cosHpgen = Hp_Gen.GetCosDecayAngle();
242  double cosHgen = H_Gen.GetCosDecayAngle();
243  double cosWgen = W_Gen.GetCosDecayAngle();
244 
245  // Reconstructed observables
246  MHp = Hp.GetMass();
247  MHpN = Hp.GetMass()/MHpgen;
248  MWN = W.GetMass()/MWgen;
249  cosHp = Hp.GetCosDecayAngle();
250  cosH = H.GetCosDecayAngle();
251  cosW = W.GetCosDecayAngle();
252  DcosHp = asin(sqrt(1.-cosHp*cosHp)*cosHpgen-sqrt(1.-cosHpgen*cosHpgen)*cosHp);
253  DcosH = asin(sqrt(1.-cosH*cosH)*cosHgen-sqrt(1.-cosHgen*cosHgen)*cosH);
254  DcosW = asin(sqrt(1.-cosW*cosW)*cosWgen-sqrt(1.-cosWgen*cosWgen)*cosW);
255 
256  histPlot->Fill(cat_list[m]);
257 
258  if(m == 0){
259  histPlot->Fill(cat_Hp);
260 
261  TVector3 Hboost = H.GetFourVector().BoostVector();
262  TLorentzVector vP_G1 = G1.GetFourVector();
263  vP_G1.Boost(-Hboost);
264  cosH = -vP_G1.Vect().Unit().Dot(Hboost.Unit());
265  DcosH = asin(sqrt(1.-cosH*cosH)*cosHgen-sqrt(1.-cosHgen*cosHgen)*cosH);
266 
267  histPlot->Fill(cat_LAB);
268  }
269  }
270 
271  LAB_Gen.PrintGeneratorEfficiency();
272  }
273 
274  histPlot->Draw();
275 
276  TFile fout(output_name.c_str(),"RECREATE");
277  fout.Close();
278  histPlot->WriteOutput(output_name);
279  histPlot->WriteHist(output_name);
280  treePlot->WriteOutput(output_name);
281 
282  g_Log << LogInfo << "Finished" << LogEnd;
283 
284 }
285 
286 # ifndef __CINT__ // main function for stand-alone compilation
287 int main(){
288  example_Hp_to_HggWlnu();
289  return 0;
290 }
291 #endif
Generated on Fri Aug 12 2016 14:55:13 for RestFrames by   doxygen 1.8.9.1