File indexing completed on 2025-01-18 09:14:27
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
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
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
0039 #include <set>
0040 #include <algorithm>
0041
0042 using namespace dd4hep::sim;
0043 using PropertyMask = dd4hep::detail::ReferenceBitMask<int>;
0044
0045
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
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);
0064 m_needsControl = true;
0065 }
0066
0067
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);
0079 m_needsControl = true;
0080 }
0081
0082
0083 Geant4ParticleHandler::~Geant4ParticleHandler() {
0084 clear();
0085 detail::releasePtr(m_userHandler);
0086 InstanceCount::decrement(this);
0087 }
0088
0089
0090 Geant4ParticleHandler& Geant4ParticleHandler::operator=(const Geant4ParticleHandler&) {
0091 return *this;
0092 }
0093
0094
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
0112 void Geant4ParticleHandler::clear() {
0113 detail::releaseObjects(m_particleMap);
0114 m_particleMap.clear();
0115
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
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
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
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
0150 void Geant4ParticleHandler::mark(const G4Track* track) {
0151 PropertyMask mask(m_currTrack.reason);
0152 mask.set(G4PARTICLE_CREATED_HIT);
0153
0154
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
0164 mask.set(G4PARTICLE_CREATED_TRACKER_HIT);
0165
0166
0167 }
0168
0169
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
0176 if ( m_userHandler ) {
0177 m_userHandler->generate(event, this);
0178 }
0179 }
0180
0181
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
0188
0189
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
0197 if ( m_userHandler ) {
0198 m_userHandler->step(step_value, mgr, m_currTrack);
0199 }
0200 }
0201
0202
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
0213 if (m_haveSuspended) {
0214
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
0221 existingParticle = m_suspendedPM.find(h.id());
0222 if(existingParticle != m_suspendedPM.end()) {
0223 m_currTrack.get_data(*(existingParticle->second));
0224
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
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
0305 if ( m_userHandler ) {
0306 m_userHandler->begin(track, m_currTrack);
0307 }
0308 }
0309
0310
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
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
0329 PropertyMask simStatus(m_currTrack.status);
0330
0331
0332 const G4Step* theLastStep = track->GetStep();
0333 G4StepPoint* theLastPostStepPoint = NULL;
0334 if(theLastStep) theLastPostStepPoint = theLastStep->GetPostStepPoint();
0335 if( theLastPostStepPoint &&
0336 ( theLastPostStepPoint->GetStepStatus() == fWorldBoundary
0337
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
0348 if ( m_userHandler ) {
0349 m_userHandler->end(track, m_currTrack);
0350 }
0351
0352
0353
0354
0355
0356
0357
0358
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
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
0380
0381
0382
0383 int pid = m_currTrack.g4Parent;
0384 m_equivalentTracks[g4_id] = pid;
0385
0386
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;
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
0402 if(m_particleMap.find(g4_id) != m_particleMap.end()) return;
0403
0404 auto iPart = m_suspendedPM.emplace(g4_id, new Particle());
0405 (iPart.first->second)->get_data(m_currTrack);
0406 return;
0407 }
0408
0409 }
0410
0411
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
0420 if ( m_userHandler ) {
0421 m_userHandler->begin(event);
0422 }
0423 }
0424
0425
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
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
0445 rebaseSimulatedTracks(0);
0446 if ( level <= VERBOSE ) dumpMap( "Rebased ");
0447
0448 checkConsistency();
0449
0450 if ( m_userHandler ) {
0451 m_userHandler->end(event);
0452 }
0453 setVertexEndpointBit();
0454
0455
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
0463 void Geant4ParticleHandler::rebaseSimulatedTracks(int ) {
0464
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
0474
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
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
0489 orgParticles[p->id] = count;
0490 finalParticles[count] = p;
0491 p->id = count;
0492 ++count;
0493 }
0494 }
0495
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;
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;
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
0525
0526
0527
0528
0529
0530
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;
0538 if ( (ipar=finalParticles.find(equiv_id)) != finalParticles.end() ) {
0539 Particle* q = (*ipar).second;
0540 bool prim = (p->reason&G4PARTICLE_PRIMARY) == G4PARTICLE_PRIMARY;
0541
0542
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
0572
0573
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
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
0600 if ( mask.isNull() || (secondaries && low_energy && !hits_produced) ) {
0601 return true;
0602 }
0603
0604 else if ( !hits_produced && low_energy ) {
0605 return true;
0606 }
0607
0608 else if ( !tracker_track && calo_track && low_energy ) {
0609 return true;
0610 }
0611 else {
0612
0613 }
0614 return false;
0615 }
0616
0617
0618
0619 int Geant4ParticleHandler::recombineParents() {
0620 std::set<int> remove;
0621
0622
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
0627
0628
0629
0630
0631
0632
0633
0634 bool remove_me = m_userHandler ? m_userHandler->keepParticle(*p) : defaultKeepParticle(*p);
0635
0636
0637 if ( mask.isNull() || mask.isSet(G4PARTICLE_FORCE_KILL) ) {
0638 remove_me = true;
0639 }
0640 else if ( mask.isSet(G4PARTICLE_KEEP_USER) ) {
0641
0642 mask.set(G4PARTICLE_KEEP_USER);
0643 continue;
0644 }
0645 else if ( mask.isSet(G4PARTICLE_PRIMARY) ) {
0646
0647 continue;
0648 }
0649 else if ( mask.isSet(G4PARTICLE_KEEP_ALWAYS) ) {
0650 continue;
0651 }
0652 else if ( mask.isSet(G4PARTICLE_KEEP_PARENT) ) {
0653
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
0665
0666 }
0667
0668
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
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
0695 void Geant4ParticleHandler::checkConsistency() const {
0696 int num_errors = 0;
0697
0698
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
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
0714
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 }