Geant4InputHandling.cpp

Go to the documentation of this file.

//==========================================================================
//  AIDA Detector description implementation 
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author     : M.Frank
//
//==========================================================================
 
// Framework include files
#include <DD4hep/Printout.h>
#include <DDG4/Geant4InputHandling.h>
#include <DDG4/Geant4Primary.h>
#include <DDG4/Geant4Context.h>
#include <DDG4/Geant4Action.h>
#include <DDG4/Geant4PrimaryHandler.h>
#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
 
// Geant4 include files
#include <G4Event.hh>
#include <G4PrimaryVertex.hh>
#include <G4PrimaryParticle.hh>
#include <G4ParticleDefinition.hh>
#include <G4VPhysicalVolume.hh>
#include <G4LogicalVolume.hh>
#include <G4VSolid.hh>
 
// C/C++ include files
#include <stdexcept>
#include <limits>
#include <cmath>
 
using namespace dd4hep::sim;
using PropertyMask = dd4hep::detail::ReferenceBitMask<int>;
 
Geant4Vertex* dd4hep::sim::createPrimary(const G4PrimaryVertex* g4)      {        
  Geant4Vertex* v = new Geant4Vertex();
  v->x = g4->GetX0();
  v->y = g4->GetY0();
  v->z = g4->GetZ0();
  v->time = g4->GetT0();
  return v;
}
 
Geant4Particle* 
dd4hep::sim::createPrimary(int particle_id, 
                           const Geant4Vertex* v,
                           const G4PrimaryParticle* g4p)
{
  Geant4ParticleHandle p = new Geant4Particle();
  p->id           = particle_id;
  p->reason       = 0;
  p->pdgID        = g4p->GetPDGcode();
  p->psx          = g4p->GetPx();
  p->psy          = g4p->GetPy();
  p->psz          = g4p->GetPz();
  p->time         = g4p->GetProperTime();
  p->properTime   = g4p->GetProperTime();
  p->vsx          = v->x;
  p->vsy          = v->y;
  p->vsz          = v->z;
  p->vex          = v->x;
  p->vey          = v->y;
  p->vez          = v->z;
  //p->definition   = g4p->GetG4code();
  //p->process      = 0;
  p->spin[0]      = 0;
  p->spin[1]      = 0;
  p->spin[2]      = 0;
  p->colorFlow[0] = 0;
  p->colorFlow[0] = 0;
  p->mass         = g4p->GetMass();
  p->charge       = int(3.0 * g4p->GetCharge());
  PropertyMask status(p->status);
  status.set(G4PARTICLE_GEN_STABLE);
  return p;
}
 
static void collectPrimaries(Geant4PrimaryMap*         pm,
                             Geant4PrimaryInteraction* interaction, 
                             Geant4Vertex*             particle_origine,
                             G4PrimaryParticle*        gp)
{
  //if the particle is in the map, we do not have to do anything
  if ( pm->get(gp) )  {
    return;
  }
 
  int pid = int(interaction->particles.size());
  Geant4Particle* p = createPrimary(pid,particle_origine,gp);
  G4PrimaryParticle* dau = gp->GetDaughter();
  PropertyMask status(p->status);
  int mask = interaction->mask;
  
  interaction->particles.emplace(p->id,p);
  status.set(G4PARTICLE_PRIMARY);
  p->mask = mask;
  particle_origine->out.insert(p->id);
  // Insert pair in map. Do not forget to increase reference count!
  pm->insert(gp,p);
 
  if ( dau )   {
    Geant4Vertex* dv = new Geant4Vertex(*particle_origine);
    PropertyMask reason(p->reason);
    reason.set(G4PARTICLE_HAS_SECONDARIES);
 
    dv->mask = mask;
    dv->in.insert(p->id);
 
    interaction->vertices[mask].emplace_back(dv) ;
 
    for(; dau; dau = dau->GetNext())
      collectPrimaries(pm, interaction, dv, dau);
  }
}
 
Geant4PrimaryInteraction* 
dd4hep::sim::createPrimary(int mask,
                           Geant4PrimaryMap* pm,
                           std::set<G4PrimaryVertex*>const& primaries)
{
  Geant4PrimaryInteraction* interaction = new Geant4PrimaryInteraction();
  interaction->locked = true;
  interaction->mask = mask;
  for (auto const& gv: primaries) {
    Geant4Vertex* v = createPrimary(gv);
    v->mask = mask;
    interaction->vertices[mask].emplace_back(v);
    for (G4PrimaryParticle *gp = gv->GetPrimary(); gp; gp = gp->GetNext()) {
      collectPrimaries(pm, interaction, v, gp);
    }
  }
  return interaction;
}
 
int dd4hep::sim::generationInitialization(const Geant4Action* /* caller */,
                                          const Geant4Context* context)
{
  context->event().addExtension(new Geant4PrimaryMap());
 
  // The final set of created particles in the simulation.
  context->event().addExtension(new Geant4ParticleMap());
 
  //
  // The Geant4PrimaryEvent extension contains a whole set of
  // Geant4PrimaryInteraction objects each may represent a complete
  // interaction. Particles and vertices may be unbiased.
  // This is the input to the translator forming the final
  // Geant4PrimaryInteraction (see below) containing rebiased particle
  // and vertex maps.
  Geant4PrimaryEvent* evt = new Geant4PrimaryEvent();
  context->event().addExtension(evt);
  //
  // The Geant4PrimaryInteraction extension contains the final
  // vertices and particles ready to be injected to Geant4.
  Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
  inter->setNextPID(-1);
  context->event().addExtension(inter);
  return 1;
}
 
static void appendInteraction(const Geant4Action* caller,
                              Geant4PrimaryInteraction* output,
                              Geant4PrimaryInteraction* input)
{
  Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
  for( ip=input->particles.begin(), ipend=input->particles.end(); ip != ipend; ++ip )  {
    Geant4Particle* p = (*ip).second;
    output->particles.emplace(p->id,p->addRef());
  }
  Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
  for( iv=input->vertices.begin(), ivend=input->vertices.end(); iv != ivend; ++iv )   {
    int theMask = input->mask;
    ivfnd = output->vertices.find(theMask);
    if ( ivfnd != output->vertices.end() )   {
      caller->abortRun("Duplicate primary interaction identifier!",
                       "Cannot handle 2 interactions with identical identifiers!");
    }
    for(Geant4Vertex* vtx :  (*iv).second )
      output->vertices[theMask].emplace_back( vtx->addRef() );
  }
}
 
static void rebaseParticles(Geant4PrimaryInteraction::ParticleMap& particles, int &offset)    {
  Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
  int mx_id = offset;
  // Now move begin and end-vertex of all primary and generator particles accordingly
  for( ip=particles.begin(), ipend=particles.end(); ip != ipend; ++ip )  {
    Geant4ParticleHandle p((*ip).second);
    p.offset(offset);
    mx_id = p->id+1 > mx_id ? p->id+1 : mx_id;
  }
  offset = mx_id;
}
 
static void rebaseVertices(Geant4PrimaryInteraction::VertexMap& vertices, int part_offset)    {
  Geant4PrimaryInteraction::VertexMap::iterator iv, ivend;
  std::set<int> in, out;
  std::set<int>::iterator i;
  // Now move begin and end-vertex of all primary vertices accordingly
  for(iv=vertices.begin(), ivend=vertices.end(); iv != ivend; ++iv)  {
    for( Geant4Vertex* v : (*iv).second ){ 
      in  = v->in;
      out = v->out;
      for(in=v->in, i=in.begin(), v->in.clear(); i != in.end(); ++i)
        v->in.insert((*i)+part_offset);
      for(out=v->out, i=out.begin(), v->out.clear(); i != out.end(); ++i)
        v->out.insert((*i)+part_offset);
    }
  }
}
int dd4hep::sim::mergeInteractions(const Geant4Action* caller,
                                   const Geant4Context* context)
{
  typedef Geant4PrimaryEvent::Interaction  Interaction;
  typedef std::vector<Interaction*>        Interactions;
  Geant4Event& event = context->event();
  Geant4PrimaryEvent* evt = event.extension<Geant4PrimaryEvent>();
  Interaction* output = event.extension<Interaction>();
  Interactions inter  = evt->interactions();
  int particle_offset = 0;
 
  for(Interactions::const_iterator i=inter.begin(); i != inter.end(); ++i)  {
    Interaction* interaction = *i;
    int vertex_offset = particle_offset;
    if ( !interaction->applyMask() )   {
      caller->abortRun("Found single interaction with multiple primary vertices!",
                       "Cannot merge individual interactions with more than one primary!");
    }
    rebaseParticles(interaction->particles,particle_offset);
    rebaseVertices(interaction->vertices,vertex_offset);
    appendInteraction(caller,output,interaction);
  }
  output->setNextPID(particle_offset);
  Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
  caller->debug("+++ Merging MC input record from %d interactions:",(int)inter.size());
  for( ip=output->particles.begin(), ipend=output->particles.end(); ip != ipend; ++ip )
    Geant4ParticleHandle((*ip).second).dump1(DEBUG,caller->name(),"Merged particles");
  return 1;
}
 
int dd4hep::sim::boostInteraction(const Geant4Action* caller,
                                  Geant4PrimaryEvent::Interaction* inter,
                                  double alpha)
{
#define SQR(x)  (x*x)
  Geant4PrimaryEvent::Interaction::VertexMap::iterator iv;
  Geant4PrimaryEvent::Interaction::ParticleMap::iterator ip;
  double gamma = std::sqrt(1 + SQR(tan(alpha)));
  double betagamma = std::tan(alpha);
 
  if ( inter->locked )  {
    caller->abortRun("Locked interactions may not be boosted!",
                     "Cannot boost interactions with a native G4 primary record!");
  }
  else if ( alpha != 0.0 )  {
    // Now move begin and end-vertex of all primary vertices accordingly
    for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv){  
      for( Geant4Vertex* v :  (*iv).second ) {
        double t = gamma * v->time + betagamma * v->x / CLHEP::c_light;
        double x = gamma * v->x + betagamma * CLHEP::c_light * v->time;
        double y = v->y;
        double z = v->z;
        v->x = x;
        v->y = y;
        v->z = z;
        v->time = t;
      }
    }
    // Now move begin and end-vertex of all primary and generator particles accordingly
    for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip)  {
      Geant4ParticleHandle p = (*ip).second;
      double t = gamma * p->time + betagamma * p->vsx / CLHEP::c_light;
      double x = gamma * p->vsx + betagamma * CLHEP::c_light * p->time;
      double y = p->vsy;
      double z = p->vsz;
      
      double m  = p->mass;
      double e2 = SQR(p->psx)+SQR(p->psy)+SQR(p->psz)+SQR(m);
      double px = betagamma * std::sqrt(e2) + gamma * p->psx;
      double py = p->psy;
      double pz = p->psz;
      
      p->vsx = x;
      p->vsy = y;
      p->vsz = z;
      p->time = t;
      
      p->psx = px;
      p->psy = py;
      p->psz = pz;
    }
  }
  return 1;
}
 
int dd4hep::sim::smearInteraction(const Geant4Action* caller,
                                  Geant4PrimaryEvent::Interaction* inter,
                                  double dx, double dy, double dz, double dt)
{
  Geant4PrimaryEvent::Interaction::VertexMap::iterator iv;
  Geant4PrimaryEvent::Interaction::ParticleMap::iterator ip;
 
  if ( inter->locked )  {
    caller->abortRun("Locked interactions may not be smeared!",
                     "Cannot smear interactions with a native G4 primary record!");
  }
  
  // Now move begin and end-vertex of all primary vertices accordingly
  for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv)  {
    for( Geant4Vertex* v : (*iv).second ){
      v->x += dx;
      v->y += dy;
      v->z += dz;
      v->time += dt;
    }
  }
  // Now move begin and end-vertex of all primary and generator particles accordingly
  for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip)  {
    Geant4Particle* p = (*ip).second;
    p->vsx += dx;
    p->vsy += dy;
    p->vsz += dz;
    p->vex += dx;
    p->vey += dy;
    p->vez += dz;
    p->time += dt;
  }
  return 1;
}
 
static G4PrimaryParticle* createG4Primary(const Geant4ParticleHandle p)  {
  G4PrimaryParticle* g4 = 0;
  const G4ParticleDefinition* def = p.definition();
  double energy = p.energy();
  double mom2   = (p->psx*p->psx) + (p->psy*p->psy) + (p->psz*p->psz);
  double mass2  = energy*energy - mom2;
  if ( mass2 < 0e0 )   {
    if ( def )   {
      mass2 = def->GetPDGMass() * def->GetPDGMass();
    }
    energy = std::sqrt(mom2 + mass2);
    if ( std::fabs(p.energy()-energy) > 0e0 /* 1e-10 */ )  {
      dd4hep::printout(dd4hep::INFO,"createG4Primary",
                       "Change particle %s energy from %10.5f MeV by %g ppm to avoid negative Energy^2",
                       (def) ? def->GetParticleName().c_str() : "???", p.energy(), std::fabs(p.energy()-energy)*1e6);
    }
  }
  if ( 0 != p->pdgID )   {
    // For ions we use the pdgID of the definition, in case we had to zero the excitation level, see Geant4Particle.cpp
    const int pdgID =  p->pdgID < 1000000000 ? p->pdgID : p.definition()->GetPDGEncoding();
    g4 = new G4PrimaryParticle(pdgID, p->psx, p->psy, p->psz, energy);
  }
  else   {
    g4 = new G4PrimaryParticle(def, p->psx, p->psy, p->psz, energy);
    g4->SetCharge(double(p.charge())/3.0);
  }
  // The particle is fully defined with the 4-vector set above, setting the mass isn't necessary, not
  // using the 4-vector, means the PDG mass is used, and the momentum is scaled if the mass is set here
  // g4->SetMass(p->mass);
  if ( p->spin[0] != 0.0 || p->spin[1] != 0.0 || p->spin[2] != 0.0 )  {
    g4->SetPolarization(p->spin[0], p->spin[1], p->spin[2]);
  }
  return g4;
}
 
static std::vector< std::pair<Geant4Particle*,G4PrimaryParticle*> >
getRelevant(std::set<int>& visited,
            std::map<int,G4PrimaryParticle*>& prim,
            Geant4PrimaryInteraction::ParticleMap& pm,
            const Geant4PrimaryConfig& primaryConfig,
            const Geant4ParticleHandle p)
{
  typedef std::vector< std::pair<Geant4Particle*,G4PrimaryParticle*> > Primaries;
  using dd4hep::printout;
  
  Primaries res;
  visited.insert(p->id);
  PropertyMask status(p->status);
  if ( status.isSet(G4PARTICLE_GEN_STABLE) )  {
    bool rejectParticle = false
      or (primaryConfig.m_rejectPDGs.count(abs(p->pdgID)) != 0) // quarks, gluon, "strings", W, Z etc.
      ;
    printout(dd4hep::DEBUG, "Input",
             "Checking rejection of stable: PDG(%-10d), Definition(%s), reject(%s)",
             p->pdgID,
             p.definition() ? "true" : "false",
             rejectParticle ? "true" : "false");
    if (not rejectParticle and prim.find(p->id) == prim.end() )  {
      G4PrimaryParticle* p4 = createG4Primary(p);
      prim[p->id] = p4;
      res.emplace_back(p,p4);
    }
  }
  else if ( p->daughters.size() > 0 )  {
    const Geant4Particle::Particles& dau = p->daughters;
    int first_daughter = *(dau.begin());
    Geant4ParticleHandle dp = pm[first_daughter];
    double en = p.energy();
    double me = en > std::numeric_limits<double>::epsilon() ? p->mass / en : 0.0;
    //  fix by S.Morozov for real != 0
    double proper_time = fabs(dp->time-p->time) * me;
    double proper_time_Precision = pow(10.,-DBL_DIG)*fabs(me)*fmax(fabs(p->time),fabs(dp->time));
    bool isProperTimeZero = (fabs(proper_time) <= fabs(proper_time_Precision));
 
    // -- remove original if ---
    bool rejectParticle = not p.definition()                    // completely unknown to geant4
      or (primaryConfig.m_rejectPDGs.count(abs(p->pdgID)) != 0) // quarks, gluon, "strings", W, Z etc.
      or (isProperTimeZero and p.definition()->GetPDGStable() ) // initial state electrons, etc.
      or (isProperTimeZero and primaryConfig.m_zeroTimePDGs.count(abs(p->pdgID)) != 0 )  // charged 'documentation' leptons, e.g. in lepton pairs w/ FSR
      or (status.isSet(G4PARTICLE_GEN_DOCUMENTATION) || status.isSet(G4PARTICLE_GEN_BEAM) || status.isSet(G4PARTICLE_GEN_OTHER))  // documentation generator status
      or false;
 
    printout(dd4hep::DEBUG, "Input",
             "Checking rejection: PDG(%-10d), Definition(%s), isProperTimeZero(%s, %3.15f), stable(%s), doc(%s), reject(%s)",
             p->pdgID,
             p.definition() ? "true" : "false",
             isProperTimeZero ? "true" : "false", proper_time,
             (bool(p.definition()) ? p.definition()->GetPDGStable() : false)  ? "true" : "false",
             status.isSet(G4PARTICLE_GEN_DOCUMENTATION) || status.isSet(G4PARTICLE_GEN_BEAM) || status.isSet(G4PARTICLE_GEN_OTHER) ? "true" : "false",
             rejectParticle ? "true" : "false");
    // end running simulation if we have a really inconsistent record, that is unrejected stable particle with children
    bool failStableWithChildren = (not rejectParticle and p.definition()->GetPDGStable());
    if (failStableWithChildren) {
      printout(dd4hep::FATAL,"Input",
               "+++ Stable particle (PDG: %-10d) with daughters! check your MC record, adapt particle.tbl file...",
               p->pdgID);
      throw std::runtime_error("Cannot Simmulate this MC Record");
    }
    if (not rejectParticle) {
      std::map<int, G4PrimaryParticle*>::iterator ip4 = prim.find(p->id);
      G4PrimaryParticle* p4 = (ip4 == prim.end()) ? 0 : (*ip4).second;
      if ( !p4 )  {
        p4 = createG4Primary(p);
        // if the user wants the particle with this PDG id to be decayed according to the lifetime distrution configured
        // in particle.tbl (or the geant4 defaults) then they have to configure this. This is needed since 0.0 is now a
        // allowed pre-defined decay time by geant4
        if(primaryConfig.m_decayByGeant.count(abs(p->pdgID))) {
          p4->SetProperTime(-1);
        } else {
          p4->SetProperTime(proper_time);
        }
        prim[p->id] = p4;
        Primaries daughters;
        for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i)  {
          if ( visited.find(*i) == visited.end() )  {
            Primaries tmp = getRelevant(visited,prim,pm,primaryConfig,pm[*i]);
            daughters.insert(daughters.end(), tmp.begin(),tmp.end());
          }
        }
        for(Primaries::iterator i=daughters.begin(); i!=daughters.end(); ++i)
          p4->SetDaughter((*i).second);
      }
      res.emplace_back(p,p4);
    }
    else  {
      for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i)  {
        if ( visited.find(*i) == visited.end() )  {
          Primaries tmp = getRelevant(visited,prim,pm,primaryConfig,pm[*i]);
          res.insert(res.end(), tmp.begin(),tmp.end());
        }
      }
    }
  }
  return res;
}
 
namespace {
 
  bool insideWorldVolume(const G4VPhysicalVolume* world, const G4ThreeVector& point)   {
    auto* logical   = world ? world->GetLogicalVolume() : nullptr;
    G4VSolid* solid = logical ? logical->GetSolid() : nullptr;
    return solid ? (solid->Inside(point) != kOutside) : true;
  }
}
 
int dd4hep::sim::generatePrimaries(const Geant4Action* caller,
                                   const Geant4Context* context,
                                   G4Event* event)
{
  typedef std::vector< std::pair<Geant4Particle*,G4PrimaryParticle*> > Primaries;
  typedef Geant4PrimaryInteraction Interaction;
  Geant4PrimaryMap* primaries   = context->event().extension<Geant4PrimaryMap>();
  Interaction*      interaction = context->event().extension<Interaction>();
  Interaction::ParticleMap& pm  = interaction->particles;
  Interaction::VertexMap&   vm  = interaction->vertices;
  std::map<int,G4PrimaryParticle*> prim;
  std::set<int> visited;
 
  auto const* primHandler = dynamic_cast<const Geant4PrimaryHandler*>(caller);
  auto const& primaryConfig = primHandler ? primHandler->m_primaryConfig : Geant4PrimaryConfig();
 
  caller->debug("PrimaryConfiguration:%s", primaryConfig.toString().c_str());
 
  if ( interaction->locked )  {
    caller->abortRun("Locked interactions may not be used to generate primaries!",
                     "Cannot handle a native G4 primary record!");
    return 0;
  }
  else   {
    Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
    for(Interaction::VertexMap::const_iterator iend=vm.end(),i=vm.begin(); i!=iend; ++i)  {
      for( Geant4Vertex* v : (*i).second ){
 
        int num_part = 0;
        if ( !insideWorldVolume(context->world(), G4ThreeVector(v->x, v->y, v->z)) )   {
          if ( primaryConfig.m_skipParticlesOutsideWorldVolume )   {
            caller->warning("+++ Dropping primary vertex at (%+.2e,%+.2e,%+.2e) [mm]: outside the world volume",
                            v->x/CLHEP::mm, v->y/CLHEP::mm, v->z/CLHEP::mm);
            continue;
          }
          caller->except("Primary vertex at (%+.2e,%+.2e,%+.2e) [mm] is outside the world volume. "
                         "Geant4 cannot track particles starting there. Enlarge the world volume, or set "
                         "SkipParticlesOutsideWorldVolume=True to drop such vertices instead.",
                         v->x/CLHEP::mm, v->y/CLHEP::mm, v->z/CLHEP::mm);
        }
        G4PrimaryVertex* v4 = new G4PrimaryVertex(v->x,v->y,v->z,v->time);
        event->AddPrimaryVertex(v4);
        caller->print("+++++ G4PrimaryVertex at (%+.2e,%+.2e,%+.2e) [mm] %+.2e [ns]",
                      v->x/CLHEP::mm,v->y/CLHEP::mm,v->z/CLHEP::mm,v->time/CLHEP::ns);
        for(Geant4Vertex::Particles::const_iterator ip=v->out.begin(); ip!=v->out.end(); ++ip)  {
          Geant4ParticleHandle p = pm[*ip];
          if ( p->daughters.size() > 0 )  {
            PropertyMask mask(p->reason);
            mask.set(G4PARTICLE_HAS_SECONDARIES);
          }
          if ( p->parents.size() == 0 )  {
            Primaries relevant = getRelevant(visited,prim,pm,primaryConfig,p);
            for(Primaries::const_iterator j=relevant.begin(); j!= relevant.end(); ++j)  {
              Geant4ParticleHandle r = (*j).first;
              G4PrimaryParticle* p4 = (*j).second;
              PropertyMask reason(r->reason);
              char text[64];
              
              reason.set(G4PARTICLE_PRIMARY);
              v4->SetPrimary(p4);
              ::snprintf(text,sizeof(text),"-> G4Primary[%3d]",num_part);
              r.dumpWithMomentum(caller->outputLevel()-1,caller->name(),text);
              ++num_part;
            }
          }
        }
        if(caller->outputLevel() <= VERBOSE){
          v4->Print();
        }
      }
    }
    for( const auto& vtx : prim )   {
      Geant4ParticleHandle p = pm[vtx.first];
      primaries->insert(vtx.second, p);
    }
  }
  return 1;
}