17 #include <TLorentzVector.h>
22 #include <gsl/gsl_randist.h>
23 #include <gsl/gsl_rng.h>
46 if (name1.compare(
"e") == 0)
50 else if (name1.compare(
"mu") == 0)
56 std::cout <<
"Invalid decay " << name1 <<
", valid is e or mu" << std::endl;
65 if ((name1.compare(
"e-") == 0 && name2.compare(
"e+") == 0) ||
66 (name1.compare(
"e+") == 0 && name2.compare(
"e-") == 0) ||
67 (name1.compare(
"mu+") == 0 && name2.compare(
"mu-") == 0) ||
68 (name1.compare(
"mu-") == 0 && name2.compare(
"mu+") == 0))
70 decay1_names.insert(std::pair<unsigned int, std::string>(decay_id, name1));
71 decay2_names.insert(std::pair<unsigned int, std::string>(decay_id, name2));
77 std::cout <<
"invalid decay channel Y --> " << name1 <<
" + " << name2 << std::endl;
165 static const double mmuon = 105.6583745e-3;
167 static const double melectron = 0.5109989461e-3;
170 if (
decay1.compare(
"e+") == 0 ||
decay1.compare(
"e-") == 0)
174 else if (
decay1.compare(
"mu+") == 0 ||
decay1.compare(
"mu-") == 0)
180 std::cout <<
"Do not recognize the decay type " <<
decay1 << std::endl;
185 if (
decay2.compare(
"e+") == 0 ||
decay2.compare(
"e-") == 0)
189 else if (
decay2.compare(
"mu+") == 0 ||
decay2.compare(
"mu-") == 0)
195 std::cout <<
"Do not recognize the decay type " <<
decay2 << std::endl;
206 std::cout <<
"PHG4ParticleGeneratorVectorMeson::InitRun started." << std::endl;
208 trand =
new TRandom3();
210 trand->SetSeed(iseed);
214 gRandom->SetSeed(iseed);
217 fsin =
new TF1(
"fsin",
"sin(x)", 0,
M_PI);
221 frap->SetParameter(0, 1.0);
222 frap->SetParameter(1, 0.0);
223 frap->SetParameter(2, 0.8749);
226 fpt =
new TF1(
"fpt",
"2.0*3.14159*x*[0]*pow((1 + x*x/(4*[1]) ),-[2])",
pt_min,
pt_max);
227 fpt->SetParameter(0, 72.1);
228 fpt->SetParameter(1, 26.516);
229 fpt->SetParameter(2, 10.6834);
231 ineve = findNode::getClass<PHG4InEvent>(topNode,
"PHG4INEVENT");
245 std::cout <<
"PHG4ParticleGeneratorVectorMeson::InitRun endeded." << std::endl;
252 if (!
ineve) std::cout <<
" G4InEvent not found " << std::endl;
262 pt =
fpt->GetRandom();
273 y =
frap->GetRandom();
288 double mt = sqrt(mnow * mnow + pt * pt);
289 double eta = asinh(sinh(y) * mt / pt);
294 vm.SetPtEtaPhiM(pt, eta, phi, mnow);
313 unsigned int decay_id =
it->first;
314 std::string decay1_name =
it->second;
315 std::string decay2_name;
316 std::map<unsigned int, std::string>::iterator jt =
decay2_names.find(decay_id);
323 decay2_name = jt->second;
329 std::cout <<
PHWHERE <<
"Decay particles && vertex info can't be found !!" << std::endl;
349 else if (decay_id == 0)
357 double E1 = (mnow * mnow -
m2 *
m2 +
m1 *
m1) / (2.0 * mnow);
358 double p1 = sqrt((mnow * mnow - (
m1 +
m2) * (
m1 +
m2)) * (mnow * mnow - (
m1 -
m2) * (
m1 -
m2))) / (2.0 * mnow);
364 double th1 =
fsin->GetRandom();
370 double px1 = p1 * sin(th1) * cos(phi1);
371 double py1 = p1 * sin(th1) * sin(phi1);
372 double pz1 = p1 * cos(th1);
374 v1.SetPxPyPzE(px1, py1, pz1, E1);
379 double betax = vm.Px() / vm.E();
380 double betay = vm.Py() / vm.E();
381 double betaz = vm.Pz() / vm.E();
382 v1.Boost(betax, betay, betaz);
386 TLorentzVector
v2 = vm -
v1;
390 AddParticle(decay1_name, v1.Px(), v1.Py(), v1.Pz());
391 AddParticle(decay2_name, v2.Px(), v2.Py(), v2.Pz());
412 std::cout << std::endl
413 <<
"Output some sanity check info from PHG4ParticleGeneratorVectorMeson:" << std::endl;
415 std::cout <<
" Vertex for this event (X,Y,Z) is (" <<
get_vtx_x() <<
", " <<
get_vtx_y() <<
", " <<
get_vtx_z() <<
")" << std::endl;
418 std::cout <<
" Decay particle 1:"
422 <<
" eta " << v1.PseudoRapidity()
423 <<
" phi " << v1.Phi()
424 <<
" theta " << v1.Theta()
426 <<
" mass " << v1.M()
430 std::cout <<
" Decay particle 2:"
434 <<
" eta " << v2.PseudoRapidity()
435 <<
" phi " << v2.Phi()
436 <<
" theta " << v2.Theta()
438 <<
" mass " << v2.M()
443 std::cout <<
" Vector meson input kinematics: mass " << vm.M()
447 <<
" eta " << vm.PseudoRapidity()
448 <<
" y " << vm.Rapidity()
455 TLorentzVector vreco = v1 +
v2;
457 std::cout <<
" Reco'd vector meson kinematics: mass " << vreco.M()
458 <<
" px " << vreco.Px()
459 <<
" py " << vreco.Py()
460 <<
" pz " << vreco.Pz()
461 <<
" eta " << vreco.PseudoRapidity()
462 <<
" y " << vreco.Rapidity()
463 <<
" pt " << vreco.Pt()
464 <<
" E " << vreco.E()
482 else if (dist ==
Gaus)