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