95 {
96
97 std::vector<ISFParticle*> selectedParticles;
98 for (const auto isfp : particles) {
100 ATH_MSG_VERBOSE(
"ISFParticle " << *isfp <<
" does not pass selection. Ignoring.");
101 continue;
102 }
103 selectedParticles.push_back(isfp);
104 }
105 if (selectedParticles.empty()) {
107 return StatusCode::SUCCESS;
108 }
109
111 CLHEP::HepRandomEngine* randomEngine = m_randomEngine->getEngine(ctx);
115 << selectedParticles.size() << " particles for simulation.");
116
117
119 auto bField = std::make_shared<ATLASMagneticFieldWrapper>();
128
137
138 simulator.charged.interactions = ActsFatras::makeStandardChargedElectroMagneticInteractions(
m_interact_minPt * Acts::UnitConstants::MeV);
139
140
142
143
147 auto anygctx = gctx.context();
148
150 for (const auto isfp : selectedParticles) {
151
152
153
154
155
156
157 ATH_MSG_DEBUG(
name() <<
" Convert ISF::Particle(mass) " << isfp->id()<<
"|" << *isfp<<
"(" << isfp->mass() <<
")");
158 std::vector<ActsFatras::Particle>
input = std::vector<ActsFatras::Particle>{
159 ActsFatras::Particle(ActsFatras::Barcode().withVertexPrimary(0).withParticle(isfp->id()), static_cast<Acts::PdgParticle>(isfp->pdgCode()),
160 isfp->charge(),isfp->mass() * Acts::UnitConstants::MeV)
161 .setDirection(Acts::makeDirectionFromPhiEta(isfp->momentum().phi(), isfp->momentum().eta()))
162 .setAbsoluteMomentum(isfp->momentum().mag() * Acts::UnitConstants::MeV)
164 ATH_MSG_DEBUG(
name() <<
" Propagating ActsFatras::Particle vertex|particle|generation|subparticle, " << input[0]);
165 std::vector<ActsFatras::Particle> simulatedInitial;
166 std::vector<ActsFatras::Particle> simulatedFinal;
167 std::vector<ActsFatras::Hit>
hits;
168
169 auto result=simulator.simulate(anygctx, mctx, generator, input, simulatedInitial, simulatedFinal, hits);
170 auto simulatedFailure=
result.value();
171 if (simulatedFailure.size()>0){
172 for (const auto& simfail : simulatedFailure){
173 auto errCode = Acts::make_error_code(Acts::PropagatorError(simfail.error.value()));
174 ATH_MSG_WARNING(
name() <<
" Particle id " <<simfail.particle.particleId()<<
": fail to be simulated during Propagation: " << errCode.message());
175 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()));
176 return StatusCode::SUCCESS;
177 }
178 }
179
183 for (const auto& hit : hits) {
186 if (i>5) break;
187 }
188 ATH_MSG_DEBUG(
name() <<
" No. of particles after ActsFatras simulator: " << simulatedFinal.size());
189 if (!simulatedFinal.empty()){
191 auto itr = simulatedFinal.begin();
192
193 std::vector<ActsFatras::Hit> particle_hits;
194 if (itr->numberOfHits() > 0) {
195 std::copy(
hits.begin(),
hits.begin()+itr->numberOfHits(), std::back_inserter(particle_hits));
197 }
198
199 auto isKilled = !itr->isAlive();
200 int maxGeneration = simulatedFinal.back().particleId().generation();
202 for (
int gen = 0;
gen <= maxGeneration; ++
gen){
203 ATH_MSG_DEBUG(
name() <<
" start with generation "<< gen <<
"|" << maxGeneration <<
": "<< *itr);
204 auto vecsecisfp = std::make_unique<ISF::ISFParticleVector>();
205 std::unique_ptr<ISF::ISFParticle> newisfp = nullptr;
206
207 while (itr != simulatedFinal.end() && static_cast<int>(itr->particleId().generation()) == gen) {
208 ATH_MSG_DEBUG(
name() <<
" genration "<< gen <<
"|" << maxGeneration <<
": "<< *itr);
209 if(itr->isSecondary()){
210
213 double mass = itr->mass() / Acts::UnitConstants::MeV;
214 double charge = itr->charge();
215 int pdgid = itr->pdg();
219 auto secisfp = std::make_unique<ISF::ISFParticle>(pos,mom,mass,
charge,pdgid,status,properTime,*isfp,
id);
220 secisfp->setNextGeoID(
m_geoIDSvc->identifyNextGeoID(*secisfp));
223 << " pos=" << secisfp->position() << " mom=" << secisfp->momentum()
224 <<
" GeoID=" <<
m_geoIDSvc->identifyNextGeoID(*secisfp));
225 vecsecisfp->push_back(secisfp.release());
226 }
227 else{
228
230
231 if (!isKilled) {
232 auto fisfp = std::make_unique<ISF::ISFParticle>(*isfp);
235 ATH_MSG_DEBUG(
name() <<
" After simulation, primary particle state: " << *fisfp);
237 ATH_MSG_VERBOSE(
"ISFParticle" << fisfp <<
" after simulation does not pass selection. Ignoring for boundary check.");
238 continue;
239 }
241 <<
" new particle GeoID: " <<
m_geoIDSvc->identifyGeoID(*fisfp)
242 <<
", nextGeoID: " <<
m_geoIDSvc->identifyNextGeoID(*fisfp));
243
244
246
247
248 Acts::BoundTrackParameters startParams = Acts::BoundTrackParameters::createCurvilinear(
249 itr->fourPosition(), itr->direction(), itr->qOverP(), std::nullopt, itr->hypothesis());
250
251
252 auto do_exit_startsurface =
checkStartSurface(mctx, anygctx, surfaceCheckPropagator, startParams);
253 ATH_MSG_DEBUG(
name() <<
" checkStartSurface returned: " << do_exit_startsurface);
254 if (do_exit_startsurface) {
255 ATH_MSG_DEBUG(
name() <<
" Particle starts at a valid surface, doing extrapolation...");
256
257
260 try {
261 auto stepsResult =
m_extrapolationTool->propagationSteps(ctx, startParams, Acts::Direction::Forward());
262 auto steps = stepsResult.value().first;
264 if (
steps.size() != 0) {
265 for (const auto& step : steps) {
267 <<
" (eta " << Acts::VectorHelpers::eta(
step.position)
268 <<
") with GeoID " <<
step.geoID);
270 nextGeoID =
m_geoIDSvc->identifyGeoID(entryPos);
274 break;
275 }
276 }
277 } else {
279 }
280 }
281 catch (const std::exception& e) {
283 break;
284 }
285
288 double mass = itr->mass() / Acts::UnitConstants::MeV;
289 double charge = itr->charge();
290 int pdgid = itr->pdg();
292
293
294 newisfp = std::make_unique<ISF::ISFParticle>(entryPos, mom, mass,
charge, pdgid, isfp->status(), properTime, *isfp, isfp->id(), isfp->barcode());
295 newisfp->setNextGeoID(nextGeoID);
296 ATH_MSG_DEBUG(
name() <<
" Truthbinding of parent ISFParticle: " << (isfp->getTruthBinding() ?
"exists" :
"null"));
297 if (isfp->getTruthBinding()) {
298 ATH_MSG_DEBUG(
name() <<
" Current GenParticle: " << isfp->getTruthBinding()->getCurrentGenParticle());
299 }
302 << "(" << nextGeoID << ")");
303 }
304
305
309 << "(" << newisfp->nextGeoID() << ")");
310
311
313
314 switch(nextGeoID) {
318 break;
322 break;
323 default:
325 break;
326 }
327
329 vecsecisfp->push_back(newisfp.release());
330 } else {
332 }
333 }
334 }
335 else {
336 ATH_MSG_DEBUG(
name() <<
" No starting surface found, skipping boundary check and extrapolation.");
337 }
338 }
339 }
340 ++itr;
341 }
342
343
344 if (!vecsecisfp->empty()) {
345
346 int processCode = 0;
349
350 if (newisfp) {
351
352 processCode = 91;
353 geoID = newisfp->nextGeoID();
355 } else {
356
361 }
362
363 ISF::ISFTruthIncident truth(*isfp,
364 *vecsecisfp,
365 processCode,
366 geoID,
367 isParentKilled);
368
369 ATH_MSG_DEBUG(
name() <<
" Truth incident parentPt2(MinPt2) " << truth.parentPt2() <<
" (100 MeV)");
370 ATH_MSG_DEBUG(
name() <<
" Truth incident ChildPt2(MinPt2) " << truth.childrenPt2Pass(300) <<
" (300 MeV)");
372 truth.updateParentAfterIncidentProperties();
373 truth.updateChildParticleProperties();
374 for (auto *secisfp : *vecsecisfp){
375 if (secisfp->getTruthBinding()) {
376 secondaries.push_back(secisfp);
378 << *isfp << ")\n Secondary (" << *secisfp <<")");
379
380 ATH_MSG_DEBUG(
"Secondary particle push back to ISF, TruthBinding: " << secisfp->getTruthBinding()->getCurrentGenParticle() <<
" (current) | " << secisfp->getTruthBinding()->getPrimaryGenParticle() <<
" (primary) | " << secisfp->getTruthBinding()->getGenerationZeroGenParticle() <<
" (zero)");
381 if (secisfp->getTruthBinding()->getCurrentGenParticle() !=
nullptr)
ATH_MSG_DEBUG(
"Secondary particle GenParticle EndVertex: " << (secisfp->getTruthBinding()->getCurrentGenParticle()->end_vertex() ?
HepMC::barcode(secisfp->getTruthBinding()->getCurrentGenParticle()->end_vertex()) : 1));
382 } else {
383 ATH_MSG_WARNING(
"Secondary particle not written out to truth.\n Parent ("
384 << *isfp << ")\n Secondary (" << *secisfp <<")");
385 delete secisfp;
386 }
387 }
388 }
389 }
390 }
393
394 std::vector<ActsFatras::Particle>().swap(input);
395 std::vector<ActsFatras::Particle>().swap(simulatedInitial);
396 std::vector<ActsFatras::Particle>().swap(simulatedFinal);
397 std::vector<ActsFatras::Hit>().swap(hits);
398 }
399 return StatusCode::SUCCESS;
400}
double charge(const T &p)
static const char * getName(int region)
const std::string process
constexpr double timeToAthena(const double actsT)
Converts a time unit from Acts to Athena units.
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...
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.
Acts::Vector4 convertPosToActs(const Amg::Vector3D &athenaPos, const double athenaTime=0.)
Converts a position vector & time from Athena units into Acts units.
AtlasRegion
A simple enum of ATLAS regions and sub-detectors.
constexpr int UNDEFINED_ID
constexpr int SIM_STATUS_THRESHOLD
Constant definiting the status threshold for simulated particles, eg. can be used to separate generat...
EntryLayer
Identifiers for the TrackRecordCollections on the boundaries between CaloEntry: Inner Detector - Calo...
generator
Configure Herwig7 These are the commands corresponding to what would go into the regular Herwig infil...