//==========================================================================
// AIDA Detector description implementation
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author : M.Frank
//
//==========================================================================
// Framework include files
#include <DD4hep/Primitives.h>
#include <DD4hep/InstanceCount.h>
#include <DDG4/Geant4StepHandler.h>
#include <DDG4/Geant4TrackHandler.h>
#include <DDG4/Geant4EventAction.h>
#include <DDG4/Geant4SensDetAction.h>
#include <DDG4/Geant4TrackingAction.h>
#include <DDG4/Geant4SteppingAction.h>
#include <DDG4/Geant4ParticleHandler.h>
#include <DDG4/Geant4ParticleInformation.h>
#include <DDG4/Geant4UserParticleHandler.h>
// Geant4 include files
#include <G4Step.hh>
#include <G4Track.hh>
#include <G4Event.hh>
#include <G4TrackStatus.hh>
#include <G4PrimaryVertex.hh>
#include <G4PrimaryParticle.hh>
#include <G4TrackingManager.hh>
#include <G4ParticleDefinition.hh>
#include <CLHEP/Units/SystemOfUnits.h>
// C/C++ include files
#include <set>
#include <stdexcept>
#include <algorithm>
using namespace std;
using namespace dd4hep;
using namespace dd4hep::sim;
typedef detail::ReferenceBitMask<int> PropertyMask;
/// Standard constructor
Geant4ParticleHandler::Geant4ParticleHandler(Geant4Context* ctxt, const string& nam)
: Geant4GeneratorAction(ctxt,nam), Geant4MonteCarloTruth(),
m_ownsParticles(false), m_userHandler(0), m_primaryMap(0)
{
InstanceCount::increment(this);
//generatorAction().adopt(this);
eventAction().callAtBegin(this, &Geant4ParticleHandler::beginEvent);
eventAction().callAtEnd(this, &Geant4ParticleHandler::endEvent);
trackingAction().callAtFinal(this, &Geant4ParticleHandler::end,CallbackSequence::FRONT);
trackingAction().callUpFront(this, &Geant4ParticleHandler::begin,CallbackSequence::FRONT);
steppingAction().call(this, &Geant4ParticleHandler::step);
m_globalParticleID = 0;
declareProperty("PrintEndTracking", m_printEndTracking = false);
declareProperty("PrintStartTracking", m_printStartTracking = false);
declareProperty("KeepAllParticles", m_keepAll = false);
declareProperty("SaveProcesses", m_processNames);
declareProperty("MinimalKineticEnergy", m_kinEnergyCut = 100e0*CLHEP::MeV);
declareProperty("MinDistToParentVertex", m_minDistToParentVertex = 2.2e-14*CLHEP::mm);//default tolerance for g4ThreeVector isNear
m_needsControl = true;
}
/// No default constructor
Geant4ParticleHandler::Geant4ParticleHandler()
: Geant4GeneratorAction(0,""), Geant4MonteCarloTruth(),
m_ownsParticles(false), m_userHandler(0), m_primaryMap(0)
{
m_globalParticleID = 0;
declareProperty("PrintEndTracking", m_printEndTracking = false);
declareProperty("PrintStartTracking", m_printStartTracking = false);
declareProperty("KeepAllParticles", m_keepAll = false);
declareProperty("SaveProcesses", m_processNames);
declareProperty("MinimalKineticEnergy", m_kinEnergyCut = 100e0*CLHEP::MeV);
declareProperty("MinDistToParentVertex", m_minDistToParentVertex = 2.2e-14*CLHEP::mm);//default tolerance for g4ThreeVector isNear
m_needsControl = true;
}
/// Default destructor
Geant4ParticleHandler::~Geant4ParticleHandler() {
clear();
detail::releasePtr(m_userHandler);
InstanceCount::decrement(this);
}
/// No assignment operator
Geant4ParticleHandler& Geant4ParticleHandler::operator=(const Geant4ParticleHandler&) {
return *this;
}
/// Adopt the user particle handler
bool Geant4ParticleHandler::adopt(Geant4Action* action) {
if ( action ) {
if ( !m_userHandler ) {
if ( Geant4UserParticleHandler* h = dynamic_cast<Geant4UserParticleHandler*>(action) ) {
m_userHandler = h;
m_userHandler->addRef();
return true;
}
except("Cannot add an invalid user particle handler object [Invalid-object-type].");
}
except("Cannot add an user particle handler object [Object-exists].");
}
except("Cannot add an invalid user particle handler object [NULL-object].");
return false;
}
/// Clear particle maps
void Geant4ParticleHandler::clear() {
detail::releaseObjects(m_particleMap);
m_particleMap.clear();
// m_suspendedPM should already be empty and cleared...
assert(m_suspendedPM.empty() && "There was something wrong with the particle record treatment, please open a bug report!");
m_equivalentTracks.clear();
}
/// Mark a Geant4 track to be kept for later MC truth analysis
void Geant4ParticleHandler::mark(const G4Track* track, int reason) {
if ( track ) {
if ( reason != 0 ) {
PropertyMask(m_currTrack.reason).set(reason);
return;
}
}
except("Cannot mark the G4Track if the pointer is invalid!");
}
/// Store a track produced in a step to be kept for later MC truth analysis
void Geant4ParticleHandler::mark(const G4Step* step_value, int reason) {
if ( step_value ) {
mark(step_value->GetTrack(),reason);
return;
}
except("Cannot mark the G4Track if the step-pointer is invalid!");
}
/// Mark a Geant4 track of the step to be kept for later MC truth analysis
void Geant4ParticleHandler::mark(const G4Step* step_value) {
if ( step_value ) {
mark(step_value->GetTrack());
return;
}
except("Cannot mark the G4Track if the step-pointer is invalid!");
}
/// Mark a Geant4 track of the step to be kept for later MC truth analysis
void Geant4ParticleHandler::mark(const G4Track* track) {
PropertyMask mask(m_currTrack.reason);
mask.set(G4PARTICLE_CREATED_HIT);
/// Check if the track origines from the calorimeter.
// If yes, flag it, because it is a candidate for removal.
G4LogicalVolume* vol = track->GetVolume()->GetLogicalVolume();
G4VSensitiveDetector* g4 = vol->GetSensitiveDetector();
Geant4ActionSD* sd = dynamic_cast<Geant4ActionSD*>(g4);
string typ = sd ? sd->sensitiveType() : string();
if ( typ == "calorimeter" )
mask.set(G4PARTICLE_CREATED_CALORIMETER_HIT);
else if ( typ == "tracker" )
mask.set(G4PARTICLE_CREATED_TRACKER_HIT);
else // Assume by default "tracker"
mask.set(G4PARTICLE_CREATED_TRACKER_HIT);
//Geant4ParticleHandle(&m_currTrack).dump4(outputLevel(),vol->GetName(),"hit created by particle");
}
/// Event generation action callback
void Geant4ParticleHandler::operator()(G4Event* event) {
typedef Geant4MonteCarloTruth _MC;
debug("+++ Event:%d Add EVENT extension of type Geant4ParticleHandler.....",event->GetEventID());
context()->event().addExtension((_MC*)this, false);
clear();
/// Call the user particle handler
if ( m_userHandler ) {
m_userHandler->generate(event, this);
}
}
/// User stepping callback
void Geant4ParticleHandler::step(const G4Step* step_value, G4SteppingManager* mgr) {
typedef vector<const G4Track*> _Sec;
++m_currTrack.steps;
if ( (m_currTrack.reason&G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ) {
//
// Tracks below the energy threshold are NOT stored.
// If these tracks produce hits or are selected due to another signature,
// this criterium will anyhow take precedence.
//
const _Sec* sec=step_value->GetSecondaryInCurrentStep();
if ( not sec->empty() ) {
PropertyMask(m_currTrack.reason).set(G4PARTICLE_HAS_SECONDARIES);
}
}
/// Update of the particle using the user handler
if ( m_userHandler ) {
m_userHandler->step(step_value, mgr, m_currTrack);
}
}
/// Pre-track action callback
void Geant4ParticleHandler::begin(const G4Track* track) {
Geant4TrackHandler h(track);
double kine = h.kineticEnergy();
G4ThreeVector mom = h.momentum();
const G4ThreeVector& v = h.vertex();
int reason = (kine > m_kinEnergyCut) ? G4PARTICLE_ABOVE_ENERGY_THRESHOLD : 0;
const G4PrimaryParticle* prim = h.primary();
Particle* prim_part = 0;
// if particles are not tracked to the end, we pick up where we stopped previously
if (m_haveSuspended) {
//primary particles are already in the particle map, we don't have to store them in another map
auto existingParticle = m_particleMap.find(h.id());
if(existingParticle != m_particleMap.end()) {
m_currTrack.get_data(*(existingParticle->second));
return;
}
//other particles might not be in the particleMap yet, so we take them from here
existingParticle = m_suspendedPM.find(h.id());
if(existingParticle != m_suspendedPM.end()) {
m_currTrack.get_data(*(existingParticle->second));
// make sure we delete a suspended particle in the map, fill it back later...
delete (*existingParticle).second;
m_suspendedPM.erase(existingParticle);
return;
}
}
if ( prim ) {
prim_part = m_primaryMap->get(prim);
if ( !prim_part ) {
except("+++ Tracking preaction: Primary particle without generator particle!");
}
reason |= (G4PARTICLE_PRIMARY|G4PARTICLE_ABOVE_ENERGY_THRESHOLD);
m_particleMap[h.id()] = prim_part->addRef();
}
if ( prim_part ) {
m_currTrack.id = prim_part->id;
m_currTrack.reason = prim_part->reason|reason;
m_currTrack.mask = prim_part->mask;
m_currTrack.status = prim_part->status;
m_currTrack.genStatus = prim_part->genStatus;
m_currTrack.spin[0] = prim_part->spin[0];
m_currTrack.spin[1] = prim_part->spin[1];
m_currTrack.spin[2] = prim_part->spin[2];
m_currTrack.colorFlow[0] = prim_part->colorFlow[0];
m_currTrack.colorFlow[1] = prim_part->colorFlow[1];
m_currTrack.parents = prim_part->parents;
m_currTrack.daughters = prim_part->daughters;
m_currTrack.pdgID = prim_part->pdgID;
m_currTrack.mass = prim_part->mass;
m_currTrack.charge = int(3.0 * h.charge());
}
else {
m_currTrack.id = m_globalParticleID;
m_currTrack.reason = reason;
m_currTrack.mask = 0;
m_currTrack.status = G4PARTICLE_SIM_CREATED;
m_currTrack.genStatus = 0;
m_currTrack.spin[0] = 0;
m_currTrack.spin[1] = 0;
m_currTrack.spin[2] = 0;
m_currTrack.colorFlow[0] = 0;
m_currTrack.colorFlow[1] = 0;
m_currTrack.parents.clear();
m_currTrack.daughters.clear();
m_currTrack.pdgID = h.pdgID();
m_currTrack.mass = h.mass();
m_currTrack.charge = int(3.0 * h.charge());
++m_globalParticleID;
}
m_currTrack.steps = 0;
m_currTrack.secondaries = 0;
m_currTrack.g4Parent = h.parent();
m_currTrack.originalG4ID= h.id();
m_currTrack.process = h.creatorProcess();
m_currTrack.time = h.globalTime();
m_currTrack.vsx = v.x();
m_currTrack.vsy = v.y();
m_currTrack.vsz = v.z();
m_currTrack.vex = 0.0;
m_currTrack.vey = 0.0;
m_currTrack.vez = 0.0;
m_currTrack.psx = mom.x();
m_currTrack.psy = mom.y();
m_currTrack.psz = mom.z();
m_currTrack.pex = 0.0;
m_currTrack.pey = 0.0;
m_currTrack.pez = 0.0;
// If the creator process of the track is in the list of process products to be kept, set the proper flag
if ( m_currTrack.process ) {
Processes::iterator i=find(m_processNames.begin(),m_processNames.end(),m_currTrack.process->GetProcessName());
if ( i != m_processNames.end() ) {
PropertyMask(m_currTrack.reason).set(G4PARTICLE_KEEP_PROCESS);
}
}
if ( m_keepAll ) {
PropertyMask(m_currTrack.reason).set(G4PARTICLE_KEEP_ALWAYS);
}
/// Initial update of the particle using the user handler
if ( m_userHandler ) {
m_userHandler->begin(track, m_currTrack);
}
}
/// Post-track action callback
void Geant4ParticleHandler::end(const G4Track* track) {
Geant4TrackHandler h(track);
Geant4ParticleHandle ph(&m_currTrack);
const int g4_id = h.id();
int track_reason = m_currTrack.reason;
PropertyMask mask(m_currTrack.reason);
// Update vertex end point and final momentum
G4ThreeVector mom = track->GetMomentum();
const G4ThreeVector& pos = track->GetPosition();
ph->pex = mom.x();
ph->pey = mom.y();
ph->pez = mom.z();
ph->vex = pos.x();
ph->vey = pos.y();
ph->vez = pos.z();
// Set the simulator status bits
PropertyMask simStatus(m_currTrack.status);
// check if the last step ended on the worldVolume boundary
const G4Step* theLastStep = track->GetStep();
G4StepPoint* theLastPostStepPoint = NULL;
if(theLastStep) theLastPostStepPoint = theLastStep->GetPostStepPoint();
if( theLastPostStepPoint &&
( theLastPostStepPoint->GetStepStatus() == fWorldBoundary //particle left world volume
//|| theLastPostStepPoint->GetStepStatus() == fGeomBoundary
)
) {
simStatus.set(G4PARTICLE_SIM_LEFT_DETECTOR);
}
if(track->GetKineticEnergy() <= 0.) {
simStatus.set(G4PARTICLE_SIM_STOPPED);
}
/// Final update of the particle using the user handler
if ( m_userHandler ) {
m_userHandler->end(track, m_currTrack);
}
// These are candidate tracks with a probability to be stored due to their properties:
// - primary particle
// - hits created
// - secondaries
// - above energy threshold
// - to be kept due to creator process
// - to be kept due to user information of type 'Geant4ParticleInformation' stored in the G4Track
//
Geant4ParticleInformation* track_info =
dynamic_cast<Geant4ParticleInformation*>(track->GetUserInformation());
if ( !mask.isNull() || track_info ) {
m_equivalentTracks[g4_id] = g4_id;
ParticleMap::iterator ip = m_particleMap.find(g4_id);
if ( mask.isSet(G4PARTICLE_PRIMARY) ) {
ph.dump2(outputLevel()-1,name(),"Add Primary",h.id(),ip!=m_particleMap.end());
}
// Create a new MC particle from the current track information saved in the pre-tracking action
Particle* part = 0;
if ( ip==m_particleMap.end() ) part = m_particleMap[g4_id] = new Particle();
else part = (*ip).second;
if ( track_info ) {
mask.set(G4PARTICLE_KEEP_USER);
part->extension.reset(track_info->release());
}
part->get_data(m_currTrack);
}
else {
// These are tracks without any special properties.
//
// We will not store them on the record, but have to memorise the
// track identifier in order to restore the history for the created hits.
int pid = m_currTrack.g4Parent;
m_equivalentTracks[g4_id] = pid;
// Need to find the last stored particle and OR this particle's mask
// with the mask of the last stored particle
auto iend = m_equivalentTracks.end(), iequiv=m_equivalentTracks.end();
ParticleMap::iterator ip;
for(ip=m_particleMap.find(pid); ip == m_particleMap.end(); ip=m_particleMap.find(pid)) {
if (iequiv=m_equivalentTracks.find(pid); iequiv == iend) break; // ERROR
pid = (*iequiv).second;
}
if ( ip != m_particleMap.end() )
(*ip).second->reason |= track_reason;
else
ph.dumpWithVertex(outputLevel()+3,name(),"FATAL: No real particle parent present");
}
if(track->GetTrackStatus() == fSuspend) {
m_haveSuspended = true;
//track is already in particle map, we pick it up from there in begin again
if(m_particleMap.find(g4_id) != m_particleMap.end()) return;
//track is not already stored, keep it in special map
auto iPart = m_suspendedPM.emplace(g4_id, new Particle());
(iPart.first->second)->get_data(m_currTrack);
return; // we trust that we eventually return to this function with another status and go on then
}
}
/// Pre-event action callback
void Geant4ParticleHandler::beginEvent(const G4Event* event) {
Geant4PrimaryInteraction* interaction = context()->event().extension<Geant4PrimaryInteraction>();
info("+++ Event %d Begin event action. Access event related information.",event->GetEventID());
m_primaryMap = context()->event().extension<Geant4PrimaryMap>();
m_globalParticleID = interaction->nextPID();
m_particleMap.clear();
m_equivalentTracks.clear();
/// Call the user particle handler
if ( m_userHandler ) {
m_userHandler->begin(event);
}
}
/// Debugging: Dump Geant4 particle map
void Geant4ParticleHandler::dumpMap(const char* tag) const {
const string& n = name();
Geant4ParticleHandle::header4(INFO,n,tag);
for(ParticleMap::const_iterator iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i) {
Geant4ParticleHandle((*i).second).dump4(INFO,n,tag);
}
}
/// Post-event action callback
void Geant4ParticleHandler::endEvent(const G4Event* event) {
int count = 0;
int level = outputLevel();
do {
if ( level <= VERBOSE ) dumpMap("Particle ");
debug("+++ Iteration:%d Tracks:%d Equivalents:%d",++count,m_particleMap.size(),m_equivalentTracks.size());
} while( recombineParents() > 0 );
if ( level <= VERBOSE ) dumpMap( "Recombined");
// Rebase the simulated tracks, so that they fit to the generator particles
rebaseSimulatedTracks(0);
if ( level <= VERBOSE ) dumpMap( "Rebased ");
// Consistency check....
checkConsistency();
/// Call the user particle handler
if ( m_userHandler ) {
m_userHandler->end(event);
}
setVertexEndpointBit();
// Now export the data to the final record.
Geant4ParticleMap* part_map = context()->event().extension<Geant4ParticleMap>();
part_map->adopt(m_particleMap, m_equivalentTracks);
m_primaryMap = 0;
clear();
}
/// Rebase the simulated tracks, so that they fit to the generator particles
void Geant4ParticleHandler::rebaseSimulatedTracks(int ) {
/// No we have to update the map of equivalent tracks and assign the 'equivalentTrack' entry
TrackEquivalents equivalents, orgParticles;
ParticleMap finalParticles;
ParticleMap::const_iterator ipar, iend, i;
int count;
Geant4PrimaryInteraction* interaction = context()->event().extension<Geant4PrimaryInteraction>();
ParticleMap& pm = interaction->particles;
// (1.0) Copy the pre-defined particle mapping for the simulated tracks
// It is assumed the mapping is ZERO based without holes.
for(count = 0, iend=pm.end(), i=pm.begin(); i!=iend; ++i) {
Particle* p = (*i).second;
orgParticles[p->id] = p->id;
finalParticles[p->id] = p;
if ( p->id > count ) count = p->id;
if ( (p->reason&G4PARTICLE_PRIMARY) != G4PARTICLE_PRIMARY ) {
p->addRef();
}
}
// (1.1) Define the new particle mapping for the simulated tracks
for(++count, iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i) {
Particle* p = (*i).second;
if ( (p->reason&G4PARTICLE_PRIMARY) != G4PARTICLE_PRIMARY ) {
//if ( orgParticles.find(p->id) == orgParticles.end() ) {
orgParticles[p->id] = count;
finalParticles[count] = p;
p->id = count;
++count;
}
}
// (2) Re-evaluate the corresponding geant4 track equivalents using the new mapping
for(TrackEquivalents::iterator ie=m_equivalentTracks.begin(),ie_end=m_equivalentTracks.end(); ie!=ie_end; ++ie) {
int g4_equiv = (*ie).first;
while( (ipar=m_particleMap.find(g4_equiv)) == m_particleMap.end() ) {
TrackEquivalents::const_iterator iequiv = m_equivalentTracks.find(g4_equiv);
if ( iequiv == ie_end ) {
break; // ERROR !! Will be handled by printout below because ipar==end()
}
g4_equiv = (*iequiv).second;
}
TrackEquivalents::mapped_type equiv = (*ie).second;
if ( ipar != m_particleMap.end() ) {
Geant4ParticleHandle p = (*ipar).second;
equivalents[(*ie).first] = p->id; // requires (1) to be filled properly!
const G4ParticleDefinition* def = p.definition();
int pdg = int(fabs(def->GetPDGEncoding())+0.1);
if ( pdg != 0 && pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 ) {
error("+++ ERROR: Geant4 particle for track:%d last known is:%d -- is gluon or quark!",equiv,g4_equiv);
}
pdg = int(fabs(p->pdgID)+0.1);
if ( pdg != 0 && pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 ) {
error("+++ ERROR(2): Geant4 particle for track:%d last known is:%d -- is gluon or quark!",equiv,g4_equiv);
}
}
else {
error("+++ No Equivalent particle for track:%d last known is:%d",equiv,g4_equiv);
}
}
// (3) Compute the particle's parents and daughters.
// Replace the original Geant4 track with the
// equivalent particle still present in the record.
// Note:
// We rely here on the ordering of the particles accoding to their
// Processing by Geant4 to establish mother daughter relationships.
// == > use finalParticles map and NOT m_particleMap.
int equiv_id = -1;
for( auto& part : finalParticles ) {
auto& p = part.second;
if ( p->g4Parent > 0 ) {
TrackEquivalents::iterator iequ = equivalents.find(p->g4Parent);
if ( iequ != equivalents.end() ) {
equiv_id = (*iequ).second;//equivalents[p->g4Parent];
if ( (ipar=finalParticles.find(equiv_id)) != finalParticles.end() ) {
Particle* q = (*ipar).second;
bool prim = (p->reason&G4PARTICLE_PRIMARY) == G4PARTICLE_PRIMARY;
// We assume that the mother daughter relationship
// is filled by the event readers!
if ( !prim ) {
p->parents.insert(q->id);
}
if ( !p->parents.empty() ) {
int parent_id = (*p->parents.begin());
if ( parent_id == q->id )
q->daughters.insert(p->id);
else if ( !prim )
error("+++ Inconsistency in equivalent record! Parent: %d Daughter:%d",q->id, p->id);
}
else {
error("+++ Inconsistency in parent relashionship: %d NO parent!", p->id);
}
continue;
}
}
error("+++ Inconsistency in particle record: Geant4 parent %d "
"of particle %d not in record of final particles!",
p->g4Parent,p->id);
}
}
#if 0
for(iend=finalParticles.end(), i=finalParticles.begin(); i!=iend; ++i) {
Particle* p = (*i).second;
if ( p->g4Parent > 0 ) {
int parent_id = (*p->parents.begin());
if ( (ipar=finalParticles.find(parent_id)) != finalParticles.end() ) {
Particle* q = (*ipar).second;
// Generator particles have a proper history.
// We only deal with particles, which are not of MC origin.
//p->parents.insert(q->id);
if ( parent_id == q->id )
q->daughters.insert(p->id);
else
error("+++ Inconsistency in equivalent record! Parent: %d Daughter:%d",q->id, p->id);
continue;
}
error("+++ Inconsistency in particle record: Geant4 parent %d "
"of particle %d not in record of final particles!",
p->g4Parent,p->id);
}
}
#endif
m_equivalentTracks = equivalents;
m_particleMap = finalParticles;
}
/// Default callback to be answered if the particle should be kept if NO user handler is installed
bool Geant4ParticleHandler::defaultKeepParticle(Particle& particle) {
PropertyMask mask(particle.reason);
bool secondaries = mask.isSet(G4PARTICLE_HAS_SECONDARIES);
bool tracker_track = mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT);
bool calo_track = mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT);
bool hits_produced = mask.isSet(G4PARTICLE_CREATED_HIT);
bool low_energy = !mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD);
/// Remove this track if it has not created a hit and the energy is below threshold
if ( mask.isNull() || (secondaries && low_energy && !hits_produced) ) {
return true;
}
/// Remove this track if the energy is below threshold. Reassign hits to parent.
else if ( !hits_produced && low_energy ) {
return true;
}
/// Remove this track if the origine is in the calorimeter. Reassign hits to parent.
else if ( !tracker_track && calo_track && low_energy ) {
return true;
}
else {
//printout(INFO,name(),"+++ Track: %d should be kept for no obvious reason....",id);
}
return false;
}
/// Clean the monte carlo record. Remove all unwanted stuff.
/// This is the core of the object executed at the end of each event action.
int Geant4ParticleHandler::recombineParents() {
set<int> remove;
/// Need to start from BACK, to clean first the latest produced stuff.
for(ParticleMap::reverse_iterator i=m_particleMap.rbegin(); i!=m_particleMap.rend(); ++i) {
Particle* p = (*i).second;
PropertyMask mask(p->reason);
// Allow the user to force the particle handling either by
// or the reason mask with G4PARTICLE_KEEP_USER or
// to set the reason mask to NULL in order to drop it.
//
// If the mask entry is set to G4PARTICLE_FORCE_KILL
// or is set to NULL, the particle is ALWAYS removed
//
// Note: This may override all other decisions!
bool remove_me = m_userHandler ? m_userHandler->keepParticle(*p) : defaultKeepParticle(*p);
// Now look at the property mask of the particle
if ( mask.isNull() || mask.isSet(G4PARTICLE_FORCE_KILL) ) {
remove_me = true;
}
else if ( mask.isSet(G4PARTICLE_KEEP_USER) ) {
/// If user decides it must be kept, it MUST be kept!
mask.set(G4PARTICLE_KEEP_USER);
continue;
}
else if ( mask.isSet(G4PARTICLE_PRIMARY) ) {
/// Primary particles MUST be kept!
continue;
}
else if ( mask.isSet(G4PARTICLE_KEEP_ALWAYS) ) {
continue;
}
else if ( mask.isSet(G4PARTICLE_KEEP_PARENT) ) {
//continue;
}
else if ( mask.isSet(G4PARTICLE_KEEP_PROCESS) ) {
if(ParticleMap::iterator ip = m_particleMap.find(p->g4Parent); ip != m_particleMap.end() ) {
Particle* parent_part = (*ip).second;
PropertyMask parent_mask(parent_part->reason);
if ( parent_mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ) {
parent_mask.set(G4PARTICLE_KEEP_PARENT);
continue;
}
}
// Low energy stuff. Remove it. Reassign to parent.
//remove_me = true;
}
/// Remove this track from the list and also do the cleanup in the parent's children list
if ( remove_me ) {
int g4_id = (*i).first;
remove.insert(g4_id);
m_equivalentTracks[g4_id] = p->g4Parent;
if(ParticleMap::iterator ip = m_particleMap.find(p->g4Parent); ip != m_particleMap.end() ) {
Particle* parent_part = (*ip).second;
PropertyMask(parent_part->reason).set(mask.value());
parent_part->steps += p->steps;
parent_part->secondaries += p->secondaries;
/// Update of the particle using the user handler
if ( m_userHandler ) {
m_userHandler->combine(*p, *parent_part);
}
}
}
}
for( int r : remove ) {
if( auto ir = m_particleMap.find(r); ir != m_particleMap.end() ) {
(*ir).second->release();
m_particleMap.erase(ir);
}
}
return int(remove.size());
}
/// Check the record consistency
void Geant4ParticleHandler::checkConsistency() const {
int num_errors = 0;
/// First check the consistency of the particle map itself
for(const auto& part : m_particleMap ) {
Geant4Particle* particle = part.second;
Geant4ParticleHandle p(particle);
PropertyMask mask(p->reason);
PropertyMask status(p->status);
set<int>& daughters = p->daughters;
ParticleMap::const_iterator j;
// For all particles, the set of daughters must be contained in the record.
for( int id_dau : daughters ) {
if ( j=m_particleMap.find(id_dau); j == m_particleMap.end() ) {
++num_errors;
error("+++ Particle:%d Daughter %d is not in particle map!",p->id,id_dau);
}
}
// We assume that particles from the generator have consistent parents
// For all other particles except the primaries, the parent must be contained in the record.
if ( !mask.isSet(G4PARTICLE_PRIMARY) && !status.anySet(G4PARTICLE_GEN_STATUS) ) {
bool in_map = false, in_parent_list = false;
int parent_id = -1;
if( auto eq_it=m_equivalentTracks.find(p->g4Parent); eq_it != m_equivalentTracks.end() ) {
parent_id = (*eq_it).second;
in_map = (j=m_particleMap.find(parent_id)) != m_particleMap.end();
in_parent_list = p->parents.find(parent_id) != p->parents.end();
}
if ( !in_map || !in_parent_list ) {
char parent_list[1024];
parent_list[0] = 0;
++num_errors;
p.dumpWithMomentum(ERROR,name(),"INCONSISTENCY");
for( int ip : p->parents )
::snprintf(parent_list+strlen(parent_list),sizeof(parent_list)-strlen(parent_list),"%d ",ip);
error("+++ Particle:%d Parent %d (G4id:%d) In record:%s In parent list:%s [%s]",
p->id,parent_id,p->g4Parent,yes_no(in_map),yes_no(in_parent_list),parent_list);
}
}
}
if ( num_errors > 0 ) {
except("+++ Consistency check failed. Found %d problems.",num_errors);
}
}
void Geant4ParticleHandler::setVertexEndpointBit() {
for( auto& part : m_particleMap ) {
auto* p = part.second;
if( !p->parents.empty() ) {
Geant4Particle *parent(m_particleMap[ *p->parents.begin() ]);
const double X( parent->vex - p->vsx );
const double Y( parent->vey - p->vsy );
const double Z( parent->vez - p->vsz );
if( sqrt(X*X + Y*Y + Z*Z) > m_minDistToParentVertex ){
PropertyMask(p->status).set(G4PARTICLE_SIM_PARENT_RADIATED);
}
}
}
}