g4simulation/g4main/PHG4ParticleGeneratorBase.cc

0001 #include "PHG4ParticleGeneratorBase.h"
0002
0003 #include "PHG4Particle.h"  // for PHG4Particle
0004 #include "PHG4Particlev1.h"
0005
0006 #include "PHG4InEvent.h"
0007 #include "PHG4TruthInfoContainer.h"
0008 #include "PHG4VtxPoint.h"
0009
0010 #include <phhepmc/PHHepMCGenEvent.h>
0011 #include <phhepmc/PHHepMCGenEventMap.h>
0012
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/PHDataNode.h>      // for PHDataNode
0015 #include <phool/PHNode.h>          // for PHNode
0016 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0017 #include <phool/PHObject.h>        // for PHObject
0018 #include <phool/PHRandomSeed.h>
0019 #include <phool/getClass.h>
0020 #include <phool/phool.h>  // for PHWHERE
0021
0022 #include <TSystem.h>
0023
0024 #include <Geant4/G4ParticleDefinition.hh>
0025 #include <Geant4/G4ParticleTable.hh>
0026 #include <Geant4/G4SystemOfUnits.hh>
0027
0028 #include <HepMC/SimpleVector.h>  // for FourVector
0029
0030 #include <gsl/gsl_rng.h>
0031
0032 #include <cassert>
0033 #include <cstdlib>   // for exit
0034 #include <iostream>  // for operator<<, basic_ostream
0035 #include <iterator>  // for operator!=, reverse_iterator
0036 #include <map>       // for map<>::const_iterator, map
0037 #include <utility>   // for pair
0038
0039 PHG4ParticleGeneratorBase::PHG4ParticleGeneratorBase(const std::string &name)
0040   : SubsysReco(name)
0041   , m_RandomGenerator(gsl_rng_alloc(gsl_rng_mt19937))
0042 {
0043   m_Seed = PHRandomSeed();  // fixed seed is handled in this funtcion
0044   gsl_rng_set(m_RandomGenerator, m_Seed);
0045   return;
0046 }
0047
0048 PHG4ParticleGeneratorBase::~PHG4ParticleGeneratorBase()
0049 {
0050   while (particlelist.begin() != particlelist.end())
0051   {
0052     delete particlelist.back();
0053     particlelist.pop_back();
0054   }
0055   gsl_rng_free(m_RandomGenerator);
0056   return;
0057 }
0058
0059 int PHG4ParticleGeneratorBase::get_pdgcode(const std::string &name)
0060 {
0061   G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
0062   G4ParticleDefinition *particledef = particleTable->FindParticle(name);
0063   if (particledef)
0064   {
0065     return particledef->GetPDGEncoding();
0066   }
0067   return 0;
0068 }
0069
0070 std::string
0071 PHG4ParticleGeneratorBase::get_pdgname(const int pdgcode)
0072 {
0073   G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
0074   G4ParticleDefinition *particledef = particleTable->FindParticle(pdgcode);
0075   if (particledef)
0076   {
0077     return particledef->GetParticleName();
0078   }
0079   // if we cannot find the particle definition we'll make it ia geantino
0080   return "geantino";
0081 }
0082
0083 double
0084 PHG4ParticleGeneratorBase::get_mass(const int pdgcode)
0085 {
0086   G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
0087   G4ParticleDefinition *particledef = particleTable->FindParticle(get_pdgname(pdgcode));
0088   if (particledef)
0089   {
0090     return particledef->GetPDGMass() / GeV;
0091   }
0092   return 0;
0093 }
0094
0095 void PHG4ParticleGeneratorBase::set_name(const std::string &particle)
0096 {
0097   CheckAndCreateParticleVector();
0098   particlelist[0]->set_name(particle);
0099   particlelist[0]->set_pid(get_pdgcode(particle));
0100   return;
0101 }
0102
0103 void PHG4ParticleGeneratorBase::set_pid(const int pid)
0104 {
0105   CheckAndCreateParticleVector();
0106   particlelist[0]->set_pid(pid);
0107 }
0108
0109 void PHG4ParticleGeneratorBase::set_mom(const double x, const double y, const double z)
0110 {
0111   CheckAndCreateParticleVector();
0112   particlelist[0]->set_px(x);
0113   particlelist[0]->set_py(y);
0114   particlelist[0]->set_pz(z);
0115   return;
0116 }
0117
0118 void PHG4ParticleGeneratorBase::set_vtx(const double x, const double y, const double z)
0119 {
0120   m_Vtx_x = x;
0121   m_Vtx_y = y;
0122   m_Vtx_z = z;
0123   return;
0124 }
0125
0126 int PHG4ParticleGeneratorBase::InitRun(PHCompositeNode *topNode)
0127 {
0128   PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0129   if (!ineve)
0130   {
0131     PHNodeIterator iter(topNode);
0132     PHCompositeNode *dstNode;
0133     dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0134
0135     ineve = new PHG4InEvent();
0136     PHDataNode<PHObject> *newNode = new PHDataNode<PHObject>(ineve, "PHG4INEVENT", "PHObject");
0137     dstNode->addNode(newNode);
0138   }
0139   return 0;
0140 }
0141
0142 int PHG4ParticleGeneratorBase::process_event(PHCompositeNode * /*topNode*/)
0143 {
0144   std::cout << PHWHERE << " " << Name() << " using empty process_event" << std::endl;
0145   return 0;
0146 }
0147
0148 void PHG4ParticleGeneratorBase::PrintParticles(const std::string & /*what*/) const
0149 {
0150   std::vector<PHG4Particle *>::const_iterator iter;
0151   int i = 0;
0152   for (iter = particlelist.begin(); iter != particlelist.end(); ++iter)
0153   {
0154     std::cout << "particle " << i << std::endl;
0155     (*iter)->identify();
0156     i++;
0157   }
0158 }
0159
0160 void PHG4ParticleGeneratorBase::AddParticle(const std::string &particle, const double x, const double y, const double z)
0161 {
0162   PHG4Particle *part = new PHG4Particlev1(particle, get_pdgcode(particle), x, y, z);
0163   particlelist.push_back(part);
0164 }
0165
0166 void PHG4ParticleGeneratorBase::AddParticle(const int pid, const double x, const double y, const double z)
0167 {
0168   PHG4Particle *particle = new PHG4Particlev1();
0169   particle->set_pid(pid);
0170   particle->set_px(x);
0171   particle->set_py(y);
0172   particle->set_pz(z);
0173   particlelist.push_back(particle);
0174 }
0175
0176 void PHG4ParticleGeneratorBase::CheckAndCreateParticleVector()
0177 {
0178   if (particlelist.empty())
0179   {
0180     PHG4Particle *part = new PHG4Particlev1();
0181     particlelist.push_back(part);
0182   }
0183   return;
0184 }
0185
0186 void PHG4ParticleGeneratorBase::SetParticleId(PHG4Particle *particle, PHG4InEvent *ineve) const
0187 {
0188   if (particle->get_name().empty())  // no size -> empty name string
0189   {
0190     particle->set_name(get_pdgname(particle->get_pid()));
0191   }
0192   if (particle->get_pid() == 0)
0193   {
0194     particle->set_pid(get_pdgcode(particle->get_name()));
0195   }
0196   if (m_EmbedFlag)
0197   {
0198     ineve->AddEmbeddedParticle(particle, m_EmbedFlag);
0199   }
0200   return;
0201 }
0202
0203 void PHG4ParticleGeneratorBase::set_seed(const unsigned int iseed)
0204 {
0205   m_Seed = iseed;
0206   std::cout << Name() << " random seed: " << m_Seed << std::endl;
0207   gsl_rng_set(m_RandomGenerator, m_Seed);
0208 }
0209
0210 int PHG4ParticleGeneratorBase::ReuseExistingVertex(PHCompositeNode *topNode)
0211 {
0212   if (!m_ReUseExistingVertexFlag)
0213   {
0214     return 0;
0215   }
0216
0217   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0218
0219   if (genevtmap)
0220   {
0221     // use the highest priority HepMC subevent's vertex, ie that with highest embedding ID
0222
0223     PHHepMCGenEventMap::ConstReverseIter iter =
0224         genevtmap->rbegin();
0225
0226     if (iter != genevtmap->rend())
0227     {
0228       const PHHepMCGenEvent *hepmc_evt = iter->second;
0229
0230       assert(hepmc_evt);
0231
0232       const HepMC::FourVector &vtx = hepmc_evt->get_collision_vertex();
0233
0234       set_vtx(vtx.x(), vtx.y(), vtx.z());
0235
0236       if (Verbosity() > 0)
0237       {
0238         std::cout << "PHG4ParticleGeneratorBase::ReuseExistingVertex - reuse PHHepMCGenEventMap vertex "
0239                   << vtx.x() << ", " << vtx.y() << ", " << vtx.z() << " cm. Source event:"
0240                   << std::endl;
0241         hepmc_evt->identify();
0242       }
0243
0244       return 1;
0245     }
0246   }
0247
0248   // try PHG4INEVENT
0249   PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0250
0251   if (ineve->GetNVtx() > 0)
0252   {
0253     std::pair<std::map<int, PHG4VtxPoint *>::const_iterator,
0254               std::map<int, PHG4VtxPoint *>::const_iterator>
0255         range = ineve->GetVertices();
0256     std::map<int, PHG4VtxPoint *>::const_iterator iter = range.first;
0257     PHG4VtxPoint *vtx = iter->second;
0258
0259     if (!vtx)
0260     {
0261       std::cout << PHWHERE << "::Error - PHG4SimpleEventGenerator expects an existing vertex in PHG4InEvent, but none exists" << std::endl;
0262       exit(1);
0263     }
0264     if (Verbosity() > 0)
0265     {
0266       std::cout << PHWHERE << "::Info - use this primary vertex from PHG4InEvent:" << std::endl;
0267       vtx->identify();
0268     }
0269
0270     set_vtx(vtx->get_x(), vtx->get_y(), vtx->get_z());
0271     return 1;
0272
0273   }  // if (_ineve->GetNVtx() > 0) {
0274   // try the next option for getting primary vertex
0275
0276   PHG4TruthInfoContainer *truthInfoList =  //
0277       findNode::getClass<PHG4TruthInfoContainer>(topNode,
0278                                                  "G4TruthInfo");
0279   if (truthInfoList)
0280   {
0281     // embed to vertexes as set by primary vertex from the truth container, e.g. when embedding into DST
0282
0283     PHG4VtxPoint *vtx = truthInfoList->GetPrimaryVtx(1);
0284
0285     if (!vtx)
0286     {
0287       truthInfoList->identify();
0288       std::cout << PHWHERE << "::Error - PHG4SimpleEventGenerator expects an existing truth vertex in PHG4TruthInfoContainer, but none exists"
0289                 << std::endl;
0290       exit(1);
0291     }
0292
0293     if (Verbosity() > 0)
0294     {
0295       std::cout << PHWHERE << "::Info - use this primary vertex from PHG4TruthInfoContainer:" << std::endl;
0296       vtx->identify();
0297     }
0298
0299     set_vtx(vtx->get_x(), vtx->get_y(), vtx->get_z());
0300     return 1;
0301   }
0302
0303   // I am out of options.....
0304
0305   std::cout << PHWHERE << "::Error we expect an existing truth vertex, but none exists"
0306             << std::endl;
0307   gSystem->Exit(1);
0308
0309   return 0;
0310 }
0311
0312 void PHG4ParticleGeneratorBase::ResetParticleList()
0313 {
0314   while (particlelist_begin() != particlelist_end())
0315   {
0316     delete particlelist.back();
0317     particlelist.pop_back();
0318   }
0319   return;
0320 }