DDG4/src/Geant4ParticleGenerator.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 // Framework include files
0015 #include <DD4hep/Printout.h>
0016 #include <DD4hep/InstanceCount.h>
0017 #include <DDG4/Geant4Context.h>
0018 #include <DDG4/Geant4Primary.h>
0019 #include <DDG4/Geant4ParticleGenerator.h>
0020 #include <DDG4/Geant4Random.h>
0021 #include <CLHEP/Units/SystemOfUnits.h>
0022
0023 // Geant4 include files
0024 #include <G4ParticleTable.hh>
0025 #include <G4ParticleDefinition.hh>
0026
0027 // C/C++ include files
0028 #include <stdexcept>
0029 #include <cmath>
0030
0031 using namespace dd4hep::sim;
0032
0033 /// Standard constructor
0034 Geant4ParticleGenerator::Geant4ParticleGenerator(Geant4Context* ctxt, const std::string& nam)
0035   : Geant4GeneratorAction(ctxt, nam), m_direction(0,1,0), m_position(0,0,0), m_particle(0)
0036 {
0037   InstanceCount::increment(this);
0038   m_needsControl = true;
0039   declareProperty("Particle",      m_particleName = "e-");
0040   declareProperty("Energy",        m_energy = -1);
0041   declareProperty("energy",        m_energy = -1);
0042   declareProperty("MomentumMin",   m_momentumMin = 0.0);
0043   declareProperty("MomentumMax",   m_momentumMax = 50 * CLHEP::MeV);
0044   declareProperty("Multiplicity",  m_multiplicity = 1);
0045   declareProperty("Mask",          m_mask = 0);
0046   declareProperty("Position",      m_position = ROOT::Math::XYZVector(0.,0.,0.));
0047   declareProperty("Direction",     m_direction = ROOT::Math::XYZVector(1.,1.,1.));
0048 }
0049
0050 /// Default destructor
0051 Geant4ParticleGenerator::~Geant4ParticleGenerator() {
0052   InstanceCount::decrement(this);
0053 }
0054
0055 /// Particle modification. Caller presets defaults to: ( direction = m_direction, momentum = [mMin, mMax])
0056 void Geant4ParticleGenerator::getParticleDirection(int , ROOT::Math::XYZVector& , double& momentum) const {
0057   getParticleMomentumUniform(momentum);
0058 }
0059
0060 /// Uniform momentum distribution
0061 void Geant4ParticleGenerator::getParticleMomentumUniform(double& momentum) const {
0062   if (m_energy != -1) {
0063     momentum = m_energy;
0064     return;
0065   }
0066   Geant4Event&  evt = context()->event();
0067   Geant4Random& rnd = evt.random();
0068   if (m_momentumMax < m_momentumMin)
0069     // We no longer set m_momentumMax to -1 so, not entirely sure a) this will still happen, b) actually work
0070     momentum = m_momentumMin+(momentum-m_momentumMin)*rnd.rndm();
0071   else
0072     momentum = m_momentumMin+(m_momentumMax-m_momentumMin)*rnd.rndm();
0073 }
0074
0075 /// Particle modification. Caller presets defaults to: (multiplicity=m_multiplicity)
0076 void Geant4ParticleGenerator::getParticleMultiplicity(int& ) const   {
0077 }
0078
0079 /// Particle modification. Caller presets defaults to: (multiplicity=m_multiplicity)
0080 void Geant4ParticleGenerator::getVertexPosition(ROOT::Math::XYZVector& ) const   {
0081 }
0082
0083 /// Print single particle interaction identified by its mask
0084 void Geant4ParticleGenerator::printInteraction(int mask)  const  {
0085   Geant4PrimaryEvent* prim = context()->event().extension<Geant4PrimaryEvent>();
0086   if ( !prim )   {
0087     warning("printInteraction: Bad primary event [NULL-Pointer].");
0088     return;
0089   }
0090   Geant4PrimaryInteraction* inter = prim->get(mask);
0091   if ( !inter )   {
0092     warning("printInteraction: Bad interaction identifier 0x%08X [Unknown Mask].",mask);
0093     return;
0094   }
0095   printInteraction(inter);
0096 }
0097
0098 /// Print single particle interaction identified by its reference
0099 void Geant4ParticleGenerator::printInteraction(Geant4PrimaryInteraction* inter)  const  {
0100   int count = 0;
0101   if ( !inter )   {
0102     warning("printInteraction: Invalid interaction pointer [NULL-Pointer].");
0103     return;
0104   }
0105   for(const auto& iv : inter->vertices )   {
0106     for( Geant4Vertex* v : iv.second ){
0107       print("+-> Interaction [%d] [%.3f , %.3f] GeV %s pos:(%.3f %.3f %.3f)[mm]",
0108             count, m_momentumMin/CLHEP::GeV, m_momentumMax/CLHEP::GeV, m_particleName.c_str(),
0109             v->x/CLHEP::mm, v->y/CLHEP::mm, v->z/CLHEP::mm);
0110       ++count;
0111       for ( int i : v->out )  {
0112         Geant4ParticleHandle p = inter->particles[i];
0113         p.dumpWithVertex(outputLevel(),name(),"  +->");
0114       }
0115     }
0116   }
0117 }
0118
0119 /// Callback to generate primary particles
0120 void Geant4ParticleGenerator::operator()(G4Event*) {
0121   typedef Geant4Particle Particle;
0122
0123   if (0 == m_particle || m_particle->GetParticleName() != m_particleName.c_str()) {
0124     G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0125     m_particle = particleTable->FindParticle(m_particleName);
0126     if (0 == m_particle) {
0127       except("Geant4ParticleGenerator: Bad particle type: %s!", m_particleName.c_str());
0128     }
0129   }
0130   Geant4Event& evt = context()->event();
0131   Geant4PrimaryEvent* prim = evt.extension<Geant4PrimaryEvent>();
0132   Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
0133   prim->add(m_mask, inter);
0134
0135   Geant4Vertex* vtx = new Geant4Vertex();
0136   int multiplicity = m_multiplicity;
0137   ROOT::Math::XYZVector unit_direction, position = m_position;
0138   getVertexPosition(position);
0139   getParticleMultiplicity(multiplicity);
0140   vtx->mask = m_mask;
0141   vtx->x = position.X();
0142   vtx->y = position.Y();
0143   vtx->z = position.Z();
0144   inter->vertices[m_mask].emplace_back( vtx );
0145   for(int i=0; i<m_multiplicity; ++i)   {
0146     double momentum = 0.0;
0147     ROOT::Math::XYZVector direction = m_direction;
0148     Particle* p = new Particle();
0149     getParticleDirection(i, direction, momentum);
0150     unit_direction  = direction.unit();
0151     p->id           = inter->nextPID();
0152     p->status      |= G4PARTICLE_GEN_STABLE;
0153     p->mask         = m_mask;
0154     p->pdgID        = m_particle->GetPDGEncoding();
0155
0156     p->psx          = unit_direction.X()*momentum;
0157     p->psy          = unit_direction.Y()*momentum;
0158     p->psz          = unit_direction.Z()*momentum;
0159     p->mass         = m_particle->GetPDGMass();
0160     p->charge       = 3 * m_particle->GetPDGCharge();
0161     p->spin[0]      = 0;
0162     p->spin[1]      = 0;
0163     p->spin[2]      = 0;
0164     p->colorFlow[0] = 0;
0165     p->colorFlow[1] = 0;
0166     p->vsx        = vtx->x;
0167     p->vsy        = vtx->y;
0168     p->vsz        = vtx->z;
0169     //fg: do not set the endpoint to the start point of the particle
0170     // p->vex        = vtx->x;
0171     // p->vey        = vtx->y;
0172     // p->vez        = vtx->z;
0173     inter->particles.emplace(p->id,p);
0174     vtx->out.insert(p->id);
0175     printout(INFO,name(),"Particle [%d] %-12s Mom:%.3f GeV vertex:(%6.3f %6.3f %6.3f)[mm] direction:(%6.3f %6.3f %6.3f)",
0176              p->id, m_particleName.c_str(), momentum/CLHEP::GeV,
0177          vtx->x/CLHEP::mm, vtx->y/CLHEP::mm, vtx->z/CLHEP::mm,
0178          direction.X(), direction.Y(), direction.Z());
0179   }
0180 }