8 #include "Acts/ActsVersion.hpp"
10 #include "CLHEP/Random/RandFlat.h"
11 #include "CLHEP/Random/RandomEngine.h"
16 using namespace Acts::UnitLiterals;
19 const std::string&
name,
27 ATH_MSG_INFO(
"ISF::ActsFatrasSimTool update with ACTS version: v"
28 << Acts::VersionMajor <<
"." << Acts::VersionMinor <<
"."
29 << Acts::VersionPatch <<
" [" << Acts::CommitHash <<
"]");
31 if (!m_particleFilter.empty())
ATH_CHECK(m_particleFilter.retrieve());
34 m_logger =
makeActsAthenaLogger(
this, std::string(
"ActsFatras"),std::string(
"ActsFatrasSimTool"));
37 ATH_CHECK(m_trackingGeometryTool.retrieve());
38 m_trackingGeometry = m_trackingGeometryTool->trackingGeometry();
42 ATH_CHECK( m_fieldCacheCondObjInputKey.initialize());
45 if (m_rngSvc.retrieve().isFailure()) {
47 return StatusCode::FAILURE;
50 m_randomEngine = m_rngSvc->getEngine(
this, m_randomEngineName.value());
51 if (!m_randomEngine) {
52 ATH_MSG_FATAL(
"Could not get random engine '" << m_randomEngineName.value() <<
"'");
53 return StatusCode::FAILURE;
58 ATH_MSG_DEBUG(
"- Using ISF TruthRecordSvc : " << m_truthRecordSvc.typeAndName() );
59 return StatusCode::SUCCESS;
67 if (!m_particleFilter.empty() && !m_particleFilter->passFilter(isp)) {
68 ATH_MSG_VERBOSE(
"ISFParticle " << isp <<
" does not pass selection. Ignoring.");
69 return StatusCode::SUCCESS;
74 ATH_CHECK(this->simulateVector(ctx, ispVector, secondaries, mcEventCollection));
76 return StatusCode::SUCCESS;
80 const EventContext& ctx,
85 m_randomEngine->setSeed(m_randomEngineName, ctx);
86 CLHEP::HepRandomEngine* randomEngine = m_randomEngine->getEngine(ctx);
90 <<
particles.size() <<
" particles for simulation.");
93 Acts::Navigator navigator( Acts::Navigator::Config{ m_trackingGeometry }, m_logger);
94 auto bField = std::make_shared<ATLASMagneticFieldWrapper>();
104 simulator.charged.maxStepSize = m_maxStepSize;
105 simulator.charged.maxStep = m_maxStep;
106 simulator.charged.pathLimit = m_pathLimit;
107 simulator.charged.maxRungeKuttaStepTrials = m_maxRungeKuttaStepTrials;
108 simulator.charged.loopProtection = m_loopProtection;
109 simulator.charged.loopFraction = m_loopFraction;
110 simulator.charged.targetTolerance = m_tolerance;
111 simulator.charged.stepSizeCutOff = m_stepSizeCutOff;
113 simulator.charged.interactions = ActsFatras::makeStandardChargedElectroMagneticInteractions(m_interact_minPt *
Acts::UnitConstants::MeV);
116 Acts::MagneticFieldContext mctx = getMagneticFieldContext(ctx);
131 ATH_MSG_DEBUG(
name() <<
" Convert ISF::Particle(mass) " << isfp->id()<<
"|" << isfp<<
"(" << isfp->mass() <<
")");
132 std::vector<ActsFatras::Particle> input = std::vector<ActsFatras::Particle>{
133 ActsFatras::Particle(ActsFatras::Barcode().setVertexPrimary(0).setParticle(isfp->id()),
static_cast<Acts::PdgParticle
>(isfp->pdgCode()),
135 .setDirection(Acts::makeDirectionFromPhiEta(isfp->momentum().phi(), isfp->momentum().eta()))
138 ATH_MSG_DEBUG(
name() <<
" Propagating ActsFatras::Particle vertex|particle|generation|subparticle, " << input[0]);
139 std::vector<ActsFatras::Particle> simulatedInitial;
140 std::vector<ActsFatras::Particle> simulatedFinal;
141 std::vector<ActsFatras::Hit>
hits;
143 auto result=simulator.simulate(anygctx, mctx,
generator, input, simulatedInitial, simulatedFinal,
hits);
144 auto simulatedFailure=
result.value();
145 if (simulatedFailure.size()>0){
146 for (
const auto& simfail : simulatedFailure){
147 auto errCode = Acts::make_error_code(Acts::PropagatorError(simfail.error.value()));
148 ATH_MSG_WARNING(
name() <<
" Particle id " <<simfail.particle.particleId()<<
": fail to be simulated during Propagation: " << errCode.message());
149 ATH_MSG_WARNING(
name() <<
" Particle vertex|particle|generation|subparticle"<<simfail.particle <<
" starts from position" <<
Acts::toString(simfail.particle.position()) <<
" and direction " <<
Acts::toString(simfail.particle.direction()));
150 return StatusCode::SUCCESS;
157 for (
const auto& hit :
hits) {
162 ATH_MSG_DEBUG(
name() <<
" No. of particles after ActsFatras simulator: " << simulatedFinal.size());
163 if (!simulatedFinal.empty()){
165 auto itr = simulatedFinal.begin();
167 std::vector<ActsFatras::Hit> particle_hits;
168 std::copy(
hits.begin(),
hits.begin()+itr->numberOfHits(), std::back_inserter(particle_hits));
169 m_ActsFatrasWriteHandler->createHits(*isfp, m_trackingGeometry,particle_hits,m_pixelSiHits,m_sctSiHits);
171 auto isKilled = !itr->isAlive();
172 int maxGeneration = (simulatedFinal.back()).particleId().generation();
174 for (
int gen = 0;
gen <= maxGeneration; ++
gen){
176 auto vecsecisfp = std::make_unique<ISF::ISFParticleVector>();
177 while (
static_cast<int>(itr->particleId().generation()) ==
gen){
179 if(itr->isSecondary()){
184 double charge = itr->charge();
185 int pdgid = itr->pdg();
190 ATH_MSG_DEBUG(
name() <<
" secondaries particle (ACTS): "<<*itr<<
"("<<itr->momentum()<<
")|time "<<itr->time()<<
"|process "<< getATLASProcessCode(itr->process()));
191 ATH_MSG_DEBUG(
name() <<
" secondaries particle (ISF): " << *secisfp <<
" time "<<secisfp->timeStamp());
192 vecsecisfp->push_back(secisfp);
199 if (!vecsecisfp->empty()) {
202 getATLASProcessCode((itr-1)->
process()),
208 m_truthRecordSvc->registerTruthIncident(truth,
true);
211 for (
auto *secisfp : *vecsecisfp){
212 if (secisfp->getTruthBinding()) {
213 secondaries.push_back(secisfp);
216 ATH_MSG_WARNING(
"Secondary particle not written out to truth.\n Parent (" << isfp <<
")\n Secondary (" << *secisfp <<
")");
224 m_ActsFatrasWriteHandler->createHits(*isfp, m_trackingGeometry,
hits,pixelSiHits,sctSiHits);
226 std::vector<ActsFatras::Particle>().swap(input);
227 std::vector<ActsFatras::Particle>().swap(simulatedInitial);
228 std::vector<ActsFatras::Particle>().swap(simulatedFinal);
229 std::vector<ActsFatras::Hit>().swap(
hits);
231 return StatusCode::SUCCESS;
236 if (!readHandle.isValid()) {
237 ATH_MSG_ERROR(
name() +
": Failed to retrieve magnetic field condition data " + m_fieldCacheCondObjInputKey.key() +
".");
239 else ATH_MSG_DEBUG(
name() <<
"retrieved magnetic field condition data "<< m_fieldCacheCondObjInputKey.key());
242 return Acts::MagneticFieldContext(fieldCondObj);