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/PhysicalConstants.h>
0022 #include <CLHEP/Units/SystemOfUnits.h>
0023
0024 // Geant4 include files
0025 #include <G4ParticleTable.hh>
0026 #include <G4ParticleDefinition.hh>
0027
0028 // C/C++ include files
0029 #include <stdexcept>
0030 #include <cmath>
0031
0032 using namespace dd4hep::sim;
0033
0034 /// Standard constructor
0035 Geant4ParticleGenerator::Geant4ParticleGenerator(Geant4Context* ctxt, const std::string& nam)
0036   : Geant4GeneratorAction(ctxt, nam), m_direction(0,1,0), m_position(0,0,0), m_particle(0)
0037 {
0038   InstanceCount::increment(this);
0039   m_needsControl = true;
0040   declareProperty("Particle",      m_particleName = "e-");
0041   declareProperty("Energy",        m_energy = -1);
0042   declareProperty("energy",        m_energy = -1);
0043   declareProperty("MomentumMin",   m_momentumMin = 0.0);
0044   declareProperty("MomentumMax",   m_momentumMax = 50 * CLHEP::MeV);
0045   declareProperty("Multiplicity",  m_multiplicity = 1);
0046   declareProperty("Mask",          m_mask = 0);
0047   declareProperty("Position",      m_position = ROOT::Math::XYZVector(0.,0.,0.));
0048   declareProperty("Direction",     m_direction = ROOT::Math::XYZVector(1.,1.,1.));
0049 }
0050
0051 /// Default destructor
0052 Geant4ParticleGenerator::~Geant4ParticleGenerator() {
0053   InstanceCount::decrement(this);
0054 }
0055
0056 /// Particle modification. Caller presets defaults to: ( direction = m_direction, momentum = [mMin, mMax])
0057 void Geant4ParticleGenerator::getParticleDirection(int , ROOT::Math::XYZVector& , double& momentum) const {
0058   getParticleMomentumUniform(momentum);
0059 }
0060
0061 /// Uniform momentum distribution
0062 void Geant4ParticleGenerator::getParticleMomentumUniform(double& momentum) const {
0063   if (m_energy != -1) {
0064     momentum = m_energy;
0065     return;
0066   }
0067   Geant4Event&  evt = context()->event();
0068   Geant4Random& rnd = evt.random();
0069   if (m_momentumMax < m_momentumMin)
0070     // We no longer set m_momentumMax to -1 so, not entirely sure a) this will still happen, b) actually work
0071     momentum = m_momentumMin+(momentum-m_momentumMin)*rnd.rndm();
0072   else
0073     momentum = m_momentumMin+(m_momentumMax-m_momentumMin)*rnd.rndm();
0074 }
0075
0076 /// Particle modification. Caller presets defaults to: (multiplicity=m_multiplicity)
0077 void Geant4ParticleGenerator::getParticleMultiplicity(int& ) const   {
0078 }
0079
0080 /// Particle modification. Caller presets defaults to: (multiplicity=m_multiplicity)
0081 void Geant4ParticleGenerator::getVertexPosition(ROOT::Math::XYZVector& ) const   {
0082 }
0083
0084 /// Print single particle interaction identified by its mask
0085 void Geant4ParticleGenerator::printInteraction(int mask)  const  {
0086   Geant4PrimaryEvent* prim = context()->event().extension<Geant4PrimaryEvent>();
0087   if ( !prim )   {
0088     warning("printInteraction: Bad primary event [NULL-Pointer].");
0089     return;
0090   }
0091   Geant4PrimaryInteraction* inter = prim->get(mask);
0092   if ( !inter )   {
0093     warning("printInteraction: Bad interaction identifier 0x%08X [Unknown Mask].",mask);
0094     return;
0095   }
0096   printInteraction(inter);
0097 }
0098
0099 /// Print single particle interaction identified by its reference
0100 void Geant4ParticleGenerator::printInteraction(Geant4PrimaryInteraction* inter)  const  {
0101   int count = 0;
0102   if ( !inter )   {
0103     warning("printInteraction: Invalid interaction pointer [NULL-Pointer].");
0104     return;
0105   }
0106   for(const auto& iv : inter->vertices )   {
0107     for( Geant4Vertex* v : iv.second ){
0108       print("+-> Interaction [%d] [%.3f , %.3f] GeV %s pos:(%.3f %.3f %.3f)[mm]",
0109             count, m_momentumMin/CLHEP::GeV, m_momentumMax/CLHEP::GeV, m_particleName.c_str(),
0110             v->x/CLHEP::mm, v->y/CLHEP::mm, v->z/CLHEP::mm);
0111       ++count;
0112       for ( int i : v->out )  {
0113         Geant4ParticleHandle p = inter->particles[i];
0114         p.dumpWithVertex(outputLevel(),name(),"  +->");
0115       }
0116     }
0117   }
0118 }
0119
0120 /// Callback to generate primary particles
0121 void Geant4ParticleGenerator::operator()(G4Event*) {
0122   typedef Geant4Particle Particle;
0123
0124   if (0 == m_particle || m_particle->GetParticleName() != m_particleName.c_str()) {
0125     G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0126     m_particle = particleTable->FindParticle(m_particleName);
0127     if (0 == m_particle) {
0128       except("Geant4ParticleGenerator: Bad particle type: %s!", m_particleName.c_str());
0129     }
0130   }
0131   Geant4Event& evt = context()->event();
0132   Geant4PrimaryEvent* prim = evt.extension<Geant4PrimaryEvent>();
0133   Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
0134   prim->add(m_mask, inter);
0135
0136   Geant4Vertex* vtx = new Geant4Vertex();
0137   int multiplicity = m_multiplicity;
0138   ROOT::Math::XYZVector unit_direction, position = m_position;
0139   getVertexPosition(position);
0140   getParticleMultiplicity(multiplicity);
0141   vtx->mask = m_mask;
0142   vtx->x = position.X();
0143   vtx->y = position.Y();
0144   vtx->z = position.Z();
0145   inter->vertices[m_mask].emplace_back( vtx );
0146   for(int i=0; i<m_multiplicity; ++i)   {
0147     double momentum = 0.0;
0148     ROOT::Math::XYZVector direction = m_direction;
0149     Particle* p = new Particle();
0150     getParticleDirection(i, direction, momentum);
0151     unit_direction  = direction.unit();
0152     p->id           = inter->nextPID();
0153     p->status      |= G4PARTICLE_GEN_STABLE;
0154     p->mask         = m_mask;
0155     p->pdgID        = m_particle->GetPDGEncoding();
0156
0157     p->psx          = unit_direction.X()*momentum;
0158     p->psy          = unit_direction.Y()*momentum;
0159     p->psz          = unit_direction.Z()*momentum;
0160     p->mass         = m_particle->GetPDGMass();
0161     p->charge       = 3 * m_particle->GetPDGCharge();
0162     // Optical photons must have a defined polarization; assign a random one
0163     // perpendicular to the momentum direction (physical requirement for Geant4 optics).
0164     if ( m_particle->GetParticleName() == "opticalphoton" )  {
0165       if ( momentum == 0.0 )  {
0166         except("Geant4ParticleGenerator: Cannot assign polarization to optical photon with zero momentum!");
0167       }
0168       // Build two unit vectors orthogonal to the photon direction using Gram-Schmidt.
0169       ROOT::Math::XYZVector e1, e2;
0170       if ( std::fabs(unit_direction.X()) < 0.9 )
0171         e1 = ROOT::Math::XYZVector(1, 0, 0);
0172       else
0173         e1 = ROOT::Math::XYZVector(0, 1, 0);
0174       e2  = unit_direction.Cross(e1).unit();
0175       e1  = e2.Cross(unit_direction).unit();    // normalize defensively
0176       Geant4Random& rnd = evt.random();
0177       double angle = CLHEP::twopi * rnd.rndm();
0178       ROOT::Math::XYZVector pol = (std::cos(angle)*e1 + std::sin(angle)*e2).unit();
0179       p->spin[0] = pol.X();
0180       p->spin[1] = pol.Y();
0181       p->spin[2] = pol.Z();
0182     }
0183     else  {
0184       p->spin[0] = 0;
0185       p->spin[1] = 0;
0186       p->spin[2] = 0;
0187     }
0188     p->colorFlow[0] = 0;
0189     p->colorFlow[1] = 0;
0190     p->vsx        = vtx->x;
0191     p->vsy        = vtx->y;
0192     p->vsz        = vtx->z;
0193     //fg: do not set the endpoint to the start point of the particle
0194     // p->vex        = vtx->x;
0195     // p->vey        = vtx->y;
0196     // p->vez        = vtx->z;
0197     inter->particles.emplace(p->id,p);
0198     vtx->out.insert(p->id);
0199     printout(INFO,name(),"Particle [%d] %-12s Mom:%.3f GeV vertex:(%6.3f %6.3f %6.3f)[mm] direction:(%6.3f %6.3f %6.3f)",
0200              p->id, m_particleName.c_str(), momentum/CLHEP::GeV,
0201              vtx->x/CLHEP::mm, vtx->y/CLHEP::mm, vtx->z/CLHEP::mm,
0202              direction.X(), direction.Y(), direction.Z());
0203   }
0204 }