generators/ReactionPlaneAfterburner/ReactionPlaneAfterburner.cc

0001
0002 #include "ReactionPlaneAfterburner.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/PHRandomSeed.h>
0008 #include <phool/getClass.h>
0009 // phhepmc
0010 #include <phhepmc/PHHepMCGenEvent.h>
0011 #include <phhepmc/PHHepMCGenEventMap.h>
0012
0013 // hepmc
0014 #pragma GCC diagnostic push
0015 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0016 #include <HepMC/GenEvent.h>
0017 #include <HepMC/GenParticle.h>  // for GenParticle
0018 #include <HepMC/GenVertex.h>    // for GenVertex, GenVertex::part...
0019 #include <HepMC/HeavyIon.h>     // for HeavyIon
0020 #include <HepMC/SimpleVector.h>
0021 #pragma GCC diagnostic pop
0022
0023 #include <cassert>
0024
0025 //____________________________________________________________________________..
0026 ReactionPlaneAfterburner::ReactionPlaneAfterburner(const std::string &name)
0027   : SubsysReco(name)
0028   , RandomGenerator(gsl_rng_alloc(gsl_rng_mt19937))
0029 {
0030 }
0031
0032 //____________________________________________________________________________..
0033 ReactionPlaneAfterburner::~ReactionPlaneAfterburner()
0034 {
0035   gsl_rng_free(RandomGenerator);
0036 }
0037
0038 //____________________________________________________________________________..
0039 int ReactionPlaneAfterburner::Init(PHCompositeNode * /*topNode*/)
0040 {
0041   unsigned int seed = PHRandomSeed();
0042   gsl_rng_set(RandomGenerator, seed);
0043   return Fun4AllReturnCodes::EVENT_OK;
0044 }
0045
0046 //____________________________________________________________________________..
0047 int ReactionPlaneAfterburner::process_event(PHCompositeNode *topNode)
0048 {
0049   // get event
0050   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0051   for (auto &iter : *genevtmap)
0052   {
0053     PHHepMCGenEvent *genevt = iter.second;
0054     HepMC::GenEvent *evt = genevt->getEvent();
0055     assert(evt);
0056     HepMC::HeavyIon *hi = evt->heavy_ion();
0057     if (hi)
0058     {
0059       double psi = hi->event_plane_angle();
0060       if (psi != 0)
0061       {
0062         std::cout << "ReactionPlaneAfterburner::process_event(PHCompositeNode *topNode) psi = " << psi << std::endl;
0063         std::cout << "non-zero psi found, skipping" << std::endl;
0064         return Fun4AllReturnCodes::EVENT_OK;
0065       }
0066       psi = gsl_rng_uniform_pos(RandomGenerator) * 2 * M_PI;
0067       hi->set_event_plane_angle(psi);
0068       // rotate all particles and vertices
0069       for (HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p)
0070       {
0071         HepMC::FourVector v = (*p)->momentum();
0072         double px = v.px();
0073         double py = v.py();
0074         v.setPx((px * cos(psi)) - (py * sin(psi)));
0075         v.setPy((px * sin(psi)) + (py * cos(psi)));
0076         (*p)->set_momentum(v);
0077       }
0078       for (HepMC::GenEvent::vertex_iterator v = evt->vertices_begin(); v != evt->vertices_end(); ++v)
0079       {
0080         HepMC::FourVector pos = (*v)->position();
0081         double x = pos.x();
0082         double y = pos.y();
0083         pos.setX((x * cos(psi)) - (y * sin(psi)));
0084         pos.setY((x * sin(psi)) + (y * cos(psi)));
0085         (*v)->set_position(pos);
0086       }
0087     }
0088     else
0089     {
0090       std::cout << "ReactionPlaneAfterburner::process_event: no heavy ion info found, exiting" << std::endl;
0091       return Fun4AllReturnCodes::EVENT_OK;
0092     }
0093   }
0094   return Fun4AllReturnCodes::EVENT_OK;
0095 }