28 #include <HepMC/SimpleVector.h>
30 #include <CLHEP/Units/PhysicalConstants.h>
31 #include <CLHEP/Units/SystemOfUnits.h>
32 #include <CLHEP/Vector/Boost.h>
33 #include <CLHEP/Vector/LorentzRotation.h>
34 #include <CLHEP/Vector/LorentzVector.h>
35 #include <CLHEP/Vector/Rotation.h>
37 #include <gsl/gsl_randist.h>
38 #include <gsl/gsl_rng.h>
49 : _vertex_func_x(Gaus)
50 , _vertex_func_y(Gaus)
51 , _vertex_func_z(Gaus)
52 , _vertex_func_t(Gaus)
62 , _reuse_vertex(
false)
63 , _reuse_vertex_embedding_id(numeric_limits<int>::
min())
64 , _geneventmap(nullptr)
83 cout <<
PHWHERE <<
"DST Node missing doing nothing" << endl;
87 _geneventmap = findNode::getClass<PHHepMCGenEventMap>(dstNode,
"PHHepMCGenEventMap");
92 dstNode->addNode(newmapnode);
106 return &mc_evnet_template;
147 assert((hv_index == 0) or (hv_index == 1));
152 return widthA * widthB / sqrt(widthA * widthA + widthB * widthB);
160 const pair<double, double> bunch_zs(
173 const static CLHEP::Hep3Vector y_axis(0, 1, 0);
176 CLHEP::Hep3Vector beamCenterDiffAxis = (beamA_center - beamB_center);
177 assert(beamCenterDiffAxis.mag() > CLHEP::Hep3Vector::getTolerance());
178 beamCenterDiffAxis = beamCenterDiffAxis / beamCenterDiffAxis.mag();
180 CLHEP::Hep3Vector vec_crossing = beamA_center - 0.5 * (beamA_center - beamB_center);
182 CLHEP::Hep3Vector vec_longitudinal_collision = beamCenterDiffAxis * (bunch_zs.first + bunch_zs.second) / 2.;
183 double ct_collision = 0.5 * (-bunch_zs.first + bunch_zs.second) / beamCenterDiffAxis.dot(beamA_center);
185 CLHEP::Hep3Vector vec_crossing_collision = ct_collision * vec_crossing;
187 CLHEP::Hep3Vector horizontal_axis = y_axis.cross(beamCenterDiffAxis);
188 assert(horizontal_axis.mag() > CLHEP::Hep3Vector::getTolerance());
189 horizontal_axis = horizontal_axis / horizontal_axis.mag();
191 CLHEP::Hep3Vector vertical_axis = beamCenterDiffAxis.cross(horizontal_axis);
192 assert(vertical_axis.mag() > CLHEP::Hep3Vector::getTolerance());
193 vertical_axis = vertical_axis / vertical_axis.mag();
195 CLHEP::Hep3Vector vec_horizontal_collision_vertex_smear = horizontal_axis *
200 CLHEP::Hep3Vector vec_vertical_collision_vertex_smear = vertical_axis *
206 CLHEP::Hep3Vector vec_collision_vertex =
207 vec_horizontal_collision_vertex_smear +
208 vec_vertical_collision_vertex_smear +
209 vec_crossing_collision + vec_longitudinal_collision;
212 vec_collision_vertex.x(),
213 vec_collision_vertex.y(),
214 vec_collision_vertex.z(),
219 cout << __PRETTY_FUNCTION__
221 <<
"bunch_zs.first = " << bunch_zs.first <<
", "
222 <<
"bunch_zs.second = " << bunch_zs.second <<
", "
223 <<
"cos(theta/2) = " << beamCenterDiffAxis.dot(beamA_center) <<
", " << endl
225 <<
"beamCenterDiffAxis = " << beamCenterDiffAxis <<
", "
226 <<
"vec_crossing = " << vec_crossing <<
", "
227 <<
"horizontal_axis = " << horizontal_axis <<
", "
228 <<
"vertical_axis = " << vertical_axis <<
", " << endl
230 <<
"vec_longitudinal_collision = " << vec_longitudinal_collision <<
", "
231 <<
"vec_crossing_collision = " << vec_crossing_collision <<
", "
232 <<
"vec_vertical_collision_vertex_smear = " << vec_vertical_collision_vertex_smear <<
", "
233 <<
"vec_horizontal_collision_vertex_smear = " << vec_horizontal_collision_vertex_smear <<
", " << endl
234 <<
"vec_collision_vertex = " << vec_collision_vertex <<
", " << endl
236 <<
"ct_collision = " << ct_collision <<
", "
237 <<
"t_collision = " << t_collision <<
", "
246 const double &
theta = theta_phi.first;
247 const double &
phi = theta_phi.second;
249 return CLHEP::Hep3Vector(
250 sin(theta) * cos(phi),
251 sin(theta) * sin(phi),
276 cout <<
"PHHepMCGenHelper::HepMC2Lab_boost_rotation_translation - Fatal Error - the requested source subevent with embedding ID "
296 cout << __PRETTY_FUNCTION__ <<
": copied boost rotation shift of the collision" << endl;
304 pair<double, double> beam_bunch_zs;
315 beam_bunch_zs.first = beam_bunch_zs.second = init_vertex_longitudinal;
320 const static CLHEP::Hep3Vector z_axis(0, 0, 1);
327 cout << __PRETTY_FUNCTION__ <<
": " << endl;
328 cout <<
"beamA_center = " << beamA_center << endl;
329 cout <<
"beamB_center = " << beamB_center << endl;
332 assert(fabs(beamB_center.mag2() - 1) < CLHEP::Hep3Vector::getTolerance());
333 assert(fabs(beamB_center.mag2() - 1) < CLHEP::Hep3Vector::getTolerance());
335 if (beamA_center.dot(beamB_center) > -0.5)
337 cout <<
"PHHepMCGenHelper::HepMC2Lab_boost_rotation_translation - WARNING -"
338 <<
"Beam A and Beam B are not near back to back. "
339 <<
"Please double check beam direction setting at set_beam_direction_theta_phi()."
340 <<
"beamA_center = " << beamA_center <<
","
341 <<
"beamB_center = " << beamB_center <<
","
342 <<
" Current setting:";
348 auto smear_beam_divergence = [&,
this](
349 const CLHEP::Hep3Vector &beam_center,
350 const std::pair<double, double> &divergence_hv,
351 const std::pair<double, double> &beam_angular_z_coefficient_hv) {
352 const double &x_divergence = divergence_hv.first;
353 const double &y_divergence = divergence_hv.second;
356 static const CLHEP::Hep3Vector accelerator_plane(0, 1, 0);
358 CLHEP::Hep3Vector beam_direction(beam_center);
359 CLHEP::HepRotation x_smear_in_accelerator_plane(
362 beam_bunch_zs.first * beam_angular_z_coefficient_hv.first,
365 CLHEP::HepRotation y_smear_out_accelerator_plane(
366 accelerator_plane.cross(beam_center),
368 beam_bunch_zs.second * beam_angular_z_coefficient_hv.second,
372 return y_smear_out_accelerator_plane * x_smear_in_accelerator_plane * beam_center;
375 CLHEP::Hep3Vector beamA_vec = smear_beam_divergence(beamA_center,
378 CLHEP::Hep3Vector beamB_vec = smear_beam_divergence(beamB_center,
384 cout << __PRETTY_FUNCTION__ <<
": " << endl;
385 cout <<
"beamA_vec = " << beamA_vec << endl;
386 cout <<
"beamB_vec = " << beamB_vec << endl;
389 assert(fabs(beamA_vec.mag2() - 1) < CLHEP::Hep3Vector::getTolerance());
390 assert(fabs(beamB_vec.mag2() - 1) < CLHEP::Hep3Vector::getTolerance());
393 CLHEP::Hep3Vector boost_axis = beamA_vec + beamB_vec;
394 if (boost_axis.mag2() > CLHEP::Hep3Vector::getTolerance())
403 cout << __PRETTY_FUNCTION__ <<
": non-zero boost " << endl;
411 cout << __PRETTY_FUNCTION__ <<
": zero boost " << endl;
416 CLHEP::Hep3Vector beamDiffAxis = (beamA_vec - beamB_vec);
417 if (beamDiffAxis.mag2() < CLHEP::Hep3Vector::getTolerance())
419 cout <<
"PHHepMCGenHelper::HepMC2Lab_boost_rotation_translation - Fatal error -"
420 <<
"Beam A and Beam B are too close to each other in direction "
421 <<
"Please double check beam direction and divergence setting. "
422 <<
"beamA_vec = " << beamA_vec <<
","
423 <<
"beamB_vec = " << beamB_vec <<
","
424 <<
" Current setting:";
431 beamDiffAxis = beamDiffAxis / beamDiffAxis.mag();
432 double cos_rotation_angle_to_z = beamDiffAxis.dot(z_axis);
435 cout << __PRETTY_FUNCTION__ <<
": check rotation ";
436 cout <<
"cos_rotation_angle_to_z= " << cos_rotation_angle_to_z << endl;
439 if (1 - cos_rotation_angle_to_z < CLHEP::Hep3Vector::getTolerance())
447 cout << __PRETTY_FUNCTION__ <<
": no rotation " << endl;
450 else if (cos_rotation_angle_to_z + 1 < CLHEP::Hep3Vector::getTolerance())
457 cout << __PRETTY_FUNCTION__ <<
": reverse beam direction " << endl;
463 CLHEP::Hep3Vector rotation_axis = (beamA_vec - beamB_vec).
cross(z_axis);
464 const double rotation_angle_to_z = -acos(cos_rotation_angle_to_z);
471 cout << __PRETTY_FUNCTION__ <<
": has rotation " << endl;
477 cout << __PRETTY_FUNCTION__ <<
": final boost rotation shift of the collision" << endl;
486 cout << __PRETTY_FUNCTION__ <<
" Fatal Error: "
487 <<
"m_use_beam_bunch_sim = " <<
m_use_beam_bunch_sim <<
". Expect to simulate bunch interaction instead of applying vertex distributions"
502 cout << __PRETTY_FUNCTION__ <<
" Fatal Error: "
503 <<
"m_use_beam_bunch_sim = " <<
m_use_beam_bunch_sim <<
". Expect to simulate bunch interaction instead of applying vertex distributions"
519 cout << __PRETTY_FUNCTION__ <<
" Fatal Error: "
520 <<
"m_use_beam_bunch_sim = " <<
m_use_beam_bunch_sim <<
". Expect to simulate bunch interaction instead of applying vertex distributions"
536 cout << __PRETTY_FUNCTION__ <<
" Fatal Error: "
537 <<
"m_use_beam_bunch_sim = " << m_use_beam_bunch_sim <<
". Expect not to simulate bunch interaction but applying vertex distributions"
561 else if (dist ==
Gaus)
567 cout <<
"PHHepMCGenHelper::smear - FATAL Error - unknown vertex function " << dist << endl;
619 cout <<
"PHHepMCGenHelper::CopySettings - fatal error - invalid input class helper_dest which is nullptr!" << endl;
630 cout <<
"PHHepMCGenHelper::CopyHelperSettings - fatal error - invalid input class helper_src which is nullptr!" << endl;
637 static map<VTXFUNC, string> vtxfunc = {{VTXFUNC::Uniform,
"Uniform"}, {VTXFUNC::Gaus,
"Gaus"}};
645 cout <<
"Vertex distribution function x: " << vtxfunc[
_vertex_func_x]
668 cout <<
"Beam bunch A width = ["
670 cout <<
"Beam bunch B width = ["