Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-20 08:25:03 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!