30 #if (!defined(__CINT__) && !defined(__CLING__))
33 #if defined(__MAKECINT__) || defined(__ROOTCLING__) || defined(COMPILER)
39 using namespace RestFrames;
41 void example_DiStop_to_bXp_bXm_to_blNblN(std::string output_name =
42 "output_DiStop_to_bXp_bXm_to_blNblN.root"){
58 g_Log << LogInfo <<
"Initializing generator frames and tree..." << LogEnd;
75 LAB_Gen.SetChildFrame(CM_Gen);
76 CM_Gen.AddChildFrame(Ta_Gen);
77 CM_Gen.AddChildFrame(Tb_Gen);
78 Ta_Gen.AddChildFrame(Ba_Gen);
79 Ta_Gen.AddChildFrame(Xa_Gen);
80 Tb_Gen.AddChildFrame(Bb_Gen);
81 Tb_Gen.AddChildFrame(Xb_Gen);
82 Xa_Gen.AddChildFrame(La_Gen);
83 Xa_Gen.AddChildFrame(Na_Gen);
84 Xb_Gen.AddChildFrame(Lb_Gen);
85 Xb_Gen.AddChildFrame(Nb_Gen);
87 if(LAB_Gen.InitializeTree())
88 g_Log << LogInfo <<
"...Successfully initialized generator tree" << LogEnd;
90 g_Log << LogError <<
"...Failed initializing generator tree" << LogEnd;
97 CM_Gen.SetVariableMass();
99 Ta_Gen.SetMass(mT); Tb_Gen.SetMass(mT);
101 Xa_Gen.SetMass(0.5*(mT+mN)); Xb_Gen.SetMass(0.5*(mT+mN));
103 Ba_Gen.SetMass(mB); Bb_Gen.SetMass(mB);
105 La_Gen.SetMass(mL); Lb_Gen.SetMass(mL);
107 Na_Gen.SetMass(mN); Nb_Gen.SetMass(mN);
110 Ba_Gen.SetPtCut(20.); Bb_Gen.SetPtCut(20.);
111 Ba_Gen.SetEtaCut(2.5); Bb_Gen.SetEtaCut(2.5);
112 La_Gen.SetPtCut(15.); Lb_Gen.SetPtCut(15.);
113 La_Gen.SetEtaCut(2.5); Lb_Gen.SetEtaCut(2.5);
115 if(LAB_Gen.InitializeAnalysis())
116 g_Log << LogInfo <<
"...Successfully initialized generator analysis" << LogEnd;
118 g_Log << LogError <<
"...Failed initializing generator analysis" << LogEnd;
122 g_Log << LogInfo <<
"Initializing reconstruction frames and tree..." << LogEnd;
140 LAB.SetChildFrame(CM);
141 CM.AddChildFrame(Ta);
142 CM.AddChildFrame(Tb);
143 Ta.AddChildFrame(Ba);
144 Ta.AddChildFrame(Xa);
145 Tb.AddChildFrame(Bb);
146 Tb.AddChildFrame(Xb);
147 Xa.AddChildFrame(La);
148 Xa.AddChildFrame(Na);
149 Xb.AddChildFrame(Lb);
150 Xb.AddChildFrame(Nb);
152 if(LAB.InitializeTree())
153 g_Log << LogInfo <<
"...Successfully initialized reconstruction tree" << LogEnd;
155 g_Log << LogError <<
"...Failed initializing reconstruction tree" << LogEnd;
160 std::string group_name;
163 group_name =
"#tilde{#nu} #tilde{#nu} Jigsaws";
170 B.SetNElementsForFrame(Ba, 1);
171 B.SetNElementsForFrame(Bb, 1);
174 std::string jigsaw_name;
177 jigsaw_name =
"M_{#tilde{#nu} #tilde{#nu}} ~ m_{#it{l} #it{l}}";
179 INV.AddJigsaw(NuNuM);
181 jigsaw_name =
"#eta_{#tilde{#nu} #tilde{#nu}} = #eta_{b #it{l} b #it{l}}";
183 INV.AddJigsaw(NuNuR);
184 NuNuR.AddVisibleFrames(CM.GetListVisibleFrames());
186 jigsaw_name =
"min ( M_{#tilde{t} a}- M_{#tilde{t} b} )^{2}";
189 INV.AddJigsaw(MinDeltaMt);
190 MinDeltaMt.AddInvisibleFrame(Na, 0);
191 MinDeltaMt.AddInvisibleFrame(Nb, 1);
192 MinDeltaMt.AddVisibleFrames(La+Ba, 0);
193 MinDeltaMt.AddVisibleFrames(Lb+Bb, 1);
194 MinDeltaMt.AddMassFrame(La, 0);
195 MinDeltaMt.AddMassFrame(Lb, 1);
198 jigsaw_name =
"Minimize M(b #it{l} )_{a} , M(b #it{l} )_{b}";
202 MinBL.AddFrames(La+Ba,0);
203 MinBL.AddFrames(Lb+Bb,1);
205 if(LAB.InitializeAnalysis())
206 g_Log << LogInfo <<
"...Successfully initialized analysis" << LogEnd;
208 g_Log << LogError <<
"...Failed initializing analysis" << LogEnd;
216 treePlot->
Draw(
"GenTree",
"Generator Tree",
true);
219 treePlot->
Draw(
"RecoTree",
"Reconstruction Tree");
222 treePlot->
Draw(
"VisTree",
"b-jet Jigsaws",
true);
225 treePlot->
Draw(
"InvTree",
"Inivisible Jigsaws",
true);
229 std::string plot_title =
230 "#tilde{t} #tilde{t} #rightarrow 2x [ b #tilde{#chi}^{ #pm}(#it{l} #tilde{#nu}) ]";
232 sprintf(smasses,
"m_{#tilde{t}}= %.0f, m_{#tilde{#nu}}= %.0f", mT, mN);
233 HistPlot* histPlot =
new HistPlot(
"HistPlot", plot_title+
" , "+std::string(smasses));
235 std::string R_MXN =
"#frac{m_{#tilde{#chi}^{ #pm}} - m_{#tilde{#nu}}}";
236 std::string R_MXD =
"{m_{#tilde{t}} - m_{#tilde{#nu}}}";
238 for(
int m = 0; m < NMX; m++){
239 char smass[200], scat[50];
240 sprintf(scat,
"MX%d_%d", m+1, NMX+1);
241 sprintf(smass,
" = #frac{%d}{%d}", m+1, NMX+1);
243 cat_list += histPlot->
GetNewCategory(scat,
"#scale[0.7]{"+R_MXN+R_MXD+std::string(smass)+
"}");
248 plot_title+
" , "+std::string(smasses));
252 "E_{b a}^{#tilde{t} a} / E_{b a}^{#tilde{t} a gen}", 0., 2.);
254 "E_{b b}^{#tilde{t} b} / E_{b b}^{#tilde{t} b gen}", 0., 2.);
256 "E_{#it{l} a}^{#tilde{#chi}^{ #pm} a} / E_{#it{l} a}^{#tilde{#chi}^{ #pm} a true}", 0., 2.);
258 "E_{#it{l} b}^{#tilde{#chi}^{ #pm} b} / E_{#it{l} b}^{#tilde{#chi}^{ #pm} b true}", 0., 2.);
262 const HistPlotVar& cosXa = histPlot->
GetNewVar(
"cosXa",
"cos #theta_{#tilde{#chi}^{ #pm} a}", -1., 1.);
263 const HistPlotVar& cosXb = histPlot->
GetNewVar(
"cosXb",
"cos #theta_{#tilde{#chi}^{ #pm} b}", -1., 1.);
265 -acos(-1.)/2., acos(-1.)/2.);
267 "#theta_{#tilde{t} a} - #theta_{#tilde{t} a}^{true}", -acos(-1.)/2., acos(-1.)/2.);
269 "#theta_{#tilde{t} b} - #theta_{#tilde{t} b}^{true}", -acos(-1.)/2., acos(-1.)/2.);
271 "#theta_{#tilde{#chi}^{ #pm} a} - #theta_{#tilde{#chi}^{ #pm} a}^{true}", -acos(-1.)/2., acos(-1.)/2.);
273 "#theta_{#tilde{#chi}^{ #pm} b} - #theta_{#tilde{#chi}^{ #pm} b}^{true}", -acos(-1.)/2., acos(-1.)/2.);
276 histPlot->
AddPlot(MCM, cat_list);
277 histPlot->
AddPlot(Eb_Ta, cat_list);
278 histPlot->
AddPlot(El_Xa, cat_list);
279 histPlot->
AddPlot(DcosCM, cat_list);
280 histPlot->
AddPlot(DcosTa, cat_list);
281 histPlot->
AddPlot(DcosXa, cat_list);
283 histPlot->
AddPlot(MCM, Eb_Ta, cat_list[NMX/2]);
284 histPlot->
AddPlot(MCM, El_Xa, cat_list[NMX/2]);
285 histPlot->
AddPlot(Eb_Ta, Eb_Tb, cat_list[NMX/2]);
286 histPlot->
AddPlot(El_Xa, El_Xb, cat_list[NMX/2]);
287 histPlot->
AddPlot(Eb_Ta, El_Xa, cat_list[NMX/2]);
288 histPlot->
AddPlot(DcosCM, MCM, cat_list[NMX/2]);
289 histPlot->
AddPlot(DcosTa, Eb_Ta, cat_list[NMX/2]);
290 histPlot->
AddPlot(DcosXa, El_Xa, cat_list[NMX/2]);
291 histPlot->
AddPlot(DcosCM, DcosTa, cat_list[NMX/2]);
292 histPlot->
AddPlot(DcosTa, DcosTb, cat_list[NMX/2]);
293 histPlot->
AddPlot(DcosXa, DcosXb, cat_list[NMX/2]);
296 histPlot->
AddPlot(MCM, RMX, cat_Reco);
297 histPlot->
AddPlot(Eb_Ta, RMX, cat_Reco);
298 histPlot->
AddPlot(El_Xa, RMX, cat_Reco);
304 for(
int m = 0; m < NMX; m++){
305 g_Log << LogInfo <<
"Generating events for ";
306 g_Log <<
"(mXp - mN) / (mT - mN)";
307 g_Log <<
" = " << m+1 <<
" / " << NMX+1 << LogEnd;
309 RMX = double(m+1)/double(NMX+1);
310 Xa_Gen.SetMass(RMX*mT + (1.-RMX)*mN);
315 Xb_Gen.SetMass(RMX*mT + (1.-RMX)*mN);
318 LAB_Gen.InitializeAnalysis();
320 for(
int igen = 0; igen < Ngen; igen++){
321 if(igen%((std::max(Ngen,10))/10) == 0)
322 g_Log << LogInfo <<
"Generating event " << igen <<
" of " << Ngen << LogEnd;
325 LAB_Gen.ClearEvent();
327 LAB_Gen.AnalyzeEvent();
329 TVector3 MET = LAB_Gen.GetInvisibleMomentum();
335 INV.SetLabFrameThreeVector(MET);
337 La.SetLabFrameFourVector(La_Gen.GetFourVector());
338 Lb.SetLabFrameFourVector(Lb_Gen.GetFourVector());
340 B.AddLabFrameFourVector(Ba_Gen.GetFourVector());
341 B.AddLabFrameFourVector(Bb_Gen.GetFourVector());
349 double MCMgen = CM_Gen.GetMass();
350 double Eb_Tagen = Ba_Gen.GetFourVector(Ta_Gen).E();
351 double Eb_Tbgen = Bb_Gen.GetFourVector(Tb_Gen).E();
352 double El_Xagen = La_Gen.GetFourVector(Xa_Gen).E();
353 double El_Xbgen = Lb_Gen.GetFourVector(Xb_Gen).E();
354 double cosCMgen = CM_Gen.GetCosDecayAngle();
355 double cosTagen = Ta_Gen.GetCosDecayAngle();
356 double cosTbgen = Tb_Gen.GetCosDecayAngle();
357 double cosXagen = Xa_Gen.GetCosDecayAngle();
358 double cosXbgen = Xb_Gen.GetCosDecayAngle();
360 MCM = CM.GetMass()/MCMgen;
361 Eb_Ta = Ba.GetFourVector(Ta).E()/Eb_Tagen;
362 Eb_Tb = Bb.GetFourVector(Tb).E()/Eb_Tbgen;
363 El_Xa = La.GetFourVector(Xa).E()/El_Xagen;
364 El_Xb = Lb.GetFourVector(Xb).E()/El_Xbgen;
365 cosCM = CM.GetCosDecayAngle();
366 cosTa = Ta.GetCosDecayAngle();
367 cosTb = Tb.GetCosDecayAngle();
368 cosXa = Xa.GetCosDecayAngle();
369 cosXb = Xb.GetCosDecayAngle();
370 DcosCM = asin(sqrt(1.-cosCM*cosCM)*cosCMgen-sqrt(1.-cosCMgen*cosCMgen)*cosCM);
371 DcosTa = asin(sqrt(1.-cosTa*cosTa)*cosTagen-sqrt(1.-cosTagen*cosTagen)*cosTa);
372 DcosTb = asin(sqrt(1.-cosTb*cosTb)*cosTbgen-sqrt(1.-cosTbgen*cosTbgen)*cosTb);
373 DcosXa = asin(sqrt(1.-cosXa*cosXa)*cosXagen-sqrt(1.-cosXagen*cosXagen)*cosXa);
374 DcosXb = asin(sqrt(1.-cosXb*cosXb)*cosXbgen-sqrt(1.-cosXbgen*cosXbgen)*cosXb);
376 histPlot->
Fill(cat_list[m]);
377 histPlot->
Fill(cat_Reco);
380 LAB_Gen.PrintGeneratorEfficiency();
385 TFile fout(output_name.c_str(),
"RECREATE");
391 g_Log << LogInfo <<
"Finished" << LogEnd;
394 # ifndef __CINT__ // main function for stand-alone compilation
396 example_DiStop_to_bXp_bXm_to_blNblN();