84 {
85
87 CLHEP::HepRandomEngine* randomEngine = m_randomEngine->getEngine(ctx);
91 <<
particles.size() <<
" particles for simulation.");
92
93
95 auto bField = std::make_shared<ATLASMagneticFieldWrapper>();
104
113
114 simulator.charged.interactions = ActsFatras::makeStandardChargedElectroMagneticInteractions(
m_interact_minPt * Acts::UnitConstants::MeV);
115
120
122 for (const auto isfp : particles) {
123
124
125
126
127
128
129 ATH_MSG_DEBUG(
name() <<
" Convert ISF::Particle(mass) " << isfp->id()<<
"|" << isfp<<
"(" << isfp->mass() <<
")");
130 std::vector<ActsFatras::Particle>
input = std::vector<ActsFatras::Particle>{
131 ActsFatras::Particle(ActsFatras::Barcode().withVertexPrimary(0).withParticle(isfp->id()), static_cast<Acts::PdgParticle>(isfp->pdgCode()),
132 isfp->charge(),isfp->mass() * Acts::UnitConstants::MeV)
133 .setDirection(Acts::makeDirectionFromPhiEta(isfp->momentum().phi(), isfp->momentum().eta()))
134 .setAbsoluteMomentum(isfp->momentum().mag() * Acts::UnitConstants::MeV)
136 ATH_MSG_DEBUG(
name() <<
" Propagating ActsFatras::Particle vertex|particle|generation|subparticle, " << input[0]);
137 std::vector<ActsFatras::Particle> simulatedInitial;
138 std::vector<ActsFatras::Particle> simulatedFinal;
139 std::vector<ActsFatras::Hit>
hits;
140
141 auto result=simulator.simulate(anygctx, mctx, generator, input, simulatedInitial, simulatedFinal, hits);
142 auto simulatedFailure=
result.value();
143 if (simulatedFailure.size()>0){
144 for (const auto& simfail : simulatedFailure){
145 auto errCode = Acts::make_error_code(Acts::PropagatorError(simfail.error.value()));
146 ATH_MSG_WARNING(
name() <<
" Particle id " <<simfail.particle.particleId()<<
": fail to be simulated during Propagation: " << errCode.message());
147 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()));
148 return StatusCode::SUCCESS;
149 }
150 }
151
155 for (const auto& hit : hits) {
158 if (i>5) break;
159 }
160 ATH_MSG_DEBUG(
name() <<
" No. of particles after ActsFatras simulator: " << simulatedFinal.size());
161 if (!simulatedFinal.empty()){
163 auto itr = simulatedFinal.begin();
164
165 std::vector<ActsFatras::Hit> particle_hits;
166 if (itr->numberOfHits() > 0) {
167 std::copy(
hits.begin(),
hits.begin()+itr->numberOfHits(), std::back_inserter(particle_hits));
169 }
170
171 auto isKilled = !itr->isAlive();
172 int maxGeneration = simulatedFinal.back().particleId().generation();
174 for (
int gen = 0;
gen <= maxGeneration; ++
gen){
175 ATH_MSG_DEBUG(
name() <<
" start with generation "<< gen <<
"|" << maxGeneration <<
": "<< *itr);
176 auto vecsecisfp = std::make_unique<ISF::ISFParticleVector>();
177 while (itr != simulatedFinal.end() && static_cast<int>(itr->particleId().generation()) == gen) {
178 ATH_MSG_DEBUG(
name() <<
" genration "<< gen <<
"|" << maxGeneration <<
": "<< *itr);
179 if(itr->isSecondary()){
180
183 double mass = itr->mass() / Acts::UnitConstants::MeV;
184 double charge = itr->charge();
185 int pdgid = itr->pdg();
189 auto secisfp =
new ISF::ISFParticle (pos,mom,mass,
charge,pdgid,status,properTime,*isfp,
id);
191 ATH_MSG_DEBUG(
name() <<
" secondaries particle (ISF): " << *secisfp <<
" time "<<secisfp->timeStamp());
192 vecsecisfp->push_back(secisfp);
193 }
194 else{
196 }
197 ++itr;
198 }
199 if (!vecsecisfp->empty()) {
200 ISF::ISFTruthIncident truth(*isfp,
201 *vecsecisfp,
203 isfp->nextGeoID(),
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)");
209 truth.updateParentAfterIncidentProperties();
210 truth.updateChildParticleProperties();
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 }
219 }
220 }
221 }
224
225 std::vector<ActsFatras::Particle>().swap(input);
226 std::vector<ActsFatras::Particle>().swap(simulatedInitial);
227 std::vector<ActsFatras::Particle>().swap(simulatedFinal);
228 std::vector<ActsFatras::Hit>().swap(hits);
229 }
230 return StatusCode::SUCCESS;
231}
#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...