Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:14:27

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/Primitives.h>
0016 #include <DD4hep/InstanceCount.h>
0017 #include <DDG4/Geant4StepHandler.h>
0018 #include <DDG4/Geant4TrackHandler.h>
0019 #include <DDG4/Geant4EventAction.h>
0020 #include <DDG4/Geant4SensDetAction.h>
0021 #include <DDG4/Geant4TrackingAction.h>
0022 #include <DDG4/Geant4SteppingAction.h>
0023 #include <DDG4/Geant4ParticleHandler.h>
0024 #include <DDG4/Geant4ParticleInformation.h>
0025 #include <DDG4/Geant4UserParticleHandler.h>
0026 
0027 // Geant4 include files
0028 #include <G4Step.hh>
0029 #include <G4Track.hh>
0030 #include <G4Event.hh>
0031 #include <G4TrackStatus.hh>
0032 #include <G4PrimaryVertex.hh>
0033 #include <G4PrimaryParticle.hh>
0034 #include <G4TrackingManager.hh>
0035 #include <G4ParticleDefinition.hh>
0036 #include <CLHEP/Units/SystemOfUnits.h>
0037 
0038 // C/C++ include files
0039 #include <set>
0040 #include <algorithm>
0041 
0042 using namespace dd4hep::sim;
0043 using PropertyMask = dd4hep::detail::ReferenceBitMask<int>;
0044 
0045 /// Standard constructor
0046 Geant4ParticleHandler::Geant4ParticleHandler(Geant4Context* ctxt, const std::string& nam)
0047   : Geant4GeneratorAction(ctxt,nam), Geant4MonteCarloTruth(),
0048     m_ownsParticles(false), m_userHandler(0), m_primaryMap(0)
0049 {
0050   InstanceCount::increment(this);
0051   //generatorAction().adopt(this);
0052   eventAction().callAtBegin(this,    &Geant4ParticleHandler::beginEvent);
0053   eventAction().callAtEnd(this,      &Geant4ParticleHandler::endEvent);
0054   trackingAction().callAtFinal(this, &Geant4ParticleHandler::end,CallbackSequence::FRONT);
0055   trackingAction().callUpFront(this, &Geant4ParticleHandler::begin,CallbackSequence::FRONT);
0056   steppingAction().call(this,        &Geant4ParticleHandler::step);
0057   m_globalParticleID = 0;
0058   declareProperty("PrintEndTracking",      m_printEndTracking = false);
0059   declareProperty("PrintStartTracking",    m_printStartTracking = false);
0060   declareProperty("KeepAllParticles",      m_keepAll = false);
0061   declareProperty("SaveProcesses",         m_processNames);
0062   declareProperty("MinimalKineticEnergy",  m_kinEnergyCut = 100e0*CLHEP::MeV);
0063   declareProperty("MinDistToParentVertex", m_minDistToParentVertex = 2.2e-14*CLHEP::mm);//default tolerance for g4ThreeVector isNear
0064   m_needsControl = true;
0065 }
0066 
0067 /// No default constructor
0068 Geant4ParticleHandler::Geant4ParticleHandler()
0069   : Geant4GeneratorAction(0,""), Geant4MonteCarloTruth(),
0070     m_ownsParticles(false), m_userHandler(0), m_primaryMap(0)
0071 {
0072   m_globalParticleID = 0;
0073   declareProperty("PrintEndTracking",      m_printEndTracking = false);
0074   declareProperty("PrintStartTracking",    m_printStartTracking = false);
0075   declareProperty("KeepAllParticles",      m_keepAll = false);
0076   declareProperty("SaveProcesses",         m_processNames);
0077   declareProperty("MinimalKineticEnergy",  m_kinEnergyCut = 100e0*CLHEP::MeV);
0078   declareProperty("MinDistToParentVertex", m_minDistToParentVertex = 2.2e-14*CLHEP::mm);//default tolerance for g4ThreeVector isNear
0079   m_needsControl = true;
0080 }
0081 
0082 /// Default destructor
0083 Geant4ParticleHandler::~Geant4ParticleHandler()  {
0084   clear();
0085   detail::releasePtr(m_userHandler);
0086   InstanceCount::decrement(this);
0087 }
0088 
0089 /// No assignment operator
0090 Geant4ParticleHandler& Geant4ParticleHandler::operator=(const Geant4ParticleHandler&) {
0091   return *this;
0092 }
0093 
0094 /// Adopt the user particle handler
0095 bool Geant4ParticleHandler::adopt(Geant4Action* action)    {
0096   if ( action )   {
0097     if ( !m_userHandler )  {
0098       if ( Geant4UserParticleHandler* h = dynamic_cast<Geant4UserParticleHandler*>(action) )  {
0099         m_userHandler = h;
0100         m_userHandler->addRef();
0101         return true;
0102       }
0103       except("Cannot add an invalid user particle handler object [Invalid-object-type].");
0104     }
0105     except("Cannot add an user particle handler object [Object-exists].");
0106   }
0107   except("Cannot add an invalid user particle handler object [NULL-object].");
0108   return false;
0109 }
0110 
0111 /// Clear particle maps
0112 void Geant4ParticleHandler::clear()  {
0113   detail::releaseObjects(m_particleMap);
0114   m_particleMap.clear();
0115   // m_suspendedPM should already be empty and cleared...
0116   assert(m_suspendedPM.empty() && "There was something wrong with the particle record treatment, please open a bug report!");
0117   m_equivalentTracks.clear();
0118 }
0119 
0120 /// Mark a Geant4 track to be kept for later MC truth analysis
0121 void Geant4ParticleHandler::mark(const G4Track* track, int reason)   {
0122   if ( track )   {
0123     if ( reason != 0 )  {
0124       PropertyMask(m_currTrack.reason).set(reason);
0125       return;
0126     }
0127   }
0128   except("Cannot mark the G4Track if the pointer is invalid!");
0129 }
0130 
0131 /// Store a track produced in a step to be kept for later MC truth analysis
0132 void Geant4ParticleHandler::mark(const G4Step* step_value, int reason)   {
0133   if ( step_value )  {
0134     mark(step_value->GetTrack(),reason);
0135     return;
0136   }
0137   except("Cannot mark the G4Track if the step-pointer is invalid!");
0138 }
0139 
0140 /// Mark a Geant4 track of the step to be kept for later MC truth analysis
0141 void Geant4ParticleHandler::mark(const G4Step* step_value)   {
0142   if ( step_value )  {
0143     mark(step_value->GetTrack());
0144     return;
0145   }
0146   except("Cannot mark the G4Track if the step-pointer is invalid!");
0147 }
0148 
0149 /// Mark a Geant4 track of the step to be kept for later MC truth analysis
0150 void Geant4ParticleHandler::mark(const G4Track* track)   {
0151   PropertyMask mask(m_currTrack.reason);
0152   mask.set(G4PARTICLE_CREATED_HIT);
0153   /// Check if the track origines from the calorimeter.
0154   // If yes, flag it, because it is a candidate for removal.
0155   G4LogicalVolume*       vol = track->GetVolume()->GetLogicalVolume();
0156   G4VSensitiveDetector*   g4 = vol->GetSensitiveDetector();
0157   Geant4ActionSD* sd = dynamic_cast<Geant4ActionSD*>(g4);
0158   std::string typ = sd ? sd->sensitiveType() : std::string();
0159   if ( typ == "calorimeter" )
0160     mask.set(G4PARTICLE_CREATED_CALORIMETER_HIT);
0161   else if ( typ == "tracker" )
0162     mask.set(G4PARTICLE_CREATED_TRACKER_HIT);
0163   else // Assume by default "tracker"
0164     mask.set(G4PARTICLE_CREATED_TRACKER_HIT);
0165 
0166   //Geant4ParticleHandle(&m_currTrack).dump4(outputLevel(),vol->GetName(),"hit created by particle");
0167 }
0168 
0169 /// Event generation action callback
0170 void Geant4ParticleHandler::operator()(G4Event* event)  {
0171   typedef Geant4MonteCarloTruth _MC;
0172   debug("+++ Event:%d Add EVENT extension of type Geant4ParticleHandler.....",event->GetEventID());
0173   context()->event().addExtension((_MC*)this, false);
0174   clear();
0175   /// Call the user particle handler
0176   if ( m_userHandler )  {
0177     m_userHandler->generate(event, this);
0178   }
0179 }
0180 
0181 /// User stepping callback
0182 void Geant4ParticleHandler::step(const G4Step* step_value, G4SteppingManager* mgr)   {
0183   typedef std::vector<const G4Track*> _Sec;
0184   ++m_currTrack.steps;
0185   if ( (m_currTrack.reason&G4PARTICLE_ABOVE_ENERGY_THRESHOLD) )  {
0186     //
0187     // Tracks below the energy threshold are NOT stored.
0188     // If these tracks produce hits or are selected due to another signature,
0189     // this criterium will anyhow take precedence.
0190     //
0191     const _Sec* sec=step_value->GetSecondaryInCurrentStep();
0192     if ( not sec->empty() )  {
0193       PropertyMask(m_currTrack.reason).set(G4PARTICLE_HAS_SECONDARIES);
0194     }
0195   }
0196   /// Update of the particle using the user handler
0197   if ( m_userHandler )  {
0198     m_userHandler->step(step_value, mgr, m_currTrack);
0199   }
0200 }
0201 
0202 /// Pre-track action callback
0203 void Geant4ParticleHandler::begin(const G4Track* track)   {
0204   Geant4TrackHandler   h(track);
0205   double               kine = h.kineticEnergy();
0206   G4ThreeVector        mom = h.momentum();
0207   const G4ThreeVector& v = h.vertex();
0208   int                  reason = (kine > m_kinEnergyCut) ? G4PARTICLE_ABOVE_ENERGY_THRESHOLD : 0;
0209   const G4PrimaryParticle* prim = h.primary();
0210   Particle* prim_part = 0;
0211 
0212   // if particles are not tracked to the end, we pick up where we stopped previously
0213   if (m_haveSuspended) {
0214     //primary particles are already in the particle map, we don't have to store them in another map
0215     auto existingParticle = m_particleMap.find(h.id());
0216     if(existingParticle != m_particleMap.end()) {
0217       m_currTrack.get_data(*(existingParticle->second));
0218       return;
0219     }
0220     //other particles might not be in the particleMap yet, so we take them from here
0221     existingParticle = m_suspendedPM.find(h.id());
0222     if(existingParticle != m_suspendedPM.end()) {
0223       m_currTrack.get_data(*(existingParticle->second));
0224       // make sure we delete a suspended particle in the map, fill it back later...
0225       delete (*existingParticle).second;
0226       m_suspendedPM.erase(existingParticle);
0227       return;
0228     }
0229   }
0230 
0231   if ( prim )   {
0232     prim_part = m_primaryMap->get(prim);
0233     if ( !prim_part )  {
0234       except("+++ Tracking preaction: Primary particle without generator particle!");
0235     }
0236     reason |= (G4PARTICLE_PRIMARY|G4PARTICLE_ABOVE_ENERGY_THRESHOLD);
0237     m_particleMap[h.id()] = prim_part->addRef();
0238   }
0239 
0240   if ( prim_part )   {
0241     m_currTrack.id           = prim_part->id;
0242     m_currTrack.reason       = prim_part->reason|reason;
0243     m_currTrack.mask         = prim_part->mask;
0244     m_currTrack.status       = prim_part->status;
0245     m_currTrack.genStatus    = prim_part->genStatus;
0246     m_currTrack.spin[0]      = prim_part->spin[0];
0247     m_currTrack.spin[1]      = prim_part->spin[1];
0248     m_currTrack.spin[2]      = prim_part->spin[2];
0249     m_currTrack.colorFlow[0] = prim_part->colorFlow[0];
0250     m_currTrack.colorFlow[1] = prim_part->colorFlow[1];
0251     m_currTrack.parents      = prim_part->parents;
0252     m_currTrack.daughters    = prim_part->daughters;
0253     m_currTrack.pdgID        = prim_part->pdgID;
0254     m_currTrack.mass         = prim_part->mass;
0255     m_currTrack.charge       = int(3.0 * h.charge());
0256   }
0257   else  {
0258     m_currTrack.id           = m_globalParticleID;
0259     m_currTrack.reason       = reason;
0260     m_currTrack.mask         = 0;
0261     m_currTrack.status       = G4PARTICLE_SIM_CREATED;
0262     m_currTrack.genStatus    = 0;
0263     m_currTrack.spin[0]      = 0;
0264     m_currTrack.spin[1]      = 0;
0265     m_currTrack.spin[2]      = 0;
0266     m_currTrack.colorFlow[0] = 0;
0267     m_currTrack.colorFlow[1] = 0;
0268     m_currTrack.parents.clear();
0269     m_currTrack.daughters.clear();
0270     m_currTrack.pdgID        = h.pdgID();
0271     m_currTrack.mass         = h.mass();
0272     m_currTrack.charge       = int(3.0 * h.charge());
0273     ++m_globalParticleID;
0274   }
0275   m_currTrack.steps       = 0;
0276   m_currTrack.secondaries = 0;
0277   m_currTrack.g4Parent    = h.parent();
0278   m_currTrack.originalG4ID= h.id();
0279   m_currTrack.process     = h.creatorProcess();
0280   m_currTrack.time        = h.globalTime();
0281   m_currTrack.vsx         = v.x();
0282   m_currTrack.vsy         = v.y();
0283   m_currTrack.vsz         = v.z();
0284   m_currTrack.vex         = 0.0;
0285   m_currTrack.vey         = 0.0;
0286   m_currTrack.vez         = 0.0;
0287   m_currTrack.psx         = mom.x();
0288   m_currTrack.psy         = mom.y();
0289   m_currTrack.psz         = mom.z();
0290   m_currTrack.pex         = 0.0;
0291   m_currTrack.pey         = 0.0;
0292   m_currTrack.pez         = 0.0;
0293   // If the creator process of the track is in the list of process products to be kept, set the proper flag
0294   if ( m_currTrack.process )  {
0295     Processes::iterator i=find(m_processNames.begin(),m_processNames.end(),m_currTrack.process->GetProcessName());
0296     if ( i != m_processNames.end() )  {
0297       PropertyMask(m_currTrack.reason).set(G4PARTICLE_KEEP_PROCESS);
0298     }
0299   }
0300   if ( m_keepAll )  {
0301     PropertyMask(m_currTrack.reason).set(G4PARTICLE_KEEP_ALWAYS);
0302   }
0303 
0304   /// Initial update of the particle using the user handler
0305   if ( m_userHandler )  {
0306     m_userHandler->begin(track, m_currTrack);
0307   }
0308 }
0309 
0310 /// Post-track action callback
0311 void Geant4ParticleHandler::end(const G4Track* track)   {
0312   Geant4TrackHandler h(track);
0313   Geant4ParticleHandle ph(&m_currTrack);
0314   const int g4_id = h.id();
0315 
0316   int track_reason = m_currTrack.reason;
0317   PropertyMask mask(m_currTrack.reason);
0318   // Update vertex end point and final momentum
0319   G4ThreeVector mom = track->GetMomentum();
0320   const G4ThreeVector& pos = track->GetPosition();
0321   ph->pex = mom.x();
0322   ph->pey = mom.y();
0323   ph->pez = mom.z();
0324   ph->vex = pos.x();
0325   ph->vey = pos.y();
0326   ph->vez = pos.z();
0327 
0328   // Set the simulator status bits
0329   PropertyMask simStatus(m_currTrack.status);
0330 
0331   // check if the last step ended on the worldVolume boundary
0332   const G4Step* theLastStep = track->GetStep();
0333   G4StepPoint* theLastPostStepPoint = NULL;
0334   if(theLastStep) theLastPostStepPoint = theLastStep->GetPostStepPoint();
0335   if( theLastPostStepPoint &&
0336       ( theLastPostStepPoint->GetStepStatus() == fWorldBoundary //particle left world volume
0337     //|| theLastPostStepPoint->GetStepStatus() == fGeomBoundary
0338       )
0339     ) {
0340     simStatus.set(G4PARTICLE_SIM_LEFT_DETECTOR);
0341   }
0342 
0343   if(track->GetKineticEnergy() <= 0.) {
0344     simStatus.set(G4PARTICLE_SIM_STOPPED);
0345   }
0346 
0347   /// Final update of the particle using the user handler
0348   if ( m_userHandler )  {
0349     m_userHandler->end(track, m_currTrack);
0350   }
0351  
0352   // These are candidate tracks with a probability to be stored due to their properties:
0353   // - primary particle
0354   // - hits created
0355   // - secondaries
0356   // - above energy threshold
0357   // - to be kept due to creator process
0358   // - to be kept due to user information of type 'Geant4ParticleInformation' stored in the G4Track
0359   //
0360   Geant4ParticleInformation* track_info =
0361     dynamic_cast<Geant4ParticleInformation*>(track->GetUserInformation());
0362   if ( !mask.isNull() || track_info )   {
0363     m_equivalentTracks[g4_id] = g4_id;
0364     ParticleMap::iterator ip = m_particleMap.find(g4_id);
0365     if ( mask.isSet(G4PARTICLE_PRIMARY) )   {
0366       ph.dump2(outputLevel()-1,name(),"Add Primary",h.id(),ip!=m_particleMap.end());
0367     }
0368     // Create a new MC particle from the current track information saved in the pre-tracking action
0369     Particle* part = 0;
0370     if ( ip==m_particleMap.end() ) part = m_particleMap[g4_id] = new Particle();
0371     else part = (*ip).second;
0372     if ( track_info )  {
0373       mask.set(G4PARTICLE_KEEP_USER);
0374       part->extension.reset(track_info->release());
0375     }
0376     part->get_data(m_currTrack);
0377   }
0378   else   {
0379     // These are tracks without any special properties.
0380     //
0381     // We will not store them on the record, but have to memorise the
0382     // track identifier in order to restore the history for the created hits.
0383     int pid = m_currTrack.g4Parent;
0384     m_equivalentTracks[g4_id] = pid;
0385     // Need to find the last stored particle and OR this particle's mask
0386     // with the mask of the last stored particle
0387     auto iend = m_equivalentTracks.end(), iequiv=m_equivalentTracks.end();
0388     ParticleMap::iterator ip;
0389     for(ip=m_particleMap.find(pid); ip == m_particleMap.end(); ip=m_particleMap.find(pid))  {
0390       if (iequiv=m_equivalentTracks.find(pid); iequiv == iend) break;  // ERROR
0391       pid = (*iequiv).second;
0392     }
0393     if ( ip != m_particleMap.end() )
0394       (*ip).second->reason |= track_reason;
0395     else
0396       ph.dumpWithVertex(outputLevel()+3,name(),"FATAL: No real particle parent present");
0397   }
0398 
0399   if(track->GetTrackStatus() == fSuspend) {
0400     m_haveSuspended = true;
0401     //track is already in particle map, we pick it up from there in begin again
0402     if(m_particleMap.find(g4_id) != m_particleMap.end()) return;
0403     //track is not already stored, keep it in special map
0404     auto iPart = m_suspendedPM.emplace(g4_id, new Particle());
0405     (iPart.first->second)->get_data(m_currTrack);
0406     return; // we trust that we eventually return to this function with another status and go on then
0407   }
0408 
0409 }
0410 
0411 /// Pre-event action callback
0412 void Geant4ParticleHandler::beginEvent(const G4Event* event)  {
0413   Geant4PrimaryInteraction* interaction = context()->event().extension<Geant4PrimaryInteraction>();
0414   info("+++ Event %d Begin event action. Access event related information.",event->GetEventID());
0415   m_primaryMap = context()->event().extension<Geant4PrimaryMap>();
0416   m_globalParticleID = interaction->nextPID();
0417   m_particleMap.clear();
0418   m_equivalentTracks.clear();
0419   /// Call the user particle handler
0420   if ( m_userHandler )  {
0421     m_userHandler->begin(event);
0422   }
0423 }
0424 
0425 /// Debugging: Dump Geant4 particle map
0426 void Geant4ParticleHandler::dumpMap(const char* tag)  const  {
0427   const std::string& n = name();
0428   Geant4ParticleHandle::header4(INFO,n,tag);
0429   for(ParticleMap::const_iterator iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i)  {
0430     Geant4ParticleHandle((*i).second).dump4(INFO,n,tag);
0431   }
0432 }
0433 
0434 /// Post-event action callback
0435 void Geant4ParticleHandler::endEvent(const G4Event* event)  {
0436   int count = 0;
0437   int level = outputLevel();
0438   do {
0439     if ( level <= VERBOSE ) dumpMap("Particle  ");
0440     debug("+++ Iteration:%d Tracks:%d Equivalents:%d",++count,m_particleMap.size(),m_equivalentTracks.size());
0441   } while( recombineParents() > 0 );
0442 
0443   if ( level <= VERBOSE ) dumpMap(  "Recombined");
0444   // Rebase the simulated tracks, so that they fit to the generator particles
0445   rebaseSimulatedTracks(0);
0446   if ( level <= VERBOSE ) dumpMap(  "Rebased   ");
0447   // Consistency check....
0448   checkConsistency();
0449   /// Call the user particle handler
0450   if ( m_userHandler )  {
0451     m_userHandler->end(event);
0452   }
0453   setVertexEndpointBit();
0454 
0455   // Now export the data to the final record.
0456   Geant4ParticleMap* part_map = context()->event().extension<Geant4ParticleMap>();
0457   part_map->adopt(m_particleMap, m_equivalentTracks);
0458   m_primaryMap = 0;
0459   clear();
0460 }
0461 
0462 /// Rebase the simulated tracks, so that they fit to the generator particles
0463 void Geant4ParticleHandler::rebaseSimulatedTracks(int )   {
0464   /// No we have to update the map of equivalent tracks and assign the 'equivalentTrack' entry
0465   TrackEquivalents equivalents, orgParticles;
0466   ParticleMap      finalParticles;
0467   ParticleMap::const_iterator ipar, iend, i;
0468   int count;
0469 
0470   Geant4PrimaryInteraction* interaction = context()->event().extension<Geant4PrimaryInteraction>();
0471   ParticleMap& pm = interaction->particles;
0472 
0473   // (1.0) Copy the pre-defined particle mapping for the simulated tracks
0474   //       It is assumed the mapping is ZERO based without holes.
0475   for(count = 0, iend=pm.end(), i=pm.begin(); i!=iend; ++i)  {
0476     Particle* p = (*i).second;
0477     orgParticles[p->id] = p->id;
0478     finalParticles[p->id] = p;
0479     if ( p->id > count ) count = p->id;
0480     if ( (p->reason&G4PARTICLE_PRIMARY) != G4PARTICLE_PRIMARY )  {
0481       p->addRef();
0482     }
0483   }
0484   // (1.1) Define the new particle mapping for the simulated tracks
0485   for(++count, iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i)  {
0486     Particle* p = (*i).second;
0487     if ( (p->reason&G4PARTICLE_PRIMARY) != G4PARTICLE_PRIMARY )  {
0488       //if ( orgParticles.find(p->id) == orgParticles.end() )  {
0489       orgParticles[p->id] = count;
0490       finalParticles[count] = p;
0491       p->id = count;
0492       ++count;
0493     }
0494   }
0495   // (2) Re-evaluate the corresponding geant4 track equivalents using the new mapping
0496   for(TrackEquivalents::iterator ie=m_equivalentTracks.begin(),ie_end=m_equivalentTracks.end(); ie!=ie_end; ++ie)  {
0497     int g4_equiv = (*ie).first;
0498     while( (ipar=m_particleMap.find(g4_equiv)) == m_particleMap.end() )  {
0499       TrackEquivalents::const_iterator iequiv = m_equivalentTracks.find(g4_equiv);
0500       if ( iequiv == ie_end )  {
0501         break;  // ERROR !! Will be handled by printout below because ipar==end()
0502       }
0503       g4_equiv = (*iequiv).second;
0504     }
0505     TrackEquivalents::mapped_type equiv = (*ie).second;
0506     if ( ipar != m_particleMap.end() )   {
0507       Geant4ParticleHandle p = (*ipar).second;
0508       equivalents[(*ie).first] = p->id;  // requires (1) to be filled properly!
0509       const G4ParticleDefinition* def = p.definition();
0510       int pdg = int(std::abs(def->GetPDGEncoding())+0.1);
0511       if ( pdg != 0 && pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 )  {
0512         error("+++ ERROR: Geant4 particle for track:%d last known is:%d -- is gluon or quark!",equiv,g4_equiv);
0513       }
0514       pdg = int(std::abs(p->pdgID)+0.1);
0515       if ( pdg != 0 && pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 )  {
0516         error("+++ ERROR(2): Geant4 particle for track:%d last known is:%d -- is gluon or quark!",equiv,g4_equiv);
0517       }
0518     }
0519     else   {
0520       error("+++ No Equivalent particle for track:%d last known is:%d",equiv,g4_equiv);
0521     }
0522   }
0523 
0524   // (3) Compute the particle's parents and daughters.
0525   //     Replace the original Geant4 track with the
0526   //     equivalent particle still present in the record.
0527   // Note:
0528   //     We rely here on the ordering of the particles accoding to their
0529   //     Processing by Geant4 to establish mother daughter relationships.
0530   //     == > use finalParticles map and NOT m_particleMap.
0531   int equiv_id = -1;
0532   for( auto& part : finalParticles )  {
0533     auto& p = part.second;
0534     if ( p->g4Parent > 0 )  {
0535       TrackEquivalents::iterator iequ = equivalents.find(p->g4Parent);
0536       if ( iequ != equivalents.end() )  {
0537         equiv_id = (*iequ).second;//equivalents[p->g4Parent];
0538         if ( (ipar=finalParticles.find(equiv_id)) != finalParticles.end() )  {
0539           Particle* q = (*ipar).second;
0540           bool      prim = (p->reason&G4PARTICLE_PRIMARY) == G4PARTICLE_PRIMARY;
0541           // We assume that the mother daughter relationship
0542           // is filled by the event readers!
0543           if ( !prim )  {
0544             p->parents.insert(q->id);
0545           }
0546           if ( !p->parents.empty() )  {
0547             int parent_id = (*p->parents.begin());
0548             if ( parent_id == q->id )
0549               q->daughters.insert(p->id);
0550             else if ( !prim )
0551               error("+++ Inconsistency in equivalent record! Parent: %d Daughter:%d",q->id, p->id);
0552           }
0553           else   {
0554             error("+++ Inconsistency in parent relashionship: %d NO parent!", p->id);
0555           }
0556           continue;
0557         }
0558       }
0559       error("+++ Inconsistency in particle record: Geant4 parent %d "
0560             "of particle %d not in record of final particles!",
0561             p->g4Parent,p->id);
0562     }
0563   }
0564 #if 0
0565   for(iend=finalParticles.end(), i=finalParticles.begin(); i!=iend; ++i)  {
0566     Particle* p = (*i).second;
0567     if ( p->g4Parent > 0 )  {
0568       int parent_id = (*p->parents.begin());
0569       if ( (ipar=finalParticles.find(parent_id)) != finalParticles.end() )  {
0570         Particle* q = (*ipar).second;
0571         // Generator particles have a proper history.
0572         // We only deal with particles, which are not of MC origin.
0573         //p->parents.insert(q->id);
0574         if ( parent_id == q->id )
0575           q->daughters.insert(p->id);
0576         else
0577           error("+++ Inconsistency in equivalent record! Parent: %d Daughter:%d",q->id, p->id);
0578         continue;
0579       }
0580       error("+++ Inconsistency in particle record: Geant4 parent %d "
0581             "of particle %d not in record of final particles!",
0582             p->g4Parent,p->id);
0583     }
0584   }
0585 #endif
0586   m_equivalentTracks = std::move(equivalents);
0587   m_particleMap = std::move(finalParticles);
0588 }
0589 
0590 /// Default callback to be answered if the particle should be kept if NO user handler is installed
0591 bool Geant4ParticleHandler::defaultKeepParticle(Particle& particle)   {
0592   PropertyMask mask(particle.reason);
0593   bool secondaries    =  mask.isSet(G4PARTICLE_HAS_SECONDARIES);
0594   bool tracker_track  =  mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT);
0595   bool calo_track     =  mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT);
0596   bool hits_produced  =  mask.isSet(G4PARTICLE_CREATED_HIT);
0597   bool low_energy     = !mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD);
0598 
0599   /// Remove this track if it has not created a hit and the energy is below threshold
0600   if ( mask.isNull() || (secondaries && low_energy && !hits_produced) )  {
0601     return true;
0602   }
0603   /// Remove this track if the energy is below threshold. Reassign hits to parent.
0604   else if ( !hits_produced && low_energy )  {
0605     return true;
0606   }
0607   /// Remove this track if the origine is in the calorimeter. Reassign hits to parent.
0608   else if ( !tracker_track && calo_track && low_energy )  {
0609     return true;
0610   }
0611   else  {
0612     //printout(INFO,name(),"+++ Track: %d should be kept for no obvious reason....",id);
0613   }
0614   return false;
0615 }
0616 
0617 /// Clean the monte carlo record. Remove all unwanted stuff.
0618 /// This is the core of the object executed at the end of each event action.
0619 int Geant4ParticleHandler::recombineParents()  {
0620   std::set<int> remove;
0621 
0622   /// Need to start from BACK, to clean first the latest produced stuff.
0623   for(ParticleMap::reverse_iterator i=m_particleMap.rbegin(); i!=m_particleMap.rend(); ++i)  {
0624     Particle* p = (*i).second;
0625     PropertyMask mask(p->reason);
0626     // Allow the user to force the particle handling either by
0627     // or the reason mask with G4PARTICLE_KEEP_USER or
0628     // to set the reason mask to NULL in order to drop it.
0629     //
0630     // If the mask entry is set to G4PARTICLE_FORCE_KILL
0631     // or is set to NULL, the particle is ALWAYS removed
0632     //
0633     // Note: This may override all other decisions!
0634     bool remove_me = m_userHandler ? m_userHandler->keepParticle(*p) : defaultKeepParticle(*p);
0635 
0636     // Now look at the property mask of the particle
0637     if ( mask.isNull() || mask.isSet(G4PARTICLE_FORCE_KILL) )  {
0638       remove_me = true;
0639     }
0640     else if ( mask.isSet(G4PARTICLE_KEEP_USER) )  {
0641       /// If user decides it must be kept, it MUST be kept!
0642       mask.set(G4PARTICLE_KEEP_USER);
0643       continue;
0644     }
0645     else if ( mask.isSet(G4PARTICLE_PRIMARY) )   {
0646       /// Primary particles MUST be kept!
0647       continue;
0648     }
0649     else if ( mask.isSet(G4PARTICLE_KEEP_ALWAYS) )   {
0650       continue;
0651     }
0652     else if ( mask.isSet(G4PARTICLE_KEEP_PARENT) )  {
0653       //continue;
0654     }
0655     else if ( mask.isSet(G4PARTICLE_KEEP_PROCESS) )  {
0656       if(ParticleMap::iterator ip = m_particleMap.find(p->g4Parent); ip != m_particleMap.end() )   {
0657         Particle* parent_part = (*ip).second;
0658         PropertyMask parent_mask(parent_part->reason);
0659         if ( parent_mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD) )   {
0660           parent_mask.set(G4PARTICLE_KEEP_PARENT);
0661           continue;
0662         }
0663       }
0664       // Low energy stuff. Remove it. Reassign to parent.
0665       //remove_me = true;
0666     }
0667 
0668     /// Remove this track from the list and also do the cleanup in the parent's children list
0669     if ( remove_me )  {
0670       int g4_id = (*i).first;
0671       remove.insert(g4_id);
0672       m_equivalentTracks[g4_id] = p->g4Parent;
0673       if(ParticleMap::iterator ip = m_particleMap.find(p->g4Parent); ip != m_particleMap.end() )   {
0674         Particle* parent_part = (*ip).second;
0675         PropertyMask(parent_part->reason).set(mask.value());
0676         parent_part->steps += p->steps;
0677         parent_part->secondaries += p->secondaries;
0678         /// Update of the particle using the user handler
0679         if ( m_userHandler )  {
0680           m_userHandler->combine(*p, *parent_part);
0681         }
0682       }
0683     }
0684   }
0685   for( int r : remove )  {
0686     if( auto ir = m_particleMap.find(r); ir != m_particleMap.end() )  {
0687       (*ir).second->release();
0688       m_particleMap.erase(ir);
0689     }
0690   }
0691   return int(remove.size());
0692 }
0693 
0694 /// Check the record consistency
0695 void Geant4ParticleHandler::checkConsistency()  const   {
0696   int num_errors = 0;
0697 
0698   /// First check the consistency of the particle map itself
0699   for(const auto& part : m_particleMap )  {
0700     Geant4Particle* particle = part.second;
0701     Geant4ParticleHandle p(particle);
0702     PropertyMask mask(p->reason);
0703     PropertyMask status(p->status);
0704     std::set<int>& daughters = p->daughters;
0705     ParticleMap::const_iterator j;
0706     // For all particles, the set of daughters must be contained in the record.
0707     for( int id_dau : daughters )   {
0708       if ( j=m_particleMap.find(id_dau); j == m_particleMap.end() )   {
0709         ++num_errors;
0710         error("+++ Particle:%d Daughter %d is not in particle map!",p->id,id_dau);
0711       }
0712     }
0713     // We assume that particles from the generator have consistent parents
0714     // For all other particles except the primaries, the parent must be contained in the record.
0715     if ( !mask.isSet(G4PARTICLE_PRIMARY) && !status.anySet(G4PARTICLE_GEN_STATUS) )  {
0716       bool in_map = false, in_parent_list = false;
0717       int  parent_id = -1;
0718       if( auto eq_it=m_equivalentTracks.find(p->g4Parent); eq_it != m_equivalentTracks.end() )   {
0719         parent_id = (*eq_it).second;
0720         in_map    = (j=m_particleMap.find(parent_id)) != m_particleMap.end();
0721         in_parent_list = p->parents.find(parent_id) != p->parents.end();
0722       }
0723       if ( !in_map || !in_parent_list )  {
0724         char parent_list[1024];
0725         parent_list[0] = 0;
0726         ++num_errors;
0727         p.dumpWithMomentum(ERROR,name(),"INCONSISTENCY");
0728         for( int ip : p->parents )
0729           ::snprintf(parent_list+strlen(parent_list),sizeof(parent_list)-strlen(parent_list),"%d ",ip);
0730         error("+++ Particle:%d Parent %d (G4id:%d)  In record:%s In parent list:%s [%s]",
0731               p->id,parent_id,p->g4Parent,yes_no(in_map),yes_no(in_parent_list),parent_list);
0732       }
0733     }
0734   }
0735 
0736   if ( num_errors > 0 )  {
0737     except("+++ Consistency check failed. Found %d problems.",num_errors);
0738   }
0739 }
0740 
0741 void Geant4ParticleHandler::setVertexEndpointBit() {
0742   for( auto& part : m_particleMap )   {
0743     auto* p = part.second;
0744     if( !p->parents.empty() )   {
0745       Geant4Particle *parent(m_particleMap[ *p->parents.begin() ]);
0746       const double X( parent->vex - p->vsx );
0747       const double Y( parent->vey - p->vsy );
0748       const double Z( parent->vez - p->vsz );
0749       if( sqrt(X*X + Y*Y + Z*Z) > m_minDistToParentVertex ){
0750         PropertyMask(p->status).set(G4PARTICLE_SIM_PARENT_RADIATED);
0751       }
0752     }
0753   }
0754 }