generators/phhepmc/PHHepMCGenHelper.cc

0001 // $Id: $
0002
0003 /*!
0004  * \file PHHepMCGenHelper.cc
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include "PHHepMCGenHelper.h"
0012
0013 #include "PHHepMCGenEvent.h"
0014 #include "PHHepMCGenEventMap.h"
0015 #include "PHHepMCGenEventv1.h"
0016
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/PHIODataNode.h>    // for PHIODataNode
0021 #include <phool/PHNode.h>          // for PHNode
0022 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0023 #include <phool/PHObject.h>        // for PHObject
0024 #include <phool/PHRandomSeed.h>
0025 #include <phool/getClass.h>
0026 #include <phool/phool.h>  // for PHWHERE
0027
0028 #include <HepMC/SimpleVector.h>  // for FourVector
0029
0030 #include <CLHEP/Units/PhysicalConstants.h>
0031 #include <CLHEP/Units/SystemOfUnits.h>
0032 #include <CLHEP/Vector/Boost.h>
0033 #include <CLHEP/Vector/LorentzRotation.h>
0034 #include <CLHEP/Vector/LorentzVector.h>
0035 #include <CLHEP/Vector/Rotation.h>
0036
0037 #include <gsl/gsl_randist.h>
0038 #include <gsl/gsl_rng.h>
0039
0040 #include <cassert>
0041 #include <cmath>
0042 #include <cstdlib>  // for exit
0043 #include <iostream>
0044 #include <limits>
0045
0046 PHHepMCGenHelper::PHHepMCGenHelper()
0047   : RandomGenerator(gsl_rng_alloc(gsl_rng_mt19937))
0048 {
0049   unsigned int seed = PHRandomSeed();  // fixed seed is handled in this function
0050   gsl_rng_set(RandomGenerator, seed);
0051 }
0052
0053 PHHepMCGenHelper::~PHHepMCGenHelper()
0054 {
0055   gsl_rng_free(RandomGenerator);
0056 }
0057
0058 //! init interface nodes
0059 int PHHepMCGenHelper::create_node_tree(PHCompositeNode *topNode)
0060 {
0061   PHNodeIterator iter(topNode);
0062   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0063   if (!dstNode)
0064   {
0065     std::cout << PHWHERE << "DST Node missing doing nothing" << std::endl;
0066     return Fun4AllReturnCodes::ABORTRUN;
0067   }
0068
0069   _geneventmap = findNode::getClass<PHHepMCGenEventMap>(dstNode, "PHHepMCGenEventMap");
0070   if (!_geneventmap)
0071   {
0072     _geneventmap = new PHHepMCGenEventMap();
0073     PHIODataNode<PHObject> *newmapnode = new PHIODataNode<PHObject>(_geneventmap, "PHHepMCGenEventMap", "PHObject");
0074     dstNode->addNode(newmapnode);
0075   }
0076
0077   assert(_geneventmap);
0078
0079   return Fun4AllReturnCodes::EVENT_OK;
0080 }
0081
0082 //! choice of reference version of the PHHepMCGenEvent
0083 const PHHepMCGenEvent *PHHepMCGenHelper::get_PHHepMCGenEvent_template()
0084 {
0085   // choice of version of PHHepMCGenEvent
0086   const static PHHepMCGenEventv1 mc_evnet_template;
0087
0088   return &mc_evnet_template;
0089 }
0090
0091 //! send HepMC::GenEvent to DST tree. This function takes ownership of evt
0092 PHHepMCGenEvent *PHHepMCGenHelper::insert_event(HepMC::GenEvent *evt)
0093 {
0094   assert(evt);
0095   assert(_geneventmap);
0096
0097   PHHepMCGenEvent *genevent = _geneventmap->insert_event(_embedding_id, get_PHHepMCGenEvent_template());
0098
0099   genevent->addEvent(evt);
0100   HepMC2Lab_boost_rotation_translation(genevent);
0101
0102   return genevent;
0103 }
0104
0105 void PHHepMCGenHelper::move_vertex(PHHepMCGenEvent *genevent)
0106 {
0107   assert(genevent);
0108
0109   assert(_vertex_width_x >= 0);
0110   assert(_vertex_width_y >= 0);
0111   assert(_vertex_width_z >= 0);
0112   assert(_vertex_width_t >= 0);
0113
0114   assert(not _reuse_vertex);  // logic check
0115
0116   // not reusing vertex so smear with the vertex parameters
0117   genevent->moveVertex(
0118       (smear(_vertex_x, _vertex_width_x, _vertex_func_x)),
0119       (smear(_vertex_y, _vertex_width_y, _vertex_func_y)),
0120       (smear(_vertex_z, _vertex_width_z, _vertex_func_z)),
0121       (smear(_vertex_t, _vertex_width_t, _vertex_func_t)));
0122 }
0123
0124 //! use m_beam_bunch_width to calculate horizontal and vertical collision width
0125 //! \param[in] hv_index 0: horizontal. 1: vertical
0126 //! https://github.com/eic/documents/blob/d06b5597a0a89dcad215bab50fe3eefa17a097a5/reports/general/Note-Simulations-BeamEffects.pdf
0127 double PHHepMCGenHelper::get_collision_width(unsigned int hv_index)
0128 {
0129   assert((hv_index == 0) or (hv_index == 1));
0130
0131   const double widthA = m_beam_bunch_width.first[hv_index];
0132   const double widthB = m_beam_bunch_width.second[hv_index];
0133
0134   return widthA * widthB / sqrt((widthA * widthA) + (widthB * widthB));
0135 }
0136
0137 //! generate vertx with bunch interaction according to
0138 //! https://github.com/eic/documents/blob/d06b5597a0a89dcad215bab50fe3eefa17a097a5/reports/general/Note-Simulations-BeamEffects.pdf
0139 //! \return pair of bunch local z position for beam A and beam B
0140 std::pair<double, double> PHHepMCGenHelper::generate_vertx_with_bunch_interaction(PHHepMCGenEvent *genevent)
0141 {
0142   const std::pair<double, double> bunch_zs(
0143       smear(
0144           0,                            //  central vertical angle shift
0145           m_beam_bunch_width.first[2],  // vertical angle smear
0146           Gaus),
0147       smear(
0148           0,                             //  central vertical angle shift
0149           m_beam_bunch_width.second[2],  // vertical angle smear
0150           Gaus));
0151
0152   CLHEP::Hep3Vector beamA_center = pair2Hep3Vector(m_beam_direction_theta_phi.first);
0153   CLHEP::Hep3Vector beamB_center = pair2Hep3Vector(m_beam_direction_theta_phi.second);
0154   //  const static CLHEP::Hep3Vector z_axis(0, 0, 1);
0155   const static CLHEP::Hep3Vector y_axis(0, 1, 0);
0156
0157   // the final longitudinal vertex smear axis
0158   CLHEP::Hep3Vector beamCenterDiffAxis = (beamA_center - beamB_center);
0159   assert(beamCenterDiffAxis.mag() > CLHEP::Hep3Vector::getTolerance());
0160   beamCenterDiffAxis = beamCenterDiffAxis / beamCenterDiffAxis.mag();
0161
0162   CLHEP::Hep3Vector vec_crossing = beamA_center - 0.5 * (beamA_center - beamB_center);
0163
0164   CLHEP::Hep3Vector vec_longitudinal_collision = beamCenterDiffAxis * (bunch_zs.first + bunch_zs.second) / 2.;
0165   double ct_collision = 0.5 * (-bunch_zs.first + bunch_zs.second) / beamCenterDiffAxis.dot(beamA_center);
0166   double t_collision = ct_collision * CLHEP::cm / CLHEP::c_light / CLHEP::ns;
0167   CLHEP::Hep3Vector vec_crossing_collision = ct_collision * vec_crossing;  // shift of collision to crossing dierction
0168
0169   CLHEP::Hep3Vector horizontal_axis = y_axis.cross(beamCenterDiffAxis);
0170   assert(horizontal_axis.mag() > CLHEP::Hep3Vector::getTolerance());
0171   horizontal_axis = horizontal_axis / horizontal_axis.mag();
0172
0173   CLHEP::Hep3Vector vertical_axis = beamCenterDiffAxis.cross(horizontal_axis);
0174   assert(vertical_axis.mag() > CLHEP::Hep3Vector::getTolerance());
0175   vertical_axis = vertical_axis / vertical_axis.mag();
0176
0177   CLHEP::Hep3Vector vec_horizontal_collision_vertex_smear = horizontal_axis *
0178                                                             smear(
0179                                                                 0,
0180                                                                 get_collision_width(0),
0181                                                                 Gaus);
0182   CLHEP::Hep3Vector vec_vertical_collision_vertex_smear = vertical_axis *
0183                                                           smear(
0184                                                               0,
0185                                                               get_collision_width(1),
0186                                                               Gaus);
0187
0188   CLHEP::Hep3Vector vec_collision_vertex =
0189       vec_horizontal_collision_vertex_smear +
0190       vec_vertical_collision_vertex_smear +  //
0191       vec_crossing_collision + vec_longitudinal_collision;
0192
0193   genevent->set_collision_vertex(HepMC::FourVector(
0194       vec_collision_vertex.x(),
0195       vec_collision_vertex.y(),
0196       vec_collision_vertex.z(),
0197       t_collision));
0198
0199   if (m_verbosity)
0200   {
0201     std::cout << __PRETTY_FUNCTION__
0202               << ":"
0203               << "bunch_zs.first  = " << bunch_zs.first << ", "
0204               << "bunch_zs.second = " << bunch_zs.second << ", "
0205               << "cos(theta/2) = " << beamCenterDiffAxis.dot(beamA_center) << ", " << std::endl
0206
0207               << "beamCenterDiffAxis = " << beamCenterDiffAxis << ", "
0208               << "vec_crossing = " << vec_crossing << ", "
0209               << "horizontal_axis = " << horizontal_axis << ", "
0210               << "vertical_axis = " << vertical_axis << ", " << std::endl
0211
0212               << "vec_longitudinal_collision = " << vec_longitudinal_collision << ", "
0213               << "vec_crossing_collision = " << vec_crossing_collision << ", "
0214               << "vec_vertical_collision_vertex_smear = " << vec_vertical_collision_vertex_smear << ", "
0215               << "vec_horizontal_collision_vertex_smear = " << vec_horizontal_collision_vertex_smear << ", " << std::endl
0216               << "vec_collision_vertex = " << vec_collision_vertex << ", " << std::endl
0217
0218               << "ct_collision = " << ct_collision << ", "
0219               << "t_collision = " << t_collision << ", "
0220               << std::endl;
0221   }
0222
0223   return bunch_zs;
0224 }
0225
0226 CLHEP::Hep3Vector PHHepMCGenHelper::pair2Hep3Vector(const std::pair<double, double> &theta_phi)
0227 {
0228   const double &theta = theta_phi.first;
0229   const double &phi = theta_phi.second;
0230
0231   return CLHEP::Hep3Vector(
0232       sin(theta) * cos(phi),
0233       sin(theta) * sin(phi),
0234       cos(theta));
0235 }
0236
0237 //! move vertex in translation,boost,rotation according to vertex settings
0238 void PHHepMCGenHelper::HepMC2Lab_boost_rotation_translation(PHHepMCGenEvent *genevent)
0239 {
0240   if (m_verbosity)
0241   {
0242     Print();
0243   }
0244
0245   assert(genevent);
0246
0247   if (_reuse_vertex)
0248   {
0249     // just copy over the vertex boost_rotation_translation
0250
0251     assert(_geneventmap);
0252
0253     PHHepMCGenEvent *vtx_evt =
0254         _geneventmap->get(_reuse_vertex_embedding_id);
0255
0256     if (!vtx_evt)
0257     {
0258       std::cout << "PHHepMCGenHelper::HepMC2Lab_boost_rotation_translation - Fatal Error - the requested source subevent with embedding ID "
0259                 << _reuse_vertex_embedding_id << " does not exist. Current HepMCEventMap:";
0260       _geneventmap->identify();
0261       exit(1);
0262     }
0263
0264     // copy boost_rotation_translation
0265
0266     genevent->moveVertex(
0267         vtx_evt->get_collision_vertex().x(),
0268         vtx_evt->get_collision_vertex().y(),
0269         vtx_evt->get_collision_vertex().z(),
0270         vtx_evt->get_collision_vertex().t());
0271
0272     genevent->set_boost_beta_vector(vtx_evt->get_boost_beta_vector());
0273     genevent->set_rotation_vector(vtx_evt->get_rotation_vector());
0274     genevent->set_rotation_angle(vtx_evt->get_rotation_angle());
0275
0276     if (m_verbosity)
0277     {
0278       std::cout << __PRETTY_FUNCTION__ << ": copied boost rotation shift of the collision" << std::endl;
0279       genevent->identify();
0280     }
0281     return;
0282   }  //!   if (_reuse_vertex)
0283
0284   // now handle the collision vertex first, in the head-on collision frame
0285   // this is used as input to the Crab angle correction
0286   std::pair<double, double> beam_bunch_zs;
0287   if (m_use_beam_bunch_sim)
0288   {
0289     // bunch interaction simulation
0290     beam_bunch_zs = generate_vertx_with_bunch_interaction(genevent);
0291   }
0292   else
0293   {
0294     // vertex distribution simulation
0295     move_vertex(genevent);
0296     const double init_vertex_longitudinal = genevent->get_collision_vertex().z();
0297     beam_bunch_zs.first = beam_bunch_zs.second = init_vertex_longitudinal;
0298   }
0299
0300   // boost-rotation from beam angles
0301
0302   const static CLHEP::Hep3Vector z_axis(0, 0, 1);
0303
0304   CLHEP::Hep3Vector beamA_center = pair2Hep3Vector(m_beam_direction_theta_phi.first);
0305   CLHEP::Hep3Vector beamB_center = pair2Hep3Vector(m_beam_direction_theta_phi.second);
0306
0307   if (m_verbosity)
0308   {
0309     std::cout << __PRETTY_FUNCTION__ << ": " << std::endl;
0310     std::cout << "beamA_center = " << beamA_center << std::endl;
0311     std::cout << "beamB_center = " << beamB_center << std::endl;
0312   }
0313
0314   assert(fabs(beamB_center.mag2() - 1) < CLHEP::Hep3Vector::getTolerance());
0315   assert(fabs(beamB_center.mag2() - 1) < CLHEP::Hep3Vector::getTolerance());
0316
0317   if (beamA_center.dot(beamB_center) > -0.5)
0318   {
0319     std::cout << "PHHepMCGenHelper::HepMC2Lab_boost_rotation_translation - WARNING -"
0320               << "Beam A and Beam B are not near back to back. "
0321               << "Please double check beam direction setting at set_beam_direction_theta_phi()."
0322               << "beamA_center = " << beamA_center << ","
0323               << "beamB_center = " << beamB_center << ","
0324               << " Current setting:";
0325
0326     Print();
0327   }
0328
0329   // function to do angular shifts relative to central beam angle
0330   auto smear_beam_divergence = [&, this](
0331                                    const CLHEP::Hep3Vector &beam_center,
0332                                    const std::pair<double, double> &divergence_hv,
0333                                    const std::pair<double, double> &beam_angular_z_coefficient_hv)
0334   {
0335     const double &x_divergence = divergence_hv.first;
0336     const double &y_divergence = divergence_hv.second;
0337
0338     // y direction in accelerator
0339     static const CLHEP::Hep3Vector accelerator_plane(0, 1, 0);
0340
0341     //    CLHEP::Hep3Vector beam_direction(beam_center);
0342     CLHEP::HepRotation x_smear_in_accelerator_plane(
0343         accelerator_plane,
0344         smear(
0345             beam_bunch_zs.first * beam_angular_z_coefficient_hv.first,  //  central horizontal angle shift
0346             x_divergence,                                               // horizontal angle smear
0347             Gaus));
0348     CLHEP::HepRotation y_smear_out_accelerator_plane(
0349         accelerator_plane.cross(beam_center),
0350         smear(
0351             beam_bunch_zs.second * beam_angular_z_coefficient_hv.second,  //  central vertical angle shift
0352             y_divergence,                                                 // vertical angle smear
0353             Gaus));
0354
0355     return y_smear_out_accelerator_plane * x_smear_in_accelerator_plane * beam_center;
0356   };
0357
0358   CLHEP::Hep3Vector beamA_vec = smear_beam_divergence(beamA_center,
0359                                                       m_beam_angular_divergence_hv.first,
0360                                                       m_beam_angular_z_coefficient_hv.first);
0361   CLHEP::Hep3Vector beamB_vec = smear_beam_divergence(beamB_center,
0362                                                       m_beam_angular_divergence_hv.second,
0363                                                       m_beam_angular_z_coefficient_hv.second);
0364
0365   if (m_verbosity)
0366   {
0367     std::cout << __PRETTY_FUNCTION__ << ": " << std::endl;
0368     std::cout << "beamA_vec = " << beamA_vec << std::endl;
0369     std::cout << "beamB_vec = " << beamB_vec << std::endl;
0370   }
0371
0372   assert(fabs(beamA_vec.mag2() - 1) < CLHEP::Hep3Vector::getTolerance());
0373   assert(fabs(beamB_vec.mag2() - 1) < CLHEP::Hep3Vector::getTolerance());
0374
0375   // apply minimal beam energy shift rotation and boost
0376   CLHEP::Hep3Vector boost_axis = beamA_vec + beamB_vec;
0377   if (boost_axis.mag2() > CLHEP::Hep3Vector::getTolerance())
0378   {
0379     // non-zero boost
0380
0381     // split the boost to half for each beam for minimal beam  energy shift
0382     genevent->set_boost_beta_vector(0.5 * boost_axis);
0383
0384     if (m_verbosity)
0385     {
0386       std::cout << __PRETTY_FUNCTION__ << ": non-zero boost " << std::endl;
0387     }
0388   }  //    if (cos_rotation_angle> CLHEP::Hep3Vector::getTolerance())
0389   else
0390   {
0391     genevent->set_boost_beta_vector(CLHEP::Hep3Vector(0, 0, 0));
0392     if (m_verbosity)
0393     {
0394       std::cout << __PRETTY_FUNCTION__ << ": zero boost " << std::endl;
0395     }
0396   }
0397
0398   // rotation to collision to along z-axis with beamA pointing to +z
0399   CLHEP::Hep3Vector beamDiffAxis = (beamA_vec - beamB_vec);
0400   if (beamDiffAxis.mag2() < CLHEP::Hep3Vector::getTolerance())
0401   {
0402     std::cout << "PHHepMCGenHelper::HepMC2Lab_boost_rotation_translation - Fatal error -"
0403               << "Beam A and Beam B are too close to each other in direction "
0404               << "Please double check beam direction and divergence setting. "
0405               << "beamA_vec = " << beamA_vec << ","
0406               << "beamB_vec = " << beamB_vec << ","
0407               << " Current setting:";
0408
0409     Print();
0410
0411     exit(1);
0412   }
0413
0414   beamDiffAxis = beamDiffAxis / beamDiffAxis.mag();
0415   double cos_rotation_angle_to_z = beamDiffAxis.dot(z_axis);
0416   if (m_verbosity)
0417   {
0418     std::cout << __PRETTY_FUNCTION__ << ": check rotation ";
0419     std::cout << "cos_rotation_angle_to_z= " << cos_rotation_angle_to_z << std::endl;
0420   }
0421
0422   if (1 - cos_rotation_angle_to_z < CLHEP::Hep3Vector::getTolerance())
0423   {
0424     // no rotation
0425     genevent->set_rotation_vector(z_axis);
0426     genevent->set_rotation_angle(0);
0427
0428     if (m_verbosity)
0429     {
0430       std::cout << __PRETTY_FUNCTION__ << ": no rotation " << std::endl;
0431     }
0432   }
0433   else if (cos_rotation_angle_to_z + 1 < CLHEP::Hep3Vector::getTolerance())
0434   {
0435     // you got beam flipped
0436     genevent->set_rotation_vector(CLHEP::Hep3Vector(0, 1, 0));
0437     genevent->set_rotation_angle(M_PI);
0438     if (m_verbosity)
0439     {
0440       std::cout << __PRETTY_FUNCTION__ << ": reverse beam direction " << std::endl;
0441     }
0442   }
0443   else
0444   {
0445     // need a rotation
0446     CLHEP::Hep3Vector rotation_axis = (beamA_vec - beamB_vec).cross(z_axis);
0447     const double rotation_angle_to_z = -acos(cos_rotation_angle_to_z);
0448
0449     genevent->set_rotation_vector(rotation_axis);
0450     genevent->set_rotation_angle(rotation_angle_to_z);
0451
0452     if (m_verbosity)
0453     {
0454       std::cout << __PRETTY_FUNCTION__ << ": has rotation " << std::endl;
0455     }
0456   }  //  if (boost_axis.mag2() > CLHEP::Hep3Vector::getTolerance())
0457
0458   if (m_verbosity)
0459   {
0460     std::cout << __PRETTY_FUNCTION__ << ": final boost rotation shift of the collision" << std::endl;
0461     genevent->identify();
0462   }
0463 }
0464
0465 void PHHepMCGenHelper::set_vertex_distribution_function(VTXFUNC x, VTXFUNC y, VTXFUNC z, VTXFUNC t)
0466 {
0467   if (m_use_beam_bunch_sim)
0468   {
0469     std::cout << __PRETTY_FUNCTION__ << " Fatal Error: "
0470               << "m_use_beam_bunch_sim = " << m_use_beam_bunch_sim << ". Expect to simulate bunch interaction instead of applying vertex distributions"
0471               << std::endl;
0472     exit(1);
0473   }
0474   _vertex_func_x = x;
0475   _vertex_func_y = y;
0476   _vertex_func_z = z;
0477   _vertex_func_t = t;
0478   return;
0479 }
0480
0481 void PHHepMCGenHelper::set_vertex_distribution_mean(const double x, const double y, const double z, const double t)
0482 {
0483   if (m_use_beam_bunch_sim)
0484   {
0485     std::cout << __PRETTY_FUNCTION__ << " Fatal Error: "
0486               << "m_use_beam_bunch_sim = " << m_use_beam_bunch_sim << ". Expect to simulate bunch interaction instead of applying vertex distributions"
0487               << std::endl;
0488     exit(1);
0489   }
0490
0491   _vertex_x = x;
0492   _vertex_y = y;
0493   _vertex_z = z;
0494   _vertex_t = t;
0495   return;
0496 }
0497
0498 void PHHepMCGenHelper::set_vertex_distribution_width(const double x, const double y, const double z, const double t)
0499 {
0500   if (m_use_beam_bunch_sim)
0501   {
0502     std::cout << __PRETTY_FUNCTION__ << " Fatal Error: "
0503               << "m_use_beam_bunch_sim = " << m_use_beam_bunch_sim << ". Expect to simulate bunch interaction instead of applying vertex distributions"
0504               << std::endl;
0505     exit(1);
0506   }
0507
0508   _vertex_width_x = x;
0509   _vertex_width_y = y;
0510   _vertex_width_z = z;
0511   _vertex_width_t = t;
0512   return;
0513 }
0514
0515 void PHHepMCGenHelper::set_beam_bunch_width(const std::vector<double> &beamA, const std::vector<double> &beamB)
0516 {
0517   if (not m_use_beam_bunch_sim)
0518   {
0519     std::cout << __PRETTY_FUNCTION__ << " Fatal Error: "
0520               << "m_use_beam_bunch_sim = " << m_use_beam_bunch_sim << ". Expect not to simulate bunch interaction but applying vertex distributions"
0521               << std::endl;
0522     exit(1);
0523   }
0524
0525   m_beam_bunch_width.first = beamA;
0526   m_beam_bunch_width.second = beamB;
0527 }
0528
0529 double PHHepMCGenHelper::smear(const double position,
0530                                const double width,
0531                                VTXFUNC dist) const
0532 {
0533   assert(width >= 0);
0534
0535   double res = position;
0536
0537   if (width == 0)
0538   {
0539     return res;
0540   }
0541
0542   if (dist == Uniform)
0543   {
0544     res = (position - width) + 2 * gsl_rng_uniform_pos(RandomGenerator) * width;
0545   }
0546   else if (dist == Gaus)
0547   {
0548     res = position + gsl_ran_gaussian(RandomGenerator, width);
0549   }
0550   else
0551   {
0552     std::cout << "PHHepMCGenHelper::smear - FATAL Error - unknown vertex function " << dist << std::endl;
0553     exit(10);
0554   }
0555   return res;
0556 }
0557
0558 void PHHepMCGenHelper::CopySettings(PHHepMCGenHelper &helper_dest)
0559 {
0560   // allow copy of vertex distributions
0561   helper_dest.use_beam_bunch_sim(false);
0562   helper_dest.set_vertex_distribution_width(_vertex_width_x, _vertex_width_y, _vertex_width_z, _vertex_width_t);
0563   helper_dest.set_vertex_distribution_function(_vertex_func_x, _vertex_func_y, _vertex_func_z, _vertex_func_t);
0564   helper_dest.set_vertex_distribution_mean(_vertex_x, _vertex_y, _vertex_z, _vertex_t);
0565
0566   // allow copy of bunch distributions
0567   helper_dest.use_beam_bunch_sim(true);
0568   helper_dest.set_beam_bunch_width(m_beam_bunch_width.first, m_beam_bunch_width.second);
0569
0570   // final bunch settings
0571   helper_dest.use_beam_bunch_sim(m_use_beam_bunch_sim);
0572
0573   helper_dest.set_beam_direction_theta_phi(
0574       m_beam_direction_theta_phi.first.first,
0575       m_beam_direction_theta_phi.first.second,
0576       m_beam_direction_theta_phi.second.first,
0577       m_beam_direction_theta_phi.second.second);
0578   helper_dest.set_beam_angular_divergence_hv(
0579       m_beam_angular_divergence_hv.first.first,
0580       m_beam_angular_divergence_hv.first.second,
0581       m_beam_angular_divergence_hv.second.first,
0582       m_beam_angular_divergence_hv.second.second);
0583   helper_dest.set_beam_angular_z_coefficient_hv(
0584       m_beam_angular_z_coefficient_hv.first.first,
0585       m_beam_angular_z_coefficient_hv.first.second,
0586       m_beam_angular_z_coefficient_hv.second.first,
0587       m_beam_angular_z_coefficient_hv.second.second);
0588
0589   helper_dest.set_beam_angular_z_coefficient_hv(
0590       m_beam_angular_z_coefficient_hv.first.first,
0591       m_beam_angular_z_coefficient_hv.first.second,
0592       m_beam_angular_z_coefficient_hv.second.first,
0593       m_beam_angular_z_coefficient_hv.second.second);
0594
0595   return;
0596 }
0597
0598 void PHHepMCGenHelper::CopySettings(PHHepMCGenHelper *helper_dest)
0599 {
0600   if (helper_dest)
0601   {
0602     CopySettings(*helper_dest);
0603   }
0604   else
0605   {
0606     std::cout << "PHHepMCGenHelper::CopySettings - fatal error - invalid input class helper_dest which is nullptr!" << std::endl;
0607     exit(1);
0608   }
0609 }
0610
0611 void PHHepMCGenHelper::CopyHelperSettings(PHHepMCGenHelper *helper_src)
0612 {
0613   if (helper_src)
0614   {
0615     helper_src->CopySettings(this);
0616   }
0617   else
0618   {
0619     std::cout << "PHHepMCGenHelper::CopyHelperSettings - fatal error - invalid input class helper_src which is nullptr!" << std::endl;
0620     exit(1);
0621   }
0622 }
0623
0624 void PHHepMCGenHelper::Print(const std::string & /*what*/) const
0625 {
0626   static std::map<VTXFUNC, std::string> vtxfunc = {{VTXFUNC::Uniform, "Uniform"}, {VTXFUNC::Gaus, "Gaus"}};
0627
0628   std::cout << "Vertex distribution width x: " << _vertex_width_x
0629             << ", y: " << _vertex_width_y
0630             << ", z: " << _vertex_width_z
0631             << ", t: " << _vertex_width_t
0632             << std::endl;
0633
0634   std::cout << "Vertex distribution function x: " << vtxfunc[_vertex_func_x]
0635             << ", y: " << vtxfunc[_vertex_func_y]
0636             << ", z: " << vtxfunc[_vertex_func_z]
0637             << ", t: " << vtxfunc[_vertex_func_t]
0638             << std::endl;
0639
0640   std::cout << "Beam direction: A  theta-phi = " << m_beam_direction_theta_phi.first.first
0641             << ", " << m_beam_direction_theta_phi.first.second << std::endl;
0642   std::cout << "Beam direction: B  theta-phi = " << m_beam_direction_theta_phi.second.first
0643             << ", " << m_beam_direction_theta_phi.second.second << std::endl;
0644
0645   std::cout << "Beam divergence: A X-Y = " << m_beam_angular_divergence_hv.first.first
0646             << ", " << m_beam_angular_divergence_hv.first.second << std::endl;
0647   std::cout << "Beam divergence: B X-Y = " << m_beam_angular_divergence_hv.second.first
0648             << ", " << m_beam_angular_divergence_hv.second.second << std::endl;
0649
0650   std::cout << "Beam angle shift as linear function of longitudinal vertex position : A X-Y = " << m_beam_angular_z_coefficient_hv.first.first
0651             << ", " << m_beam_angular_z_coefficient_hv.first.second << std::endl;
0652   std::cout << "Beam angle shift as linear function of longitudinal vertex position: B X-Y = " << m_beam_angular_z_coefficient_hv.second.first
0653             << ", " << m_beam_angular_z_coefficient_hv.second.second << std::endl;
0654
0655   std::cout << "m_use_beam_bunch_sim = " << m_use_beam_bunch_sim << std::endl;
0656
0657   std::cout << "Beam bunch A width = ["
0658             << m_beam_bunch_width.first[0] << ", " << m_beam_bunch_width.first[1] << ", " << m_beam_bunch_width.first[2] << "] cm" << std::endl;
0659   std::cout << "Beam bunch B width = ["
0660             << m_beam_bunch_width.second[0] << ", " << m_beam_bunch_width.second[1] << ", " << m_beam_bunch_width.second[2] << "] cm" << std::endl;
0661
0662   return;
0663 }