LOGO
RestFrames  v1.0.0
RestFrames HEP Event Analysis Software Library
example_DiStop_to_bXp_bXm_to_blNblN.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_DiStop_to_bXp_bXm_to_blNblN(std::string output_name =
40  "output_DiStop_to_bXp_bXm_to_blNblN.root"){
41 
42  double mT = 750.; // stop mass
43  double mN = 100.; // sneutrino mass
44  double mB = 4.18;
45  double mL = 0.106;
46 
47  int NMX = 8; // number of different chargino masses to evaluate
48 
49  bool fix_MXb = false; // fix the "b" chargino mass while varying "a"?
50  double mXb = 600.;
51 
52  // number of events to generate (per chargino mass)
53  int Ngen = 10000;
54 
56  g_Log << LogInfo << "Initializing generator frames and tree..." << LogEnd;
58  ppLabGenFrame LAB_Gen("LAB_Gen","LAB");
59  DecayGenFrame CM_Gen("CM_Gen","CM");
60  DecayGenFrame Ta_Gen("Ta_Gen","#tilde{t}_{a}");
61  DecayGenFrame Tb_Gen("Tb_Gen","#tilde{t}_{b}");
62  DecayGenFrame Xa_Gen("Xa_Gen","#tilde{#chi}_{a}^{ #pm}");
63  DecayGenFrame Xb_Gen("Xb_Gen","#tilde{#chi}_{b}^{ #mp}");
64  VisibleGenFrame Ba_Gen("Ba_Gen","b_{a}");
65  VisibleGenFrame La_Gen("La_Gen","#it{l}_{a}");
66  InvisibleGenFrame Na_Gen("Na_Gen","#tilde{#nu}_{a}");
67  VisibleGenFrame Bb_Gen("Bb_Gen","b_{b}");
68  VisibleGenFrame Lb_Gen("Lb_Gen","#it{l}_{b}");
69  InvisibleGenFrame Nb_Gen("Nb_Gen","#tilde{#nu}_{b}");
70 
71  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
72 
73  LAB_Gen.SetChildFrame(CM_Gen);
74  CM_Gen.AddChildFrame(Ta_Gen);
75  CM_Gen.AddChildFrame(Tb_Gen);
76  Ta_Gen.AddChildFrame(Ba_Gen);
77  Ta_Gen.AddChildFrame(Xa_Gen);
78  Tb_Gen.AddChildFrame(Bb_Gen);
79  Tb_Gen.AddChildFrame(Xb_Gen);
80  Xa_Gen.AddChildFrame(La_Gen);
81  Xa_Gen.AddChildFrame(Na_Gen);
82  Xb_Gen.AddChildFrame(Lb_Gen);
83  Xb_Gen.AddChildFrame(Nb_Gen);
84 
85  if(LAB_Gen.InitializeTree())
86  g_Log << LogInfo << "...Successfully initialized generator tree" << LogEnd;
87  else
88  g_Log << LogError << "...Failed initializing generator tree" << LogEnd;
89 
90  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
91 
92  if(NMX < 1) return;
93 
94  // non-resonant stop production
95  CM_Gen.SetVariableMass();
96  // set stop masses
97  Ta_Gen.SetMass(mT); Tb_Gen.SetMass(mT);
98  // set chargino masses
99  Xa_Gen.SetMass(0.5*(mT+mN)); Xb_Gen.SetMass(0.5*(mT+mN));
100  // set B masses
101  Ba_Gen.SetMass(mB); Bb_Gen.SetMass(mB);
102  // set : masses
103  La_Gen.SetMass(mL); Lb_Gen.SetMass(mL);
104  // set sneutrino masses
105  Na_Gen.SetMass(mN); Nb_Gen.SetMass(mN);
106 
107  // set b-jet/lepton pT/eta cuts
108  Ba_Gen.SetPtCut(20.); Bb_Gen.SetPtCut(20.);
109  Ba_Gen.SetEtaCut(2.5); Bb_Gen.SetEtaCut(2.5);
110  La_Gen.SetPtCut(15.); Lb_Gen.SetPtCut(15.);
111  La_Gen.SetEtaCut(2.5); Lb_Gen.SetEtaCut(2.5);
112 
113  if(LAB_Gen.InitializeAnalysis())
114  g_Log << LogInfo << "...Successfully initialized generator analysis" << LogEnd;
115  else
116  g_Log << LogError << "...Failed initializing generator analysis" << LogEnd;
120  g_Log << LogInfo << "Initializing reconstruction frames and tree..." << LogEnd;
122 
123  LabRecoFrame LAB("LAB","LAB");
124  DecayRecoFrame CM("CM","CM");
125  DecayRecoFrame Ta("Ta","#tilde{t}_{a}");
126  DecayRecoFrame Tb("Tb","#tilde{t}_{b}");
127  DecayRecoFrame Xa("Xa","#tilde{#chi}_{a}^{ #pm}");
128  DecayRecoFrame Xb("Xb","#tilde{#chi}_{b}^{ #mp}");
129  VisibleRecoFrame Ba("Ba","b_{a}");
130  VisibleRecoFrame La("La","#it{l}_{a}");
131  InvisibleRecoFrame Na("Na","#tilde{#nu}_{a}");
132  VisibleRecoFrame Bb("Bb","b_{b}");
133  VisibleRecoFrame Lb("Lb","#it{l}_{b}");
134  InvisibleRecoFrame Nb("Nb","#tilde{#nu}_{b}");
135 
136  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
137 
138  LAB.SetChildFrame(CM);
139  CM.AddChildFrame(Ta);
140  CM.AddChildFrame(Tb);
141  Ta.AddChildFrame(Ba);
142  Ta.AddChildFrame(Xa);
143  Tb.AddChildFrame(Bb);
144  Tb.AddChildFrame(Xb);
145  Xa.AddChildFrame(La);
146  Xa.AddChildFrame(Na);
147  Xb.AddChildFrame(Lb);
148  Xb.AddChildFrame(Nb);
149 
150  if(LAB.InitializeTree())
151  g_Log << LogInfo << "...Successfully initialized reconstruction tree" << LogEnd;
152  else
153  g_Log << LogError << "...Failed initializing reconstruction tree" << LogEnd;
154 
155  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
156 
158  std::string group_name;
159 
160  // Invisible Group
161  group_name = "#tilde{#nu} #tilde{#nu} Jigsaws";
162  InvisibleGroup INV("INV", group_name);
163  INV.AddFrame(Na);
164  INV.AddFrame(Nb);
165  CombinatoricGroup B("VIS","b-jet Jigsaws");
166  B.AddFrame(Ba);
167  B.AddFrame(Bb);
168  B.SetNElementsForFrame(Ba, 1);
169  B.SetNElementsForFrame(Bb, 1);
170 
172  std::string jigsaw_name;
173 
174  // Minimize difference Mt jigsaws
175  jigsaw_name = "M_{#tilde{#nu} #tilde{#nu}} ~ m_{#it{l} #it{l}}";
176  SetMassInvJigsaw NuNuM("NuNuM", jigsaw_name);
177  INV.AddJigsaw(NuNuM);
178 
179  jigsaw_name = "#eta_{#tilde{#nu} #tilde{#nu}} = #eta_{b #it{l} b #it{l}}";
180  SetRapidityInvJigsaw NuNuR("NuNuR", jigsaw_name);
181  INV.AddJigsaw(NuNuR);
182  NuNuR.AddVisibleFrames(CM.GetListVisibleFrames());
183 
184  jigsaw_name = "min ( M_{#tilde{t} a}- M_{#tilde{t} b} )^{2}";
185  //MinMassDiffInvJigsaw MinDeltaMt("MinDeltaMt", jigsaw_name, 2);
186  ContraBoostInvJigsaw MinDeltaMt("MinDeltaMt", jigsaw_name);
187  INV.AddJigsaw(MinDeltaMt);
188  MinDeltaMt.AddInvisibleFrame(Na, 0);
189  MinDeltaMt.AddInvisibleFrame(Nb, 1);
190  MinDeltaMt.AddVisibleFrames(La+Ba, 0);
191  MinDeltaMt.AddVisibleFrames(Lb+Bb, 1);
192  MinDeltaMt.AddMassFrame(La, 0);
193  MinDeltaMt.AddMassFrame(Lb, 1);
194 
195  // b-jet combinatoric jigsaws for all trees
196  jigsaw_name = "Minimize M(b #it{l} )_{a} , M(b #it{l} )_{b}";
197 
198  MinMassesCombJigsaw MinBL("MinBL", jigsaw_name);
199  B.AddJigsaw(MinBL);
200  MinBL.AddFrames(La+Ba,0);
201  MinBL.AddFrames(Lb+Bb,1);
202 
203  if(LAB.InitializeAnalysis())
204  g_Log << LogInfo << "...Successfully initialized analysis" << LogEnd;
205  else
206  g_Log << LogError << "...Failed initializing analysis" << LogEnd;
207 
210 
211  TreePlot* treePlot = new TreePlot("TreePlot","TreePlot");
212 
213  treePlot->SetTree(LAB_Gen);
214  treePlot->Draw("GenTree", "Generator Tree", true);
215 
216  treePlot->SetTree(LAB);
217  treePlot->Draw("RecoTree", "Reconstruction Tree");
218 
219  treePlot->SetTree(B);
220  treePlot->Draw("VisTree", "b-jet Jigsaws", true);
221 
222  treePlot->SetTree(INV);
223  treePlot->Draw("InvTree", "Inivisible Jigsaws", true);
224 
225  //-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//-//
226 
227  std::string plot_title =
228  "#tilde{t} #tilde{t} #rightarrow 2x [ b #tilde{#chi}^{ #pm}(#it{l} #tilde{#nu}) ]";
229  char smasses[200];
230  sprintf(smasses, "m_{#tilde{t}}= %.0f, m_{#tilde{#nu}}= %.0f", mT, mN);
231  HistPlot* histPlot = new HistPlot("HistPlot", plot_title+" , "+std::string(smasses));
232 
233  std::string R_MXN = "#frac{m_{#tilde{#chi}^{ #pm}} - m_{#tilde{#nu}}}";
234  std::string R_MXD = "{m_{#tilde{t}} - m_{#tilde{#nu}}}";
235  RFList<const HistPlotCategory> cat_list;
236  for(int m = 0; m < NMX; m++){
237  char smass[200], scat[50];
238  sprintf(scat, "MX%d_%d", m+1, NMX+1);
239  sprintf(smass, " = #frac{%d}{%d}", m+1, NMX+1);
240 
241  cat_list += histPlot->GetNewCategory(scat, "#scale[0.7]{"+R_MXN+R_MXD+std::string(smass)+"}");
242  }
243 
244  //const HistPlotCategory& cat_Reco = histPlot->GetNewCategory("Reco", smasses);
245  const HistPlotCategory& cat_Reco = histPlot->GetNewCategory("Reco",
246  plot_title+" , "+std::string(smasses));
247 
248  const HistPlotVar& MCM = histPlot->GetNewVar("MCM", "M_{CM} / m_{CM}^{ true}", 0., 2.);
249  const HistPlotVar& Eb_Ta = histPlot->GetNewVar("Eb_Ta",
250  "E_{b a}^{#tilde{t} a} / E_{b a}^{#tilde{t} a gen}", 0., 2.);
251  const HistPlotVar& Eb_Tb = histPlot->GetNewVar("Eb_Tb",
252  "E_{b b}^{#tilde{t} b} / E_{b b}^{#tilde{t} b gen}", 0., 2.);
253  const HistPlotVar& El_Xa = histPlot->GetNewVar("El_Xa",
254  "E_{#it{l} a}^{#tilde{#chi}^{ #pm} a} / E_{#it{l} a}^{#tilde{#chi}^{ #pm} a true}", 0., 2.);
255  const HistPlotVar& El_Xb = histPlot->GetNewVar("El_Xb",
256  "E_{#it{l} b}^{#tilde{#chi}^{ #pm} b} / E_{#it{l} b}^{#tilde{#chi}^{ #pm} b true}", 0., 2.);
257  const HistPlotVar& cosCM = histPlot->GetNewVar("cosCM","cos #theta_{CM}", -1., 1.);
258  const HistPlotVar& cosTa = histPlot->GetNewVar("cosTa","cos #theta_{#tilde{t} a}", -1., 1.);
259  const HistPlotVar& cosTb = histPlot->GetNewVar("cosTb","cos #theta_{#tilde{t} b}", -1., 1.);
260  const HistPlotVar& cosXa = histPlot->GetNewVar("cosXa","cos #theta_{#tilde{#chi}^{ #pm} a}", -1., 1.);
261  const HistPlotVar& cosXb = histPlot->GetNewVar("cosXb","cos #theta_{#tilde{#chi}^{ #pm} b}", -1., 1.);
262  const HistPlotVar& DcosCM = histPlot->GetNewVar("DcosCM","#theta_{CM} - #theta_{CM}^{true}",
263  -acos(-1.)/2., acos(-1.)/2.);
264  const HistPlotVar& DcosTa = histPlot->GetNewVar("Dcosta",
265  "#theta_{#tilde{t} a} - #theta_{#tilde{t} a}^{true}", -acos(-1.)/2., acos(-1.)/2.);
266  const HistPlotVar& DcosTb = histPlot->GetNewVar("Dcostb",
267  "#theta_{#tilde{t} b} - #theta_{#tilde{t} b}^{true}", -acos(-1.)/2., acos(-1.)/2.);
268  const HistPlotVar& DcosXa = histPlot->GetNewVar("DcosXa",
269  "#theta_{#tilde{#chi}^{ #pm} a} - #theta_{#tilde{#chi}^{ #pm} a}^{true}", -acos(-1.)/2., acos(-1.)/2.);
270  const HistPlotVar& DcosXb = histPlot->GetNewVar("DcosXb",
271  "#theta_{#tilde{#chi}^{ #pm} b} - #theta_{#tilde{#chi}^{ #pm} b}^{true}", -acos(-1.)/2., acos(-1.)/2.);
272  const HistPlotVar& RMX = histPlot->GetNewVar("RMX", "#scale[0.6]{"+R_MXN+R_MXD+"}", 0., 1.);
273 
274  histPlot->AddPlot(MCM, cat_list);
275  histPlot->AddPlot(Eb_Ta, cat_list);
276  histPlot->AddPlot(El_Xa, cat_list);
277  histPlot->AddPlot(DcosCM, cat_list);
278  histPlot->AddPlot(DcosTa, cat_list);
279  histPlot->AddPlot(DcosXa, cat_list);
280 
281  histPlot->AddPlot(MCM, Eb_Ta, cat_list[NMX/2]);
282  histPlot->AddPlot(MCM, El_Xa, cat_list[NMX/2]);
283  histPlot->AddPlot(Eb_Ta, Eb_Tb, cat_list[NMX/2]);
284  histPlot->AddPlot(El_Xa, El_Xb, cat_list[NMX/2]);
285  histPlot->AddPlot(Eb_Ta, El_Xa, cat_list[NMX/2]);
286  histPlot->AddPlot(DcosCM, MCM, cat_list[NMX/2]);
287  histPlot->AddPlot(DcosTa, Eb_Ta, cat_list[NMX/2]);
288  histPlot->AddPlot(DcosXa, El_Xa, cat_list[NMX/2]);
289  histPlot->AddPlot(DcosCM, DcosTa, cat_list[NMX/2]);
290  histPlot->AddPlot(DcosTa, DcosTb, cat_list[NMX/2]);
291  histPlot->AddPlot(DcosXa, DcosXb, cat_list[NMX/2]);
292 
293  if(NMX > 1){
294  histPlot->AddPlot(MCM, RMX, cat_Reco);
295  histPlot->AddPlot(Eb_Ta, RMX, cat_Reco);
296  histPlot->AddPlot(El_Xa, RMX, cat_Reco);
297  }
298 
301 
302  for(int m = 0; m < NMX; m++){
303  g_Log << LogInfo << "Generating events for ";
304  g_Log << "(mXp - mN) / (mT - mN)";
305  g_Log << " = " << m+1 << " / " << NMX+1 << LogEnd;
306 
307  RMX = double(m+1)/double(NMX+1);
308  Xa_Gen.SetMass(RMX*mT + (1.-RMX)*mN);
309 
310  if(fix_MXb){
311  Xb_Gen.SetMass(mXb);
312  } else {
313  Xb_Gen.SetMass(RMX*mT + (1.-RMX)*mN);
314  }
315 
316  LAB_Gen.InitializeAnalysis();
317 
318  for(int igen = 0; igen < Ngen; igen++){
319  if(igen%((std::max(Ngen,10))/10) == 0)
320  g_Log << LogInfo << "Generating event " << igen << " of " << Ngen << LogEnd;
321 
322  // generate event
323  LAB_Gen.ClearEvent(); // clear the gen tree
324 
325  LAB_Gen.AnalyzeEvent(); // generate a new event
326 
327  TVector3 MET = LAB_Gen.GetInvisibleMomentum(); // Get the MET from gen tree
328  MET.SetZ(0.);
329 
330  // analyze event
331  LAB.ClearEvent(); // clear the event
332 
333  INV.SetLabFrameThreeVector(MET); // Set MET in reco tree
334 
335  La.SetLabFrameFourVector(La_Gen.GetFourVector()); // Set lepton 4-vectors
336  Lb.SetLabFrameFourVector(Lb_Gen.GetFourVector());
337 
338  B.AddLabFrameFourVector(Ba_Gen.GetFourVector()); // Set b-jet 4-vectors
339  B.AddLabFrameFourVector(Bb_Gen.GetFourVector());
340 
341  LAB.AnalyzeEvent(); // analyze the event
342 
344  // Observable Calculations
346 
347  double MCMgen = CM_Gen.GetMass();
348  double Eb_Tagen = Ba_Gen.GetFourVector(Ta_Gen).E();
349  double Eb_Tbgen = Bb_Gen.GetFourVector(Tb_Gen).E();
350  double El_Xagen = La_Gen.GetFourVector(Xa_Gen).E();
351  double El_Xbgen = Lb_Gen.GetFourVector(Xb_Gen).E();
352  double cosCMgen = CM_Gen.GetCosDecayAngle();
353  double cosTagen = Ta_Gen.GetCosDecayAngle();
354  double cosTbgen = Tb_Gen.GetCosDecayAngle();
355  double cosXagen = Xa_Gen.GetCosDecayAngle();
356  double cosXbgen = Xb_Gen.GetCosDecayAngle();
357 
358  MCM = CM.GetMass()/MCMgen;
359  Eb_Ta = Ba.GetFourVector(Ta).E()/Eb_Tagen;
360  Eb_Tb = Bb.GetFourVector(Tb).E()/Eb_Tbgen;
361  El_Xa = La.GetFourVector(Xa).E()/El_Xagen;
362  El_Xb = Lb.GetFourVector(Xb).E()/El_Xbgen;
363  cosCM = CM.GetCosDecayAngle();
364  cosTa = Ta.GetCosDecayAngle();
365  cosTb = Tb.GetCosDecayAngle();
366  cosXa = Xa.GetCosDecayAngle();
367  cosXb = Xb.GetCosDecayAngle();
368  DcosCM = asin(sqrt(1.-cosCM*cosCM)*cosCMgen-sqrt(1.-cosCMgen*cosCMgen)*cosCM);
369  DcosTa = asin(sqrt(1.-cosTa*cosTa)*cosTagen-sqrt(1.-cosTagen*cosTagen)*cosTa);
370  DcosTb = asin(sqrt(1.-cosTb*cosTb)*cosTbgen-sqrt(1.-cosTbgen*cosTbgen)*cosTb);
371  DcosXa = asin(sqrt(1.-cosXa*cosXa)*cosXagen-sqrt(1.-cosXagen*cosXagen)*cosXa);
372  DcosXb = asin(sqrt(1.-cosXb*cosXb)*cosXbgen-sqrt(1.-cosXbgen*cosXbgen)*cosXb);
373 
374  histPlot->Fill(cat_list[m]);
375  histPlot->Fill(cat_Reco);
376  }
377 
378  LAB_Gen.PrintGeneratorEfficiency();
379  }
380 
381  histPlot->Draw();
382 
383  TFile fout(output_name.c_str(),"RECREATE");
384  fout.Close();
385  histPlot->WriteOutput(output_name);
386  histPlot->WriteHist(output_name);
387  treePlot->WriteOutput(output_name);
388 
389  g_Log << LogInfo << "Finished" << LogEnd;
390 }
391 
392 # ifndef __CINT__ // main function for stand-alone compilation
393 int main(){
394  example_DiStop_to_bXp_bXm_to_blNblN();
395  return 0;
396 }
397 #endif
Generated on Fri Aug 12 2016 14:55:13 for RestFrames by   doxygen 1.8.9.1