11 #include <TLorentzVector.h>
15 #include <gsl/gsl_const.h>
16 #include <gsl/gsl_randist.h>
17 #include <gsl/gsl_rng.h>
76 std::cout <<
Name() <<
" random seed: " << iseed << std::endl;
77 gRandom->SetSeed(iseed);
79 fsin =
new TF1(
"fsin",
"sin(x)", 0,
M_PI);
83 frap->SetParameter(0, 1.0);
84 frap->SetParameter(1, 0.0);
85 frap->SetParameter(2, 0.8749);
88 fpt =
new TF1(
"fpt",
"[0]*x*pow((1.0 + x*x/[1]/[1]),[2])",
pt_min,
pt_max);
89 fpt->SetParameter(0, 1.16930
e+04);
90 fpt->SetParameter(1, 3.03486
e+00);
91 fpt->SetParameter(2, -5.42819
e+00);
98 PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode,
"PHG4INEVENT");
100 double tau = 410.1e-15;
101 double cc = GSL_CONST_CGSM_SPEED_OF_LIGHT;
102 double ctau = tau * cc * 1.0e+04;
121 pt =
fpt->GetRandom();
133 y =
frap->GetRandom();
146 double mt = sqrt(mnow * mnow + pt * pt);
147 double eta = asinh(sinh(y) * mt / pt);
152 vd0.SetPtEtaPhiM(pt, eta, phi, mnow);
154 double beta = vd0.Beta();
155 double gamma = vd0.Gamma();
156 double lifetime = gsl_ran_exponential(
RandomGenerator(), 410.1
e-03) * 1.0e-12;
157 double lifepath = lifetime * gamma * beta * cc;
160 std::cout <<
"D0 px,py,pz: " << vd0.Px() <<
" " << vd0.Py() <<
" " << vd0.Pz() <<
" " << beta <<
" " << gamma << std::endl;
161 std::cout <<
" ctau = " << ctau <<
" " << lifetime <<
" " << lifepath <<
" " << lifepath * 1.0e+04 << std::endl;
163 set_vtx(vd0.Px() / vd0.P() * lifepath,
164 vd0.Py() / vd0.P() * lifepath,
165 get_vtx_z() + vd0.Pz() / vd0.P() * lifepath);
176 double E1 = (mnow * mnow -
m2 *
m2 +
m1 *
m1) / (2.0 * mnow);
177 double p1 = sqrt((mnow * mnow - (
m1 +
m2) * (
m1 +
m2)) * (mnow * mnow - (
m1 -
m2) * (
m1 -
m2))) / (2.0 * mnow);
183 double th1 =
fsin->GetRandom();
188 double px1 = p1 * sin(th1) * cos(phi1);
189 double py1 = p1 * sin(th1) * sin(phi1);
190 double pz1 = p1 * cos(th1);
192 v1.SetPxPyPzE(px1, py1, pz1, E1);
197 double betax = vd0.Px() / vd0.E();
198 double betay = vd0.Py() / vd0.E();
199 double betaz = vd0.Pz() / vd0.E();
200 v1.Boost(betax, betay, betaz);
204 TLorentzVector
v2 = vd0 -
v1;
231 std::cout << std::endl
232 <<
"Output some sanity check info from PHG4ParticleGeneratorD0:" << std::endl;
236 std::cout <<
"Kaon : " << v1.Pt() <<
" " << v1.PseudoRapidity() <<
" " << v1.M() << std::endl;
237 std::cout <<
"pion : " << v2.Pt() <<
" " << v2.PseudoRapidity() <<
" " << v2.M() << std::endl;
238 std::cout <<
"D0 : " << vd0.Pt() <<
" " << vd0.PseudoRapidity() <<
" " << vd0.M() << std::endl;
239 TLorentzVector vreco = v1 +
v2;
240 std::cout <<
"reconstructed D0 : " << vreco.Pt() <<
" " << vreco.PseudoRapidity() <<
" " << vreco.M() << std::endl;