30 #if (!defined(__CINT__) && !defined(__CLING__))
33 #if defined(__MAKECINT__) || defined(__ROOTCLING__) || defined(COMPILER)
39 using namespace RestFrames;
43 void example_X2X2_to_ZllXHggX(std::string output_name =
44 "output_X2X2_to_ZllXHggX.root"){
55 std::vector<double> mX1;
62 g_Log << LogInfo <<
"Initializing generator frames and tree..." << LogEnd;
65 DecayGenFrame X2X2_Gen(
"X2X2_Gen",
"#tilde{#chi}^{ 0}_{2} #tilde{#chi}^{ 0}_{2}");
79 LAB_Gen.SetChildFrame(X2X2_Gen);
80 X2X2_Gen.AddChildFrame(X2a_Gen);
81 X2X2_Gen.AddChildFrame(X2b_Gen);
82 X2a_Gen.AddChildFrame(Za_Gen);
83 X2a_Gen.AddChildFrame(X1a_Gen);
84 X2b_Gen.AddChildFrame(Hb_Gen);
85 X2b_Gen.AddChildFrame(X1b_Gen);
86 Za_Gen.AddChildFrame(L1_Gen);
87 Za_Gen.AddChildFrame(L2_Gen);
88 Hb_Gen.AddChildFrame(G1_Gen);
89 Hb_Gen.AddChildFrame(G2_Gen);
91 if(LAB_Gen.InitializeTree())
92 g_Log << LogInfo <<
"...Successfully initialized generator tree" << LogEnd;
94 g_Log << LogError <<
"...Failed initializing generator tree" << LogEnd;
98 if(NmX1 < 1 || mX1.size() < 1)
103 X2X2_Gen.SetVariableMass();
104 X2a_Gen.SetMass(mX2); X2b_Gen.SetMass(mX2);
105 X1a_Gen.SetMass(mX1[0]); X1b_Gen.SetMass(mX1[0]);
111 L1_Gen.SetPtCut(8.); L1_Gen.SetEtaCut(2.5);
112 L2_Gen.SetPtCut(8.); L2_Gen.SetEtaCut(2.5);
113 G1_Gen.SetPtCut(20.); G1_Gen.SetEtaCut(3.0);
114 G2_Gen.SetPtCut(20.); G2_Gen.SetEtaCut(3.0);
116 if(LAB_Gen.InitializeAnalysis())
117 g_Log << LogInfo <<
"...Successfully initialized generator analysis" << LogEnd;
119 g_Log << LogError <<
"...Failed initializing generator analysis" << LogEnd;
124 g_Log << LogInfo <<
"Initializing reconstruction frames and tree..." << LogEnd;
127 DecayRecoFrame X2X2(
"X2X2",
"#tilde{#chi}^{ 0}_{2} #tilde{#chi}^{ 0}_{2}");
141 LAB.SetChildFrame(X2X2);
142 X2X2.AddChildFrame(X2a);
143 X2X2.AddChildFrame(X2b);
144 X2a.AddChildFrame(Za);
145 X2a.AddChildFrame(X1a);
146 X2b.AddChildFrame(Hb);
147 X2b.AddChildFrame(X1b);
148 Za.AddChildFrame(L1);
149 Za.AddChildFrame(L2);
150 Hb.AddChildFrame(G1);
151 Hb.AddChildFrame(G2);
153 if(LAB.InitializeTree())
154 g_Log << LogInfo <<
"...Successfully initialized reconstruction trees" << LogEnd;
156 g_Log << LogError <<
"...Failed initializing reconstruction trees" << LogEnd;
162 INV.AddFrames(X1a+X1b);
170 SetRapidityInvJigsaw X1_eta(
"X1_eta",
"#eta_{#tilde{#chi}_{1}^{ 0} #tilde{#chi}_{1}^{ 0}} = #eta_{2#gamma+2#it{l}}");
171 INV.AddJigsaw(X1_eta);
172 X1_eta.AddVisibleFrames(X2X2.GetListVisibleFrames());
177 INV.AddJigsaw(X1X1_contra);
178 X1X1_contra.AddVisibleFrames(X2a.GetListVisibleFrames(), 0);
179 X1X1_contra.AddVisibleFrames(X2b.GetListVisibleFrames(), 1);
180 X1X1_contra.AddInvisibleFrame(X1a, 0);
181 X1X1_contra.AddInvisibleFrame(X1b, 1);
183 if(LAB.InitializeAnalysis())
184 g_Log << LogInfo <<
"...Successfully initialized analysis" << std::endl << LogEnd;
186 g_Log << LogError <<
"...Failed initializing analysis" << LogEnd;
193 treePlot->
Draw(
"GenTree",
"Generator Tree",
true);
196 treePlot->
Draw(
"RecoTree",
"Reconstruction Tree");
199 treePlot->
Draw(
"InvTree",
"Inivisible Jigsaws",
true);
205 std::string(
"#tilde{#chi}_{2}^{ 0} #tilde{#chi}_{2}^{ 0}") +
206 "#rightarrow Z(#it{l}#it{l}) #tilde{#chi}_{1}^{ 0}"+
207 "H(#gamma #gamma) #tilde{#chi}_{1}^{ 0}");
213 sprintf(smassX2,
"m_{#tilde{#chi}^{0}_{2}}= %.0f", mX2);
214 for(
int m = 0; m < NmX1; m++){
216 mX1.push_back(mX1[0] + m*(mX2-mX1[0])/
double(NmX1));
218 char smassX1[200], scat[50];
219 sprintf(scat,
"MX1_%.0f", mX1[m]);
220 sprintf(smassX1,
"m_{#tilde{#chi}^{0}_{1}}= %.0f", mX1[m]);
223 " , "+std::string(smassX1));
227 "M_{#tilde{#chi}_{2}^{ 0} #tilde{#chi}_{2}^{ 0}}",
229 const HistPlotVar& EZX2a = histPlot->
GetNewVar(
"EZX2a",
"E_{Z}^{ #tilde{#chi}_{2 a}^{ 0}}", 0., 2.);
230 const HistPlotVar& EHX2b = histPlot->
GetNewVar(
"EHX2b",
"E_{H}^{ #tilde{#chi}_{2 b}^{ 0}}", 0., 2.);
231 const HistPlotVar& cosX2a = histPlot->
GetNewVar(
"cosX2a",
"cos #theta_{#tilde{#chi}_{2 a}^{ 0}}", -1., 1.);
232 const HistPlotVar& cosX2b = histPlot->
GetNewVar(
"cosX2b",
"cos #theta_{#tilde{#chi}_{2 b}^{ 0}}", -1., 1.);
235 const HistPlotVar& DcosZ = histPlot->
GetNewVar(
"DcosZ",
"#theta_{Z} - #theta_{Z}^{gen}", -1., 1.);
236 const HistPlotVar& DcosH = histPlot->
GetNewVar(
"DcosH",
"#theta_{H} - #theta_{H}^{gen}", -1., 1.);
237 const HistPlotVar& DcosX2a = histPlot->
GetNewVar(
"DcosX2a",
"#theta_{X2a} - #theta_{X2a}^{gen}", -1., 1.);
238 const HistPlotVar& DcosX2b = histPlot->
GetNewVar(
"DcosX2b",
"#theta_{X2b} - #theta_{X2b}^{gen}", -1., 1.);
266 histPlot->
AddPlot(RISR, cat_list);
267 histPlot->
AddPlot(RISR, PTISR, cat_list[NmX1/2]);
270 histPlot->
AddPlot(MI, cat_list);
271 histPlot->
AddPlot(MP, cat_list);
272 histPlot->
AddPlot(RISR2, cat_list);
273 histPlot->
AddPlot(MI, PTISR, cat_list);
274 histPlot->
AddPlot(MP, PTISR, cat_list);
275 histPlot->
AddPlot(MI, RISR, cat_list);
276 histPlot->
AddPlot(MP, RISR, cat_list);
277 histPlot->
AddPlot(MI, RISR2, cat_list);
278 histPlot->
AddPlot(MP, RISR2, cat_list);
279 histPlot->
AddPlot(MP, MI, cat_list);
280 histPlot->
AddPlot(RISR2, PTISR, cat_list[NmX1/2]);
281 histPlot->
AddPlot(RISR2, RISR, cat_list[NmX1/2]);
283 for(
int m = 0; m < NmX1; m++){
284 g_Log << LogInfo <<
"Generating events for ";
285 g_Log <<
"mX2 = " << mX2 <<
" , ";
286 g_Log <<
"mX1 = " << mX1[m] << LogEnd;
288 X1a_Gen.SetMass(mX1[m]);
289 X1b_Gen.SetMass(mX1[m]);
291 LAB_Gen.InitializeAnalysis();
293 for(
int igen = 0; igen < Ngen; igen++){
294 if(igen%((std::max(Ngen,10))/10) == 0)
295 g_Log << LogInfo <<
"Generating event " << igen <<
" of " << Ngen << LogEnd;
298 LAB_Gen.ClearEvent();
301 LAB_Gen.SetTransverseMomentum(800.*gRandom->Rndm());
303 LAB_Gen.AnalyzeEvent();
306 TLorentzVector P_V = X2X2_Gen.GetListVisibleFrames().GetFourVector();
307 TLorentzVector P_I = X2X2_Gen.GetListInvisibleFrames().GetFourVector();
316 P_V.SetPxPyPzE(P_V.Px(),P_V.Py(),0.,sqrt(P_V.Pt()*P_V.Pt() + P_V.M2()));
317 P_I.SetPxPyPzE(P_I.Px(),P_I.Py(),0.,sqrt(P_I.Pt()*P_I.Pt() + P_I.M2()));
334 double a = ((P_V+P_I).P()*(P_V+P_I).P()*(P_V+P_I).P()*(P_V+P_I).P()-P_I.Vect().Dot((P_V+P_I).Vect())*P_I.Vect().Dot((P_V+P_I).Vect()));
335 double b = -P_V.E()*P_I.Vect().Dot((P_V+P_I).Vect())*P_I.Vect().Dot((P_V+P_I).Vect());
336 double c1 = - P_V.E()*P_V.E()*P_I.Vect().Dot((P_V+P_I).Vect())*P_I.Vect().Dot((P_V+P_I).Vect());
337 double c2 = (4.-1.)*(P_V.Vect().Cross(P_I.Vect()).Mag2())*(P_V+P_I).P()*(P_V+P_I).P();
338 double Ei = (-b + sqrt(b*b - a*(c1+c2))) / a;
339 if(b*b - a*(c1+c2) < 0.)
340 cout <<
" wtf " << endl;
342 cout <<
"HERE" << endl;
349 MI2 = (Ei*Ei - P_I.P()*P_I.P() - P_V.M2() +
350 pow(X2a_Gen.GetListVisibleFrames().GetFourVector().M()+X2b_Gen.GetListVisibleFrames().GetFourVector().M(),2.)) / 4.;
356 double Minv = sqrt(Ei*Ei - P_I.P()*P_I.P());
360 L1.SetLabFrameFourVector(L1_Gen.GetFourVector());
361 L2.SetLabFrameFourVector(L2_Gen.GetFourVector());
362 G1.SetLabFrameFourVector(G1_Gen.GetFourVector());
363 G2.SetLabFrameFourVector(G2_Gen.GetFourVector());
365 TVector3 vMET = LAB_Gen.GetInvisibleMomentum();
369 MET.SetVectM(vMET, Minv);
370 INV.SetLabFrameFourVector(MET);
372 X1a.SetMinimumMass(MI);
373 X1b.SetMinimumMass(MI);
378 MCM = X2X2.GetMass() / X2X2_Gen.GetMass();
380 EZX2a = Za.GetEnergy(X2a) / Za_Gen.GetEnergy(X2a_Gen);
381 EHX2b = Hb.GetEnergy(X2b) / Hb_Gen.GetEnergy(X2b_Gen);
383 cosX2a = X2a.GetCosDecayAngle();
384 double cosX2agen = X2a_Gen.GetCosDecayAngle();
385 cosX2b = X2b.GetCosDecayAngle();
386 double cosX2bgen = X2b_Gen.GetCosDecayAngle();
387 cosZ = Za.GetCosDecayAngle();
388 double cosZgen = Za_Gen.GetCosDecayAngle();
389 cosH = Hb.GetCosDecayAngle();
390 double cosHgen = Hb_Gen.GetCosDecayAngle();
392 DcosX2a = asin(sqrt(1.-cosX2a*cosX2a)*cosX2agen-sqrt(1.-cosX2agen*cosX2agen)*cosX2a);
393 DcosX2b = asin(sqrt(1.-cosX2b*cosX2b)*cosX2bgen-sqrt(1.-cosX2bgen*cosX2bgen)*cosX2b);
394 DcosZ = asin(sqrt(1.-cosZ*cosZ)*cosZgen-sqrt(1.-cosZgen*cosZgen)*cosZ);
395 DcosH = asin(sqrt(1.-cosH*cosH)*cosHgen-sqrt(1.-cosHgen*cosHgen)*cosH);
397 TVector3 vP_ISR = X2X2.GetFourVector(LAB).Vect();
398 TVector3 vP_I = X2X2.GetListInvisibleFrames().GetFourVector(LAB).Vect();
402 RISR = fabs(vP_I.Dot(vP_ISR.Unit())) / vP_ISR.Mag();
403 PTISR = vP_ISR.Mag();
410 histPlot->
Fill(cat_list[m]);
413 LAB_Gen.PrintGeneratorEfficiency();
418 TFile fout(output_name.c_str(),
"RECREATE");
424 g_Log << LogInfo <<
"Finished" << LogEnd;
427 # ifndef __CINT__ // main function for stand-alone compilation
429 example_X2X2_to_ZllXHggX();