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