ATLAS Offline Software
ActsFatrasSimTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 #include <algorithm>
5 #include <random>
6 
7 #include "ActsFatrasSimTool.h"
8 #include "Acts/ActsVersion.hpp"
9 
10 #include "CLHEP/Random/RandFlat.h"
11 #include "CLHEP/Random/RandomEngine.h"
12 
15 
16 using namespace Acts::UnitLiterals;
17 
19  const std::string& name,
20  const IInterface* parent)
22 
24 
27  ATH_MSG_INFO("ISF::ActsFatrasSimTool update with ACTS version: v"
28  << Acts::VersionMajor << "." << Acts::VersionMinor << "."
29  << Acts::VersionPatch << " [" << Acts::CommitHash << "]");
30  // Retrieve particle filter
31  if (!m_particleFilter.empty()) ATH_CHECK(m_particleFilter.retrieve());
32 
33  // setup logger
34  m_logger = makeActsAthenaLogger(this, std::string("ActsFatras"),std::string("ActsFatrasSimTool"));
35 
36  // retrive tracking geo tool
37  ATH_CHECK(m_trackingGeometryTool.retrieve());
38  m_trackingGeometry = m_trackingGeometryTool->trackingGeometry();
39 
40  //retrive Magnetfield tool
41  ATH_MSG_VERBOSE("Using ATLAS magnetic field service");
42  ATH_CHECK( m_fieldCacheCondObjInputKey.initialize());
43 
44  // Random number service
45  if (m_rngSvc.retrieve().isFailure()) {
46  ATH_MSG_FATAL("Could not retrieve " << m_rngSvc);
47  return StatusCode::FAILURE;
48  }
49  // Get own engine with own seeds
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;
54  }
55 
56  // ISF truth service
57  ATH_CHECK (m_truthRecordSvc.retrieve());
58  ATH_MSG_DEBUG( "- Using ISF TruthRecordSvc : " << m_truthRecordSvc.typeAndName() );
59  return StatusCode::SUCCESS;
60 }
61 
63  ISFParticle& isp, ISFParticleContainer& secondaries,
64  McEventCollection* mcEventCollection) {
65  ATH_MSG_VERBOSE("Particle " << isp << " received for simulation.");
66  // Check if particle passes filter, if there is one
67  if (!m_particleFilter.empty() && !m_particleFilter->passFilter(isp)) {
68  ATH_MSG_VERBOSE("ISFParticle " << isp << " does not pass selection. Ignoring.");
69  return StatusCode::SUCCESS;
70  }
71  // Process ParticleState from particle stack
72  // Wrap the input ISFParticle in an STL vector with size of 1
73  const ISF::ISFParticleVector ispVector(1, &isp);
74  ATH_CHECK(this->simulateVector(ctx, ispVector, secondaries, mcEventCollection));
75  ATH_MSG_VERBOSE("Simulation done");
76  return StatusCode::SUCCESS;
77 }
78 
80  const EventContext& ctx,
82  ISFParticleContainer& secondaries,
83  McEventCollection* /*mcEventCollection*/, McEventCollection *) {
84 
85  m_randomEngine->setSeed(m_randomEngineName, ctx);
86  CLHEP::HepRandomEngine* randomEngine = m_randomEngine->getEngine(ctx);
87  Generator generator(CLHEP::RandFlat::shoot(randomEngine->flat()));
88  ATH_MSG_VERBOSE(name() << " RNG seed " << CLHEP::RandFlat::shoot(randomEngine->flat()));
89  ATH_MSG_VERBOSE(name() << " received vector of size "
90  << particles.size() << " particles for simulation.");
91 
92  // construct the ACTS simulator
93  Acts::Navigator navigator( Acts::Navigator::Config{ m_trackingGeometry }, m_logger);
94  auto bField = std::make_shared<ATLASMagneticFieldWrapper>();
95  auto chargedStepper = ChargedStepper(std::move(bField));
96  auto neutralStepper = NeutralStepper();
97  auto chargedPropagator = ChargedPropagator(chargedStepper, navigator, m_logger);
98  auto neutralPropagator = NeutralPropagator(neutralStepper, navigator, m_logger);
99  ChargedSimulation simulatorCharged(std::move(chargedPropagator), m_logger);
100  NeutralSimulation simulatorNeutral(std::move(neutralPropagator), m_logger);
101  Simulation simulator=Simulation(std::move(simulatorCharged),std::move(simulatorNeutral));
102  ATH_MSG_VERBOSE(name() << " Min pT for interaction " << m_interact_minPt * Acts::UnitConstants::MeV << " GeV");
103  // Acts propagater options
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;
112  // Create interaction list
113  simulator.charged.interactions = ActsFatras::makeStandardChargedElectroMagneticInteractions(m_interact_minPt * Acts::UnitConstants::MeV);
114  // get Geo and Mag map
115  ATH_MSG_VERBOSE(name() << " Getting per event Geo and Mag map");
116  Acts::MagneticFieldContext mctx = getMagneticFieldContext(ctx);
117  const ActsGeometryContext& gctx = m_trackingGeometryTool->getNominalGeometryContext();
118  auto anygctx = gctx.context();
119  // Loop over ISFParticleVector and process each separately
120  ATH_MSG_VERBOSE(name() << " Processing particles in ISFParticleVector.");
121  // For sihit creation
122  SiHitCollection pixelSiHits;
123  SiHitCollection sctSiHits;
124  for (const auto isfp : particles) {
125  // ====ACTSFatras Simulation====
126  // //
127  // input/output particle and hits containers
128  // Convert to ActsFatras::Particle
129  // ISF: Energy, mass, and momentum are in MeV, position in mm
130  // Acts: Energy, mass, and momentum are in GeV, position in mm
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()),
134  isfp->charge(),isfp->mass() * Acts::UnitConstants::MeV)
135  .setDirection(Acts::makeDirectionFromPhiEta(isfp->momentum().phi(), isfp->momentum().eta()))
136  .setAbsoluteMomentum(isfp->momentum().mag() * Acts::UnitConstants::MeV)
137  .setPosition4(ActsTrk::convertPosToActs(isfp->position(), isfp->timeStamp()))};
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;
142  // simulate
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;
151  }
152  }
153 
154  ATH_MSG_DEBUG(name() << " initial particle " << simulatedInitial[0]);
155  ATH_MSG_DEBUG(name() << " ActsFatras simulator hits: " << hits.size());
156  int i = 0;
157  for (const auto& hit : hits) {
158  ATH_MSG_DEBUG(name() << " hit pos: " << hit.position() );
159  ++i;
160  if (i>5) break;
161  }
162  ATH_MSG_DEBUG(name() << " No. of particles after ActsFatras simulator: " << simulatedFinal.size());
163  if (!simulatedFinal.empty()){
164  ATH_MSG_DEBUG(name() << " start procesing secondaries");
165  auto itr = simulatedFinal.begin();
166  // Save hits of isfp
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);
170  // Process secondaries
171  auto isKilled = !itr->isAlive();
172  int maxGeneration = (simulatedFinal.back()).particleId().generation();
173  ATH_MSG_DEBUG(name() << " maxGeneration: "<< maxGeneration);
174  for (int gen = 0; gen <= maxGeneration; ++gen){
175  ATH_MSG_DEBUG(name() << " start with genration "<< gen << "|" << maxGeneration << ": "<< *itr);
176  auto vecsecisfp = std::make_unique<ISF::ISFParticleVector>();
177  while (static_cast<int>(itr->particleId().generation()) == gen){
178  ATH_MSG_DEBUG(name() << " genration "<< gen << "|" << maxGeneration << ": "<< *itr);
179  if(itr->isSecondary()){
180  // convert final particles to ISF::particle
181  const auto pos = ActsTrk::convertPosFromActs(itr->fourPosition()).first;
182  const auto mom = ActsTrk::convertMomFromActs(itr->fourMomentum()).first;
183  double mass = itr->mass() / Acts::UnitConstants::MeV;
184  double charge = itr->charge();
185  int pdgid = itr->pdg();
186  auto properTime = ActsTrk::timeToAthena(itr->time());
187  const int status = 1 + HepMC::SIM_STATUS_THRESHOLD;
188  const int id = HepMC::UNDEFINED_ID;
189  auto secisfp = new ISF::ISFParticle (pos,mom,mass,charge,pdgid,status,properTime,*isfp,id);
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);
193  }
194  else{
195  ATH_MSG_DEBUG(name() <<" primary particle found with generation "<< gen);
196  }
197  ++itr;
198  }
199  if (!vecsecisfp->empty()) {
200  ISF::ISFTruthIncident truth(*isfp,
201  *vecsecisfp,
202  getATLASProcessCode((itr-1)->process()),
203  isfp->nextGeoID(),
204  isKilled&&gen==maxGeneration ? ISF::fKillsPrimary : ISF::fPrimarySurvives
205  );
206  ATH_MSG_DEBUG(name() << " Truth incident parentPt2(MinPt2) " << truth.parentPt2() <<" (100 MeV)");
207  ATH_MSG_DEBUG(name() << " Truth incident ChildPt2(MinPt2) " << truth.childrenPt2Pass(300) <<" (300 MeV)");
208  m_truthRecordSvc->registerTruthIncident(truth, true);
211  for (auto *secisfp : *vecsecisfp){
212  if (secisfp->getTruthBinding()) {
213  secondaries.push_back(secisfp);
214  }
215  else {
216  ATH_MSG_WARNING("Secondary particle not written out to truth.\n Parent (" << isfp << ")\n Secondary (" << *secisfp <<")");
217  }
218  } // end of truth binding
219  }// end of store truth bind secondaries
220  }
221  }// end of secondaries
222  ATH_MSG_VERBOSE(name() << " No. of secondaries: " << secondaries.size());
223  ATH_MSG_DEBUG(name() << " End of particle " << isfp->barcode());
224  m_ActsFatrasWriteHandler->createHits(*isfp, m_trackingGeometry,hits,pixelSiHits,sctSiHits);
225 
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);
230  } // end of isfp loop
231  return StatusCode::SUCCESS;
232 }
233 
234 Acts::MagneticFieldContext ISF::ActsFatrasSimTool::getMagneticFieldContext(const EventContext& ctx) const {
235  SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCacheCondObjInputKey, ctx};
236  if (!readHandle.isValid()) {
237  ATH_MSG_ERROR(name() + ": Failed to retrieve magnetic field condition data " + m_fieldCacheCondObjInputKey.key() + ".");
238  }
239  else ATH_MSG_DEBUG(name() << "retrieved magnetic field condition data "<< m_fieldCacheCondObjInputKey.key());
240  const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
241 
242  return Acts::MagneticFieldContext(fieldCondObj);
243 }
ISF::ISFParticleContainer
std::list< ISF::ISFParticle * > ISFParticleContainer
generic ISFParticle container (not necessarily a std::list!)
Definition: ISFParticleContainer.h:23
ISF::ActsFatrasSimTool::ChargedPropagator
Acts::Propagator< ChargedStepper, Navigator > ChargedPropagator
Definition: ActsFatrasSimTool.h:179
ISF::ActsFatrasSimTool::Generator
std::ranlux48 Generator
Definition: ActsFatrasSimTool.h:174
ISF::ISFTruthIncident::updateChildParticleProperties
void updateChildParticleProperties()
Update the id and particleLink properties of the child particles (to be called after registerTruthInc...
Definition: ISFTruthIncident.cxx:239
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
ISF::ActsFatrasSimTool::~ActsFatrasSimTool
virtual ~ActsFatrasSimTool()
Definition: ActsFatrasSimTool.cxx:23
TRTCalib_Extractor.hits
hits
Definition: TRTCalib_Extractor.py:35
get_generator_info.result
result
Definition: get_generator_info.py:21
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
AtlasFieldCacheCondObj
Definition: AtlasFieldCacheCondObj.h:19
ISF::ActsFatrasSimTool::NeutralPropagator
Acts::Propagator< NeutralStepper, Navigator > NeutralPropagator
Definition: ActsFatrasSimTool.h:182
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
ISF::ActsFatrasSimTool::ChargedStepper
Acts::EigenStepper< Acts::EigenStepperDefaultExtension > ChargedStepper
Definition: ActsFatrasSimTool.h:178
ActsTrk::detail::Navigator
Acts::Navigator Navigator
Definition: Tracking/Acts/ActsTrackReconstruction/src/detail/Definitions.h:31
AtlasHitsVector< SiHit >
McEventCollection
McEventCollection
Definition: GeneratorObjectsTPCnv.cxx:60
ISF::ISFParticle
Definition: ISFParticle.h:42
ISFTruthIncident.h
ActsTrk::convertPosToActs
Acts::Vector4 convertPosToActs(const Amg::Vector3D &athenaPos, const double athenaTime=0.)
Converts a position vector & time from Athena units into Acts units.
Definition: UnitConverters.h:74
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
SUSY_SimplifiedModel_PostInclude.process
string process
Definition: SUSY_SimplifiedModel_PostInclude.py:43
ActsTrk::timeToAthena
constexpr double timeToAthena(const double actsT)
Converts a time unit from Acts to Athena units.
Definition: UnitConverters.h:56
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
xAOD::Particle
Particle_v1 Particle
Define the latest version of the particle class.
Definition: Event/xAOD/xAODParticleEvent/xAODParticleEvent/Particle.h:17
master.gen
gen
Definition: master.py:32
python.SystemOfUnits.MeV
float MeV
Definition: SystemOfUnits.py:172
makeActsAthenaLogger
std::unique_ptr< const Acts::Logger > makeActsAthenaLogger(IMessageSvc *svc, const std::string &name, int level, std::optional< std::string > parent_name)
Definition: Tracking/Acts/ActsInterop/src/Logger.cxx:64
ActsGeometryContext::context
Acts::GeometryContext context() const
Definition: ActsGeometryContext.h:45
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
ISF::ActsFatrasSimTool::simulateVector
virtual StatusCode simulateVector(const EventContext &ctx, const ISFParticleVector &particles, ISFParticleContainer &secondaries, McEventCollection *mcEventCollection, McEventCollection *shadowTruth=nullptr) override
Simulation call for vectors of particles.
Definition: ActsFatrasSimTool.cxx:79
ISF::ISFTruthIncident
Definition: ISFTruthIncident.h:35
lumiFormat.i
int i
Definition: lumiFormat.py:85
ISF::ActsFatrasSimTool::initialize
virtual StatusCode initialize() override
Definition: ActsFatrasSimTool.cxx:25
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
ISF::ISFParticleVector
std::vector< ISF::ISFParticle * > ISFParticleVector
ISFParticle vector.
Definition: ISFParticleContainer.h:26
ActsTrk::convertPosFromActs
std::pair< Amg::Vector3D, double > convertPosFromActs(const Acts::Vector4 &actsPos)
Converts an Acts 4-vector into a pair of an Athena spatial vector and the passed time.
Definition: UnitConverters.h:86
test_pyathena.parent
parent
Definition: test_pyathena.py:15
ISF::ActsFatrasSimTool::SingleParticleSimulation
Single particle simulation with fixed propagator, interactions, and decay.
Definition: ActsFatrasSimTool.h:90
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
ActsTrk::convertMomFromActs
std::pair< Amg::Vector3D, double > convertMomFromActs(const Acts::Vector4 &actsMom)
Converts an Acts four-momentum vector into an pair of an Athena three-momentum and the paritcle's ene...
Definition: UnitConverters.h:109
ISF::ISFTruthIncident::parentPt2
double parentPt2() const override final
Return pT^2 of the parent particle.
Definition: ISFTruthIncident.cxx:81
ISF::ActsFatrasSimTool::NeutralStepper
Acts::StraightLineStepper NeutralStepper
Definition: ActsFatrasSimTool.h:181
ISF::ITruthIncident::childrenPt2Pass
bool childrenPt2Pass(double pt2cut)
Return true if at least one child particle passes the given pT^2 cut (= at least one child with pT^2 ...
Definition: ITruthIncident.h:158
ISF::ActsFatrasSimTool::ActsFatrasSimTool
ActsFatrasSimTool(const std::string &type, const std::string &name, const IInterface *parent)
Definition: ActsFatrasSimTool.cxx:18
ActsGeometryContext
Include the GeoPrimitives which need to be put first.
Definition: ActsGeometryContext.h:27
HepMC::UNDEFINED_ID
constexpr int UNDEFINED_ID
Definition: MagicNumbers.h:56
HepMC::SIM_STATUS_THRESHOLD
constexpr int SIM_STATUS_THRESHOLD
Constant definiting the status threshold for simulated particles, eg. can be used to separate generat...
Definition: MagicNumbers.h:46
ISF::BaseSimulatorTool
Definition: BaseSimulatorTool.h:36
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
MagicNumbers.h
charge
double charge(const T &p)
Definition: AtlasPID.h:986
ISF::ActsFatrasSimTool::getMagneticFieldContext
virtual Acts::MagneticFieldContext getMagneticFieldContext(const EventContext &) const
Definition: ActsFatrasSimTool.cxx:234
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
ISF::ActsFatrasSimTool::simulate
virtual StatusCode simulate(const EventContext &ctx, ISFParticle &isp, ISFParticleContainer &, McEventCollection *) override
Definition: ActsFatrasSimTool.cxx:62
ISF::BaseSimulatorTool::initialize
virtual StatusCode initialize() override
Definition: BaseSimulatorTool.h:57
mc.generator
generator
Configure Herwig7 These are the commands corresponding to what would go into the regular Herwig infil...
Definition: mc.MGH7_FxFx_H71-DEFAULT_test.py:18
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
ISF::ISFTruthIncident::updateParentAfterIncidentProperties
void updateParentAfterIncidentProperties()
Update the id and particleLink properties of the parentAfterIncident (to be called after registerTrut...
Definition: ISFTruthIncident.cxx:211
LArG4FSStartPointFilter.particles
list particles
Definition: LArG4FSStartPointFilter.py:84
Simulation
Definition: BeamEffectsAlg.cxx:21
TileCalibBlobPython_writeDefaultCs.Simulation
Simulation
Definition: TileCalibBlobPython_writeDefaultCs.py:12
merge.status
status
Definition: merge.py:16
calibdata.copy
bool copy
Definition: calibdata.py:26
ISF::fPrimarySurvives
@ fPrimarySurvives
Definition: ISFTruthIncident.h:24
ActsFatrasSimTool.h
ISF::fKillsPrimary
@ fKillsPrimary
Definition: ISFTruthIncident.h:25