LOGO
RestFrames  v1.0.0
RestFrames HEP Event Analysis Software Library
example_X2X2_to_ZllXHggX.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_X2X2_to_ZllXHggX(std::string output_name =
40  "output_X2X2_to_ZllXHggX.root"){
41 
42  // set particle masses and widths [GeV]
43  double mX2 = 500.;
44  double mZ = 91.19; // PDG 2016
45  double wZ = 2.50;
46  double mH = 125.;
47  double wH = 0.04;
48 
49  // number of different neutralino masses to evaluate
50  int NmX1 = 4;
51  std::vector<double> mX1;
52  mX1.push_back(400.); // lightest X1 mass to evaluate
53 
54  // Number of events to generate
55  int Ngen = 10000;
56 
58  g_Log << LogInfo << "Initializing generator frames and tree..." << LogEnd;
60  ppLabGenFrame LAB_Gen("LAB_Gen","LAB");
61  DecayGenFrame X2X2_Gen("X2X2_Gen","#tilde{#chi}^{ 0}_{2} #tilde{#chi}^{ 0}_{2}");
62  DecayGenFrame X2a_Gen("X2a_Gen","#tilde{#chi}^{ 0}_{2 a}");
63  DecayGenFrame X2b_Gen("X2b_Gen","#tilde{#chi}^{ 0}_{2 b}");
64  ResonanceGenFrame Za_Gen("Za_Gen","Z_{a}");
65  ResonanceGenFrame Hb_Gen("Hb_Gen","H_{b}");
66  VisibleGenFrame L1_Gen("L1_Gen","#it{l}_{1}");
67  VisibleGenFrame L2_Gen("L2_Gen","#it{l}_{2}");
68  VisibleGenFrame G1_Gen("G1_Gen","#gamma_{1}");
69  VisibleGenFrame G2_Gen("G2_Gen","#gamma_{2}");
70  InvisibleGenFrame X1a_Gen("X1a_Gen","#tilde{#chi}^{ 0}_{1 a}");
71  InvisibleGenFrame X1b_Gen("X1b_Gen","#tilde{#chi}^{ 0}_{1 b}");
72 
73  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
74 
75  LAB_Gen.SetChildFrame(X2X2_Gen);
76  X2X2_Gen.AddChildFrame(X2a_Gen);
77  X2X2_Gen.AddChildFrame(X2b_Gen);
78  X2a_Gen.AddChildFrame(Za_Gen);
79  X2a_Gen.AddChildFrame(X1a_Gen);
80  X2b_Gen.AddChildFrame(Hb_Gen);
81  X2b_Gen.AddChildFrame(X1b_Gen);
82  Za_Gen.AddChildFrame(L1_Gen);
83  Za_Gen.AddChildFrame(L2_Gen);
84  Hb_Gen.AddChildFrame(G1_Gen);
85  Hb_Gen.AddChildFrame(G2_Gen);
86 
87  if(LAB_Gen.InitializeTree())
88  g_Log << LogInfo << "...Successfully initialized generator tree" << LogEnd;
89  else
90  g_Log << LogError << "...Failed initializing generator tree" << LogEnd;
91 
92  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
93 
94  if(NmX1 < 1 || mX1.size() < 1)
95  return;
96  if(mX1[0] >= mX2)
97  return;
98 
99  X2X2_Gen.SetVariableMass();
100  X2a_Gen.SetMass(mX2); X2b_Gen.SetMass(mX2);
101  X1a_Gen.SetMass(mX1[0]); X1b_Gen.SetMass(mX1[0]);
102  Za_Gen.SetMass(mZ);
103  Za_Gen.SetWidth(wZ);
104  Hb_Gen.SetMass(mH);
105  Hb_Gen.SetWidth(wH);
106 
107  L1_Gen.SetPtCut(8.); L1_Gen.SetEtaCut(2.5);
108  L2_Gen.SetPtCut(8.); L2_Gen.SetEtaCut(2.5);
109  G1_Gen.SetPtCut(20.); G1_Gen.SetEtaCut(3.0);
110  G2_Gen.SetPtCut(20.); G2_Gen.SetEtaCut(3.0);
111 
112  if(LAB_Gen.InitializeAnalysis())
113  g_Log << LogInfo << "...Successfully initialized generator analysis" << LogEnd;
114  else
115  g_Log << LogError << "...Failed initializing generator analysis" << LogEnd;
118 
120  g_Log << LogInfo << "Initializing reconstruction frames and tree..." << LogEnd;
122  LabRecoFrame LAB("LAB","LAB");
123  DecayRecoFrame X2X2("X2X2","#tilde{#chi}^{ 0}_{2} #tilde{#chi}^{ 0}_{2}");
124  DecayRecoFrame X2a("X2a","#tilde{#chi}^{ 0}_{2 a}");
125  DecayRecoFrame X2b("X2b","#tilde{#chi}^{ 0}_{2 b}");
126  DecayRecoFrame Za("Za","Z_{a}");
127  DecayRecoFrame Hb("Hb","H_{b}");
128  VisibleRecoFrame L1("L1","#it{l}_{1}");
129  VisibleRecoFrame L2("L2","#it{l}_{2}");
130  VisibleRecoFrame G1("G1","#gamma_{1}");
131  VisibleRecoFrame G2("G2","#gamma_{2}");
132  InvisibleRecoFrame X1a("X1a","#tilde{#chi}^{ 0}_{1 a}");
133  InvisibleRecoFrame X1b("X1b","#tilde{#chi}^{ 0}_{1 b}");
134 
135  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
136 
137  LAB.SetChildFrame(X2X2);
138  X2X2.AddChildFrame(X2a);
139  X2X2.AddChildFrame(X2b);
140  X2a.AddChildFrame(Za);
141  X2a.AddChildFrame(X1a);
142  X2b.AddChildFrame(Hb);
143  X2b.AddChildFrame(X1b);
144  Za.AddChildFrame(L1);
145  Za.AddChildFrame(L2);
146  Hb.AddChildFrame(G1);
147  Hb.AddChildFrame(G2);
148 
149  if(LAB.InitializeTree())
150  g_Log << LogInfo << "...Successfully initialized reconstruction trees" << LogEnd;
151  else
152  g_Log << LogError << "...Failed initializing reconstruction trees" << LogEnd;
153 
154  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
155 
156  // Invisible Groups
157  InvisibleGroup INV("INV","#tilde{#chi}_{1}^{ 0} Jigsaws");
158  INV.AddFrames(X1a+X1b);
159 
160  // Set di-LSP mass to minimum Lorentz-invariant expression
161  SetMassInvJigsaw X1_mass("X1_mass","Set M_{#tilde{#chi}_{1}^{ 0} #tilde{#chi}_{1}^{ 0}} to minimum");
162  INV.AddJigsaw(X1_mass);
163 
164  // Set di-LSP rapidity to that of visible particles
165  SetRapidityInvJigsaw X1_eta("X1_eta","#eta_{#tilde{#chi}_{1}^{ 0} #tilde{#chi}_{1}^{ 0}} = #eta_{2#gamma+2#it{l}}");
166  INV.AddJigsaw(X1_eta);
167  X1_eta.AddVisibleFrames(X2X2.GetListVisibleFrames());
168 
169  ContraBoostInvJigsaw X1X1_contra("X1X1_contra","Contraboost invariant Jigsaw");
170  //MinMassDiffInvJigsaw X1X1_contra("MinMh_R","min M_{h}, M_{h}^{ a}= M_{h}^{ b}",2);
171  //MinMassesSqInvJigsaw X1X1_contra("MinMWa_R","min M_{W}, M_{W}^{a,a}= M_{W}^{a,b}", 2);
172  INV.AddJigsaw(X1X1_contra);
173  X1X1_contra.AddVisibleFrames(X2a.GetListVisibleFrames(), 0);
174  X1X1_contra.AddVisibleFrames(X2b.GetListVisibleFrames(), 1);
175  X1X1_contra.AddInvisibleFrame(X1a, 0);
176  X1X1_contra.AddInvisibleFrame(X1b, 1);
177 
178  if(LAB.InitializeAnalysis())
179  g_Log << LogInfo << "...Successfully initialized analysis" << std::endl << LogEnd;
180  else
181  g_Log << LogError << "...Failed initializing analysis" << LogEnd;
184 
185  TreePlot* treePlot = new TreePlot("TreePlot","TreePlot");
186 
187  treePlot->SetTree(LAB_Gen);
188  treePlot->Draw("GenTree", "Generator Tree", true);
189 
190  treePlot->SetTree(LAB);
191  treePlot->Draw("RecoTree", "Reconstruction Tree");
192 
193  treePlot->SetTree(INV);
194  treePlot->Draw("InvTree", "Inivisible Jigsaws", true);
195 
196  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
197 
198  // Declare observables for histogram booking
199  HistPlot* histPlot = new HistPlot("Plots",
200  std::string("#tilde{#chi}_{2}^{ 0} #tilde{#chi}_{2}^{ 0}") +
201  "#rightarrow Z(#it{l}#it{l}) #tilde{#chi}_{1}^{ 0}"+
202  "H(#gamma #gamma) #tilde{#chi}_{1}^{ 0}");
203 
204  RFList<const HistPlotCategory> cat_list;
205  char smassX2[200];
206  sprintf(smassX2, "m_{#tilde{#chi}^{0}_{2}}= %.0f", mX2);
207  for(int m = 0; m < NmX1; m++){
208  if(m > 0)
209  mX1.push_back(mX1[0] + m*(mX2-mX1[0])/double(NmX1));
210 
211  char smassX1[200], scat[50];
212  sprintf(scat, "MX1_%.0f", mX1[m]);
213  sprintf(smassX1, "m_{#tilde{#chi}^{0}_{1}}= %.0f", mX1[m]);
214 
215  cat_list += histPlot->GetNewCategory(scat, std::string(smassX2)+
216  " , "+std::string(smassX1));
217  }
218 
219  const HistPlotVar& MCM = histPlot->GetNewVar("MCM",
220  "M_{#tilde{#chi}_{2}^{ 0} #tilde{#chi}_{2}^{ 0}}",
221  0., 2.);
222  const HistPlotVar& EZX2a = histPlot->GetNewVar("EZX2a", "E_{Z}^{ #tilde{#chi}_{2 a}^{ 0}}", 0., 2.);
223  const HistPlotVar& EHX2b = histPlot->GetNewVar("EHX2b", "E_{H}^{ #tilde{#chi}_{2 b}^{ 0}}", 0., 2.);
224  const HistPlotVar& cosX2a = histPlot->GetNewVar("cosX2a","cos #theta_{#tilde{#chi}_{2 a}^{ 0}}", -1., 1.);
225  const HistPlotVar& cosX2b = histPlot->GetNewVar("cosX2b","cos #theta_{#tilde{#chi}_{2 b}^{ 0}}", -1., 1.);
226  const HistPlotVar& cosZ = histPlot->GetNewVar("cosZ","cos #theta_{Z}", -1., 1.);
227  const HistPlotVar& cosH = histPlot->GetNewVar("cosH","cos #theta_{H}", -1., 1.);
228  const HistPlotVar& DcosZ = histPlot->GetNewVar("DcosZ","#theta_{Z} - #theta_{Z}^{gen}", -1., 1.);
229  const HistPlotVar& DcosH = histPlot->GetNewVar("DcosH","#theta_{H} - #theta_{H}^{gen}", -1., 1.);
230  const HistPlotVar& DcosX2a = histPlot->GetNewVar("DcosX2a","#theta_{X2a} - #theta_{X2a}^{gen}", -1., 1.);
231  const HistPlotVar& DcosX2b = histPlot->GetNewVar("DcosX2b","#theta_{X2b} - #theta_{X2b}^{gen}", -1., 1.);
232  const HistPlotVar& RISR = histPlot->GetNewVar("RISR","R_{ISR}", 0., 1.5);
233  const HistPlotVar& PTISR = histPlot->GetNewVar("PTISR","p_{T}^{ ISR}", 0., 1000., "[GeV]");
234 
235  histPlot->AddPlot(MCM, cat_list);
236  histPlot->AddPlot(EZX2a, cat_list);
237  histPlot->AddPlot(EHX2b, cat_list);
238  histPlot->AddPlot(cosZ, cat_list);
239  histPlot->AddPlot(cosH, cat_list);
240  histPlot->AddPlot(DcosZ, cat_list);
241  histPlot->AddPlot(DcosH, cat_list);
242  histPlot->AddPlot(DcosX2a, cat_list);
243  histPlot->AddPlot(DcosX2b, cat_list);
244  histPlot->AddPlot(MCM, EZX2a, cat_list[NmX1/2]);
245  histPlot->AddPlot(MCM, EHX2b, cat_list[NmX1/2]);
246  histPlot->AddPlot(EZX2a, EHX2b, cat_list[NmX1/2]);
247  histPlot->AddPlot(EZX2a, DcosX2a, cat_list[NmX1/2]);
248  histPlot->AddPlot(EHX2b, DcosX2b, cat_list[NmX1/2]);
249  histPlot->AddPlot(cosZ, DcosZ, cat_list[NmX1/2]);
250  histPlot->AddPlot(cosH, DcosH, cat_list[NmX1/2]);
251 
252  histPlot->AddPlot(RISR, cat_list);
253  histPlot->AddPlot(RISR, PTISR, cat_list[NmX1/2]);
254  histPlot->AddPlot(RISR, EZX2a, cat_list[NmX1/2]);
255 
256  for(int m = 0; m < NmX1; m++){
257  g_Log << LogInfo << "Generating events for ";
258  g_Log << "mX2 = " << mX2 << " , ";
259  g_Log << "mX1 = " << mX1[m] << LogEnd;
260 
261  X1a_Gen.SetMass(mX1[m]);
262  X1b_Gen.SetMass(mX1[m]);
263 
264  LAB_Gen.InitializeAnalysis();
265 
266  for(int igen = 0; igen < Ngen; igen++){
267  if(igen%((std::max(Ngen,10))/10) == 0)
268  g_Log << LogInfo << "Generating event " << igen << " of " << Ngen << LogEnd;
269 
270  // generate event
271  LAB_Gen.ClearEvent(); // clear the gen tree
272 
273  //LAB_Gen.SetTransverseMomentum(800.); // give X2X2 some Pt
274  LAB_Gen.SetTransverseMomentum(1000.*gRandom->Rndm()); // give X2X2 some Pt
275 
276  LAB_Gen.AnalyzeEvent(); // generate a new event
277 
278  // analyze event
279  LAB.ClearEvent(); // clear the reco tree
280  L1.SetLabFrameFourVector(L1_Gen.GetFourVector()); // Set lepton 4-vectors
281  L2.SetLabFrameFourVector(L2_Gen.GetFourVector());
282  G1.SetLabFrameFourVector(G1_Gen.GetFourVector()); // Set photon 4-vectors
283  G2.SetLabFrameFourVector(G2_Gen.GetFourVector());
284  TVector3 MET = LAB_Gen.GetInvisibleMomentum(); // Get the MET from gen tree
285  MET.SetZ(0.);
286  INV.SetLabFrameThreeVector(MET); // Set the MET in reco tree
287  LAB.AnalyzeEvent(); //analyze the event
288 
289  // calculate observables
290  MCM = X2X2.GetMass() / X2X2_Gen.GetMass();
291 
292  EZX2a = Za.GetEnergy(X2a) / Za_Gen.GetEnergy(X2a_Gen);
293  EHX2b = Hb.GetEnergy(X2b) / Hb_Gen.GetEnergy(X2b_Gen);
294 
295  cosX2a = X2a.GetCosDecayAngle();
296  double cosX2agen = X2a_Gen.GetCosDecayAngle();
297  cosX2b = X2b.GetCosDecayAngle();
298  double cosX2bgen = X2b_Gen.GetCosDecayAngle();
299  cosZ = Za.GetCosDecayAngle();
300  double cosZgen = Za_Gen.GetCosDecayAngle();
301  cosH = Hb.GetCosDecayAngle();
302  double cosHgen = Hb_Gen.GetCosDecayAngle();
303 
304  DcosX2a = asin(sqrt(1.-cosX2a*cosX2a)*cosX2agen-sqrt(1.-cosX2agen*cosX2agen)*cosX2a);
305  DcosX2b = asin(sqrt(1.-cosX2b*cosX2b)*cosX2bgen-sqrt(1.-cosX2bgen*cosX2bgen)*cosX2b);
306  DcosZ = asin(sqrt(1.-cosZ*cosZ)*cosZgen-sqrt(1.-cosZgen*cosZgen)*cosZ);
307  DcosH = asin(sqrt(1.-cosH*cosH)*cosHgen-sqrt(1.-cosHgen*cosHgen)*cosH);
308 
309  TVector3 vP_ISR = X2X2.GetFourVector(LAB).Vect();
310  TVector3 vP_I = X2X2.GetListInvisibleFrames().GetFourVector(LAB).Vect();
311  vP_ISR.SetZ(0.);
312  vP_I.SetZ(0.);
313 
314  RISR = fabs(vP_I.Dot(vP_ISR.Unit())) / vP_ISR.Mag();
315  PTISR = vP_ISR.Mag();
316 
317  histPlot->Fill(cat_list[m]);
318  }
319 
320  LAB_Gen.PrintGeneratorEfficiency();
321  }
322 
323  histPlot->Draw();
324 
325  TFile fout(output_name.c_str(),"RECREATE");
326  fout.Close();
327  histPlot->WriteOutput(output_name);
328  histPlot->WriteHist(output_name);
329  treePlot->WriteOutput(output_name);
330 
331  g_Log << LogInfo << "Finished" << LogEnd;
332 }
333 
334 # ifndef __CINT__ // main function for stand-alone compilation
335  int main(){
336  example_X2X2_to_ZllXHggX();
337  return 0;
338  }
339 #endif
Generated on Fri Aug 12 2016 14:55:13 for RestFrames by   doxygen 1.8.9.1