83 {
84
86 CLHEP::HepRandomEngine* randomEngine = m_randomEngine->getEngine(ctx);
90 <<
particles.size() <<
" particles for simulation.");
91
92
94 auto bField = std::make_shared<ATLASMagneticFieldWrapper>();
103
112
113 simulator.charged.interactions = ActsFatras::makeStandardChargedElectroMagneticInteractions(
m_interact_minPt * Acts::UnitConstants::MeV);
114
119
121 for (const auto isfp : particles) {
122
123
124
125
126
127
128 ATH_MSG_DEBUG(
name() <<
" Convert ISF::Particle(mass) " << isfp->id()<<
"|" << isfp<<
"(" << isfp->mass() <<
")");
129 std::vector<ActsFatras::Particle>
input = std::vector<ActsFatras::Particle>{
130 ActsFatras::Particle(ActsFatras::Barcode().withVertexPrimary(0).withParticle(isfp->id()), static_cast<Acts::PdgParticle>(isfp->pdgCode()),
131 isfp->charge(),isfp->mass() * Acts::UnitConstants::MeV)
132 .setDirection(Acts::makeDirectionFromPhiEta(isfp->momentum().phi(), isfp->momentum().eta()))
133 .setAbsoluteMomentum(isfp->momentum().mag() * Acts::UnitConstants::MeV)
135 ATH_MSG_DEBUG(
name() <<
" Propagating ActsFatras::Particle vertex|particle|generation|subparticle, " << input[0]);
136 std::vector<ActsFatras::Particle> simulatedInitial;
137 std::vector<ActsFatras::Particle> simulatedFinal;
138 std::vector<ActsFatras::Hit>
hits;
139
140 auto result=simulator.simulate(anygctx, mctx, generator, input, simulatedInitial, simulatedFinal, hits);
141 auto simulatedFailure=
result.value();
142 if (simulatedFailure.size()>0){
143 for (const auto& simfail : simulatedFailure){
144 auto errCode = Acts::make_error_code(Acts::PropagatorError(simfail.error.value()));
145 ATH_MSG_WARNING(
name() <<
" Particle id " <<simfail.particle.particleId()<<
": fail to be simulated during Propagation: " << errCode.message());
146 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()));
147 return StatusCode::SUCCESS;
148 }
149 }
150
154 for (const auto& hit : hits) {
157 if (i>5) break;
158 }
159 ATH_MSG_DEBUG(
name() <<
" No. of particles after ActsFatras simulator: " << simulatedFinal.size());
160 if (!simulatedFinal.empty()){
162 auto itr = simulatedFinal.begin();
163
164 std::vector<ActsFatras::Hit> particle_hits;
165 if (itr->numberOfHits() > 0) {
166 std::copy(
hits.begin(),
hits.begin()+itr->numberOfHits(), std::back_inserter(particle_hits));
168 }
169
170 auto isKilled = !itr->isAlive();
171 int maxGeneration = simulatedFinal.back().particleId().generation();
173 for (
int gen = 0;
gen <= maxGeneration; ++
gen){
174 ATH_MSG_DEBUG(
name() <<
" start with generation "<< gen <<
"|" << maxGeneration <<
": "<< *itr);
175 auto vecsecisfp = std::make_unique<ISF::ISFParticleVector>();
176 while (itr != simulatedFinal.end() && static_cast<int>(itr->particleId().generation()) == gen) {
177 ATH_MSG_DEBUG(
name() <<
" genration "<< gen <<
"|" << maxGeneration <<
": "<< *itr);
178 if(itr->isSecondary()){
179
182 double mass = itr->mass() / Acts::UnitConstants::MeV;
183 double charge = itr->charge();
184 int pdgid = itr->pdg();
188 auto secisfp =
new ISF::ISFParticle (pos,mom,mass,
charge,pdgid,status,properTime,*isfp,
id);
190 ATH_MSG_DEBUG(
name() <<
" secondaries particle (ISF): " << *secisfp <<
" time "<<secisfp->timeStamp());
191 vecsecisfp->push_back(secisfp);
192 }
193 else{
195 }
196 ++itr;
197 }
198 if (!vecsecisfp->empty()) {
199 ISF::ISFTruthIncident truth(*isfp,
200 *vecsecisfp,
202 isfp->nextGeoID(),
204 );
205 ATH_MSG_DEBUG(
name() <<
" Truth incident parentPt2(MinPt2) " << truth.parentPt2() <<
" (100 MeV)");
206 ATH_MSG_DEBUG(
name() <<
" Truth incident ChildPt2(MinPt2) " << truth.childrenPt2Pass(300) <<
" (300 MeV)");
208 truth.updateParentAfterIncidentProperties();
209 truth.updateChildParticleProperties();
210 for (auto *secisfp : *vecsecisfp){
211 if (secisfp->getTruthBinding()) {
212 secondaries.push_back(secisfp);
213 }
214 else {
215 ATH_MSG_WARNING(
"Secondary particle not written out to truth.\n Parent (" << isfp <<
")\n Secondary (" << *secisfp <<
")");
216 }
217 }
218 }
219 }
220 }
223
224 std::vector<ActsFatras::Particle>().swap(input);
225 std::vector<ActsFatras::Particle>().swap(simulatedInitial);
226 std::vector<ActsFatras::Particle>().swap(simulatedFinal);
227 std::vector<ActsFatras::Hit>().swap(hits);
228 }
229 return StatusCode::SUCCESS;
230}
#define ATH_MSG_WARNING(x)
double charge(const T &p)
Acts::GeometryContext context() const
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.
constexpr int UNDEFINED_ID
constexpr int SIM_STATUS_THRESHOLD
Constant definiting the status threshold for simulated particles, eg. can be used to separate generat...
generator
Configure Herwig7 These are the commands corresponding to what would go into the regular Herwig infil...