LOGO
RestFrames  v1.0.0
RestFrames HEP Event Analysis Software Library
example_DiGluino_to_bbXbbX.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_DiGluino_to_bbXbbX(std::string output_name =
40  "output_DiGluino_to_bbXbbX.root"){
41  double mG = 1000.;
42  double mX = 100.;
43 
44  // Number of events to generate
45  int Ngen = 10000;
46 
48  g_Log << LogInfo << "Initializing generator frames and tree..." << LogEnd;
50  ppLabGenFrame LAB_Gen("LAB_Gen","LAB");
51  DecayGenFrame GG_Gen("GG_Gen","#tilde{g}#tilde{g}");
52  DecayGenFrame Ga_Gen("Ga_Gen","#tilde{g}_{a}");
53  DecayGenFrame Gb_Gen("Gb_Gen","#tilde{g}_{b}");
54  VisibleGenFrame V1a_Gen("V1a_Gen","j_{1a}");
55  VisibleGenFrame V2a_Gen("V2a_Gen","j_{2a}");
56  InvisibleGenFrame Xa_Gen("Xa_Gen","#tilde{#chi}_{a}");
57  VisibleGenFrame V1b_Gen("V1b_Gen","j_{1b}");
58  VisibleGenFrame V2b_Gen("V2b_Gen","j_{2b}");
59  InvisibleGenFrame Xb_Gen("Xb_Gen","#tilde{#chi}_{b}");
60 
61  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
62 
63  LAB_Gen.SetChildFrame(GG_Gen);
64  GG_Gen.AddChildFrame(Ga_Gen);
65  GG_Gen.AddChildFrame(Gb_Gen);
66  Ga_Gen.AddChildFrame(V1a_Gen);
67  Ga_Gen.AddChildFrame(V2a_Gen);
68  Ga_Gen.AddChildFrame(Xa_Gen);
69  Gb_Gen.AddChildFrame(V1b_Gen);
70  Gb_Gen.AddChildFrame(V2b_Gen);
71  Gb_Gen.AddChildFrame(Xb_Gen);
72 
73  if(LAB_Gen.InitializeTree())
74  g_Log << LogInfo << "...Successfully initialized generator tree" << LogEnd;
75  else
76  g_Log << LogError << "...Failed initializing generator tree" << LogEnd;
77 
78  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
79 
80  GG_Gen.SetVariableMass();
81  // set gluino masses
82  Ga_Gen.SetMass(mG);
83  Gb_Gen.SetMass(mG);
84  // set X masses
85  Xa_Gen.SetMass(mX);
86  Xb_Gen.SetMass(mX);
87  // set jet pT and eta cuts
88  V1a_Gen.SetPtCut(30.); V1a_Gen.SetEtaCut(2.5);
89  V1b_Gen.SetPtCut(30.); V1b_Gen.SetEtaCut(2.5);
90  V2a_Gen.SetPtCut(30.); V2a_Gen.SetEtaCut(2.5);
91  V2b_Gen.SetPtCut(30.); V2b_Gen.SetEtaCut(2.5);
92 
93  if(LAB_Gen.InitializeAnalysis())
94  g_Log << LogInfo << "...Successfully initialized generator analysis" << std::endl << LogEnd;
95  else
96  g_Log << LogError << "...Failed initializing generator analysis" << LogEnd;
99 
101  g_Log << LogInfo << "Initializing reconstruction frames and trees..." << LogEnd;
103  LabRecoFrame LAB("LAB","LAB");
104  DecayRecoFrame GG("GG","#tilde{g}#tilde{g}");
105  DecayRecoFrame Ga("Ga","#tilde{g}_{a}");
106  DecayRecoFrame Gb("Gb","#tilde{g}_{b}");
107  DecayRecoFrame Ca("Ca","C_{a}");
108  DecayRecoFrame Cb("Cb","C_{b}");
109  VisibleRecoFrame V1a("V1a","j_{1a}");
110  VisibleRecoFrame V2a("V2a","j_{2a}");
111  InvisibleRecoFrame Xa("Xa","#tilde{#chi}_{a}");
112  VisibleRecoFrame V1b("V1b","j_{1b}");
113  VisibleRecoFrame V2b("V2b","j_{2b}");
114  InvisibleRecoFrame Xb("Xb","#tilde{#chi}_{b}");
115 
116  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
117 
118  LAB.SetChildFrame(GG);
119  GG.AddChildFrame(Ga);
120  GG.AddChildFrame(Gb);
121  Ga.AddChildFrame(V1a);
122  Ga.AddChildFrame(Ca);
123  Ca.AddChildFrame(V2a);
124  Ca.AddChildFrame(Xa);
125  Gb.AddChildFrame(V1b);
126  Gb.AddChildFrame(Cb);
127  Cb.AddChildFrame(V2b);
128  Cb.AddChildFrame(Xb);
129 
130  if(LAB.InitializeTree())
131  g_Log << LogInfo << "...Successfully initialized reconstruction trees" << LogEnd;
132  else
133  g_Log << LogError << "...Failed initializing reconstruction trees" << LogEnd;
134 
135  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
136 
137  // Set up 'background-like' analysis tree
138  LabRecoFrame LAB_B("LAB_B","LAB");
139  SelfAssemblingRecoFrame CM_B("CM_B","CM");
140  VisibleRecoFrame V_B("V_B","Vis");
141  InvisibleRecoFrame I_B("I_B","Inv");
142  LAB_B.SetChildFrame(CM_B);
143  CM_B.AddChildFrame(V_B);
144  CM_B.AddChildFrame(I_B);
145 
146  if(LAB.InitializeTree())
147  g_Log << LogInfo << "...Successfully initialized reconstruction trees" << LogEnd;
148  else
149  g_Log << LogError << "...Failed initializing reconstruction trees" << LogEnd;
150 
151  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
152 
153  InvisibleGroup INV("INV","WIMP Jigsaws");
154  INV.AddFrame(Xa);
155  INV.AddFrame(Xb);
156  CombinatoricGroup VIS("VIS","Visible Object Jigsaws");
157 
158  // visible frames in first decay step must always have at least one element
159  VIS.AddFrame(V1a);
160  VIS.AddFrame(V1b);
161  VIS.SetNElementsForFrame(V1a,1,false);
162  VIS.SetNElementsForFrame(V1b,1,false);
163  // visible frames in second decay step can have zero elements
164  VIS.AddFrame(V2a);
165  VIS.AddFrame(V2b);
166  VIS.SetNElementsForFrame(V2a,0,false);
167  VIS.SetNElementsForFrame(V2b,0,false);
168 
169  InvisibleGroup INV_B("INV_B","Invisible State Jigsaws");
170  INV_B.AddFrame(I_B);
171  CombinatoricGroup VIS_B("VIS_B","Visible Object Jigsaws");
172  VIS_B.AddFrame(V_B);
173  VIS_B.SetNElementsForFrame(V_B,1,false);
174 
175  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
176 
177  // signal-like jigsaws
178  SetMassInvJigsaw MinMassJigsaw("MINMASS", "Invisible system mass Jigsaw");
179  INV.AddJigsaw(MinMassJigsaw);
180  SetRapidityInvJigsaw RapidityJigsaw("RAPIDITY", "Invisible system rapidity Jigsaw");
181  INV.AddJigsaw(RapidityJigsaw);
182  RapidityJigsaw.AddVisibleFrames((LAB.GetListVisibleFrames()));
183  ContraBoostInvJigsaw ContraBoostJigsaw("CONTRA","Contraboost invariant Jigsaw");
184  INV.AddJigsaw(ContraBoostJigsaw);
185  ContraBoostJigsaw.AddVisibleFrames((Ga.GetListVisibleFrames()), 0);
186  ContraBoostJigsaw.AddVisibleFrames((Gb.GetListVisibleFrames()), 1);
187  ContraBoostJigsaw.AddInvisibleFrames((Ga.GetListInvisibleFrames()), 0);
188  ContraBoostJigsaw.AddInvisibleFrames((Gb.GetListInvisibleFrames()), 1);
189  MinMassesCombJigsaw HemiJigsaw("HEM_JIGSAW","Minimize m _{V_{a,b}} Jigsaw");
190  VIS.AddJigsaw(HemiJigsaw);
191  HemiJigsaw.AddFrame(V1a,0);
192  HemiJigsaw.AddFrame(V1b,1);
193  HemiJigsaw.AddFrame(V2a,0);
194  HemiJigsaw.AddFrame(V2b,1);
195  MinMassesCombJigsaw CaHemiJigsaw("CaHEM_JIGSAW","Minimize m _{C_{a}} Jigsaw");
196  VIS.AddJigsaw(CaHemiJigsaw);
197  CaHemiJigsaw.AddFrame(V1a,0);
198  CaHemiJigsaw.AddFrame(V2a,1);
199  CaHemiJigsaw.AddFrame(Xa,1);
200  MinMassesCombJigsaw CbHemiJigsaw("CbHEM_JIGSAW","Minimize m _{C_{b}} Jigsaw");
201  VIS.AddJigsaw(CbHemiJigsaw);
202  CbHemiJigsaw.AddFrame(V1b,0);
203  CbHemiJigsaw.AddFrame(V2b,1);
204  CbHemiJigsaw.AddFrame(Xb,1);
205 
206  // background tree jigsaws
207  SetMassInvJigsaw MinMassJigsaw_B("MINMASS_B","Zero Mass for invisible system");
208  INV_B.AddJigsaw(MinMassJigsaw_B);
209  SetRapidityInvJigsaw RapidityJigsaw_B("RAPIDITY_B","Invisible system rapidity Jigsaw");
210  INV_B.AddJigsaw(RapidityJigsaw_B);
211  RapidityJigsaw_B.AddVisibleFrames((LAB_B.GetListVisibleFrames()));
212 
213  // check reconstruction trees
214  if(LAB.InitializeAnalysis() && LAB_B.InitializeAnalysis())
215  g_Log << LogInfo << "...Successfully initialized analyses" << LogEnd;
216  else
217  g_Log << LogError << "...Failed initializing analyses" << LogEnd;
218 
220  // draw some pictures of our trees
222 
223  TreePlot* tree_plot = new TreePlot("TreePlot","TreePlot");
224 
225  // generator tree
226  tree_plot->SetTree(LAB_Gen);
227  tree_plot->Draw("GenTree", "Generator Tree", true);
228 
229  // signal reco tree
230  tree_plot->SetTree(LAB);
231  tree_plot->AddJigsaw(ContraBoostJigsaw);
232  tree_plot->AddJigsaw(HemiJigsaw);
233  tree_plot->AddJigsaw(CaHemiJigsaw);
234  tree_plot->AddJigsaw(CbHemiJigsaw);
235  tree_plot->Draw("SigRecoTree", "Signal Reconstruction Tree");
236 
237  // background reco tree
238  tree_plot->SetTree(LAB_B);
239  tree_plot->Draw("BkgRecoTree", "Background Reconstruction Tree");
240 
241  // Invisible Jigsaws
242  tree_plot->SetTree(INV);
243  tree_plot->Draw("InvTree", "Invisible Objects Jigsaws");
244 
245  // Visible Jigsaws
246  tree_plot->SetTree(VIS);
247  tree_plot->Draw("VisTree", "Visible Objects Jigsaws");
248 
249 
250  for(int igen = 0; igen < Ngen; igen++){
251  if(igen%((std::max(Ngen,10))/10) == 0)
252  g_Log << LogInfo << "Generating event " << igen << " of " << Ngen << LogEnd;
253 
254  // generate event
255  LAB_Gen.ClearEvent(); // clear the gen tree
256 
257  LAB_Gen.SetPToverM(LAB_Gen.GetRandom()); // give the di-gluinos some Pt
258 
259  LAB_Gen.AnalyzeEvent(); // generate a new event
260 
261  // analyze event
262  TVector3 MET = LAB_Gen.GetInvisibleMomentum(); // Get the MET from gen tree
263  MET.SetZ(0.);
264  std::vector<TLorentzVector> JETS; // Get the Jets from gen tree
265  JETS.push_back(V1a_Gen.GetFourVector());
266  JETS.push_back(V2a_Gen.GetFourVector());
267  JETS.push_back(V1b_Gen.GetFourVector());
268  JETS.push_back(V2b_Gen.GetFourVector());
269 
270  // give the signal-like tree the event info and analyze
271  LAB.ClearEvent(); // clear the signal-like tree
272  INV.SetLabFrameThreeVector(MET); // Set the MET in reco tree
273  std::vector<RFKey> jetID; // ID for tracking jets in tree
274  for(int i = 0; i < int(JETS.size()); i++)
275  jetID.push_back(VIS.AddLabFrameFourVector(JETS[i]));
276  LAB.AnalyzeEvent(); // analyze the event
277 
278  // give the background-like tree the event info and analyze
279  LAB_B.ClearEvent(); // clear the bkg-like tree
280  INV_B.SetLabFrameThreeVector(MET); // Set the MET in tree
281  for(int i = 0; i < int(JETS.size()); i++){
282  TLorentzVector jet = JETS[i];
283  jet.SetPtEtaPhiM(jet.Pt(), 0., jet.Phi(), jet.M()); // only pass transverse info to bkg-like tree
284  VIS_B.AddLabFrameFourVector(jet);
285  }
286  LAB_B.AnalyzeEvent(); // analyze the event
287 
288  DecayRecoFrame* G[2];
289  DecayRecoFrame* C[2];
290  VisibleRecoFrame* VS[2];
291  VisibleRecoFrame* VC[2];
292  InvisibleRecoFrame* X[2];
293  // Randomize the two hemispheres
294  int flip = (gRandom->Rndm() > 0.5);
295  G[flip] = &Ga;
296  G[(flip+1)%2] = &Gb;
297  C[flip] = &Ca;
298  C[(flip+1)%2] = &Cb;
299  VS[flip] = &V1a;
300  VS[(flip+1)%2] = &V1b;
301  VC[flip] = &V2a;
302  VC[(flip+1)%2] = &V2b;
303  X[flip] = &Xa;
304  X[(flip+1)%2] = &Xb;
305 
307  // Observable Calculations
309 
310  //
311  // signal tree observables
312  //
313 
314  //*** total CM mass
315  double shat = GG.GetMass();
316  //*** 'mass-less' gluino gamma in CM frame
317  double gaminv = GG.GetVisibleShape();
318 
319  TVector3 vPGG = GG.GetFourVector(LAB).Vect();
320 
321  //*** ratio of CM pT to CM mass
322  double RPT = vPGG.Pt() / (vPGG.Pt() + shat/4.);
323  //*** ratio of CM pz to CM mass
324  double RPZ = vPGG.Pz() / (vPGG.Pz() + shat/4.);
325  //*** cos decay angle of GG system
326  double cosGG = GG.GetCosDecayAngle();
327  //*** delta phi between lab and GG decay planes
328  double dphiLGG = LAB.GetDeltaPhiDecayPlanes(GG);
329 
330  TLorentzVector vV1 = G[0]->GetVisibleFourVector(*G[0]);
331  TLorentzVector vV2 = G[1]->GetVisibleFourVector(*G[1]);
332 
333  //*** gluino mass
334  double MG = (vV1.M2()-vV2.M2())/(2.*(vV1.E()-vV2.E()));
335 
336  double PG = G[0]->GetMomentum(GG);
337  double MGG = 2.*sqrt(PG*PG + MG*MG);
338  double gaminvGG = 2.*MG/MGG;
339  double beta = sqrt(1.- gaminv*gaminv);
340  double betaGG = sqrt(1.- gaminvGG*gaminvGG);
341 
342  //*** velocity difference between 'massive' and 'mass-less'
343  double DeltaBetaGG = -(betaGG-beta)/(1.-betaGG*beta);
344  //*** delta phi between GG visible decay products and GG decay axis
345  double dphiVG = GG.GetDeltaPhiDecayVisible();
346  //*** delta phi between GG visible decay products and GG momentum
347  double dphiVGG = GG.GetDeltaPhiBoostVisible();
348 
349  // 'hemisphere' (one for each 'gluino') observables
350 
351  //*** number of visible objects (jets) in hemisphere
352  double NV[2];
353  //*** cosine gluino decay angle
354  double cosG[2];
355  //*** cosine intermediate child decay angle
356  double cosC[2];
357  //*** delta phi between gluino and child decay planes
358  double dphiGC[2];
359  //*** ratio of child and gluino masses (w/ WIMP masses subtracted)
360  double RCG[2];
361  //*** 1st leading jet pT _associated with this hemisphere_
362  double jet1PT[2];
363  //*** 2nd leading jet pT _associated with this hemisphere_
364  double jet2PT[2];
365  //*** Pinv / HG
366  double Pinv[2];
367 
368  for(int i = 0; i < 2; i++){
369  NV[i] = VIS.GetNElementsInFrame(*VS[i]);
370  NV[i] += VIS.GetNElementsInFrame(*VC[i]);
371 
372  TVector3 vP1 = VS[i]->GetFourVector(*G[i]).Vect();
373  TVector3 vP2 = VC[i]->GetFourVector(*G[i]).Vect();
374  Pinv[i] = 2.*(vP1+vP2).Mag()/(vP1.Mag()+vP2.Mag()+(vP1+vP2).Mag());
375 
376  cosG[i] = G[i]->GetCosDecayAngle();
377 
378  int N = jetID.size();
379  double pTmax[2]; pTmax[0] = -1.; pTmax[1] = -1.;
380  for(int j = 0; j < N; j++){
381  const RestFrame& frame = VIS.GetFrame(jetID[j]);
382  if(VS[i]->IsSame(frame) || VC[i]->IsSame(frame)){
383  double pT = VIS.GetLabFrameFourVector(jetID[j]).Pt();
384  if(pT > pTmax[0]){
385  pTmax[1] = pTmax[0];
386  pTmax[0] = pT;
387  } else {
388  if(pT > pTmax[1]) pTmax[1] = pT;
389  }
390  }
391  }
392 
393  jet1PT[i] = pTmax[0];
394  jet2PT[i] = pTmax[1];
395 
396  if(NV[i] > 1){
397  cosC[i] = C[i]->GetCosDecayAngle();
398  dphiGC[i] = G[i]->GetDeltaPhiDecayPlanes(*C[i]);
399  RCG[i] = (C[i]->GetMass()-X[i]->GetMass())/(G[i]->GetMass()-X[i]->GetMass());
400  } else {
401  cosC[i] = -2.;
402  dphiGC[i] = -1.;
403  RCG[i] = -1.;
404  jet2PT[i] = -1.;
405  }
406  }
407 
408  //
409  // background tree observables
410  //
411  TLorentzVector Psib = I_B.GetSiblingFrame().GetFourVector(LAB_B);
412  TLorentzVector Pmet = I_B.GetFourVector(LAB_B);
413 
414  //***
415  double Rpsib = std::max(0.,Psib.Vect().Dot(Pmet.Vect().Unit()));
416  Rpsib = Rpsib / (Pmet.Pt() + Rpsib);
417 
418  TVector3 boostQCD = (Pmet+Psib).BoostVector();
419  Psib.Boost(-boostQCD);
420  double cosQCD = -1.*Psib.Vect().Unit().Dot(boostQCD.Unit());
421  cosQCD = (1.-cosQCD)/2.;
422 
423  //***
424  double DeltaQCD = (cosQCD-Rpsib)/(cosQCD+Rpsib);
425 
426  }
427 
428  TFile fout(output_name.c_str(),"RECREATE");
429  fout.Close();
430  tree_plot->WriteOutput(output_name);
431 }
432 
433 # ifndef __CINT__ // main function for stand-alone compilation
434 int main(){
435  example_DiGluino_to_bbXbbX();
436  return 0;
437 }
438 #endif
RestFrames::RestFrame
abstract base class for all Frame objects
Generated on Fri Aug 12 2016 14:55:13 for RestFrames by   doxygen 1.8.9.1