File indexing completed on 2023-03-17 11:00:55
0001 #include "FastSimulation/SimplifiedGeometryPropagator/interface/ParticleManager.h"
0002
0003 #include "HepMC/GenEvent.h"
0004 #include "HepMC/Units.h"
0005 #include "HepPDT/ParticleDataTable.hh"
0006
0007 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
0008 #include "FastSimulation/SimplifiedGeometryPropagator/interface/ParticleFilter.h"
0009 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
0010 #include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
0011
0012 #include "SimDataFormats/Track/interface/SimTrack.h"
0013 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0016
0017 fastsim::ParticleManager::ParticleManager(const HepMC::GenEvent& genEvent,
0018 const HepPDT::ParticleDataTable& particleDataTable,
0019 double beamPipeRadius,
0020 double deltaRchargedMother,
0021 const fastsim::ParticleFilter& particleFilter,
0022 std::vector<SimTrack>& simTracks,
0023 std::vector<SimVertex>& simVertices,
0024 bool useFastSimsDecayer)
0025 : genEvent_(&genEvent),
0026 genParticleIterator_(genEvent_->particles_begin()),
0027 genParticleEnd_(genEvent_->particles_end()),
0028 genParticleIndex_(1),
0029 particleDataTable_(&particleDataTable),
0030 beamPipeRadius2_(beamPipeRadius * beamPipeRadius),
0031 deltaRchargedMother_(deltaRchargedMother),
0032 particleFilter_(&particleFilter),
0033 simTracks_(&simTracks),
0034 simVertices_(&simVertices),
0035 useFastSimsDecayer_(useFastSimsDecayer)
0036
0037
0038
0039
0040
0041
0042
0043
0044 ,
0045 momentumUnitConversionFactor_(conversion_factor(genEvent_->momentum_unit(), HepMC::Units::GEV)),
0046 lengthUnitConversionFactor_(conversion_factor(genEvent_->length_unit(), HepMC::Units::LengthUnit::CM)),
0047 lengthUnitConversionFactor2_(lengthUnitConversionFactor_ * lengthUnitConversionFactor_),
0048 timeUnitConversionFactor_(lengthUnitConversionFactor_ / fastsim::Constants::speedOfLight)
0049
0050 {
0051
0052 if (genEvent.vertices_begin() != genEvent_->vertices_end()) {
0053 const HepMC::FourVector& position = (*genEvent.vertices_begin())->position();
0054 addSimVertex(math::XYZTLorentzVector(position.x() * lengthUnitConversionFactor_,
0055 position.y() * lengthUnitConversionFactor_,
0056 position.z() * lengthUnitConversionFactor_,
0057 position.t() * timeUnitConversionFactor_),
0058 -1);
0059 }
0060 }
0061
0062 fastsim::ParticleManager::~ParticleManager() {}
0063
0064 std::unique_ptr<fastsim::Particle> fastsim::ParticleManager::nextParticle(const RandomEngineAndDistribution& random) {
0065 std::unique_ptr<fastsim::Particle> particle;
0066
0067
0068 if (!particleBuffer_.empty()) {
0069 particle = std::move(particleBuffer_.back());
0070 particleBuffer_.pop_back();
0071 }
0072
0073 else {
0074 particle = nextGenParticle();
0075 if (!particle)
0076 return nullptr;
0077 }
0078
0079
0080 if (!particleFilter_->accepts(*particle)) {
0081 return nextParticle(random);
0082 }
0083
0084
0085 if (!particle->remainingProperLifeTimeIsSet() || !particle->chargeIsSet()) {
0086
0087 const HepPDT::ParticleData* particleData = particleDataTable_->particle(HepPDT::ParticleID(particle->pdgId()));
0088 if (!particleData) {
0089
0090
0091 particle->setRemainingProperLifeTimeC(0.);
0092 particle->setCharge(0.);
0093 }
0094
0095
0096 if (!particle->remainingProperLifeTimeIsSet()) {
0097
0098
0099 double width = particleData->totalWidth().value();
0100 if (width > 1.0e-35) {
0101 particle->setRemainingProperLifeTimeC(-log(random.flatShoot()) * 6.582119e-25 / width / 10.);
0102 } else {
0103 particle->setStable();
0104 }
0105 }
0106
0107
0108 if (!particle->chargeIsSet()) {
0109 particle->setCharge(particleData->charge());
0110 }
0111 }
0112
0113
0114 unsigned simTrackIndex = addSimTrack(particle.get());
0115 particle->setSimTrackIndex(simTrackIndex);
0116
0117
0118 return particle;
0119 }
0120
0121 void fastsim::ParticleManager::addSecondaries(const math::XYZTLorentzVector& vertexPosition,
0122 int parentSimTrackIndex,
0123 std::vector<std::unique_ptr<Particle> >& secondaries,
0124 const SimplifiedGeometry* layer) {
0125
0126 if (!particleFilter_->acceptsVtx(vertexPosition)) {
0127 return;
0128 }
0129
0130
0131 if (secondaries.empty()) {
0132 return;
0133 }
0134
0135
0136 unsigned simVertexIndex = addSimVertex(vertexPosition, parentSimTrackIndex);
0137
0138
0139
0140 double distMin = 99999.;
0141 int idx = -1;
0142 int idxMin = -1;
0143 for (auto& secondary : secondaries) {
0144 idx++;
0145 if (secondary->getMotherDeltaR() != -1) {
0146 if (secondary->getMotherDeltaR() > deltaRchargedMother_) {
0147
0148 secondary->resetMother();
0149 } else {
0150 if (secondary->getMotherDeltaR() < distMin) {
0151 distMin = secondary->getMotherDeltaR();
0152 idxMin = idx;
0153 }
0154 }
0155 }
0156 }
0157
0158
0159 idx = -1;
0160 for (auto& secondary : secondaries) {
0161 idx++;
0162 if (idxMin != -1) {
0163
0164 if (secondary->getMotherDeltaR() != -1 && idx != idxMin) {
0165 secondary->resetMother();
0166 }
0167 }
0168
0169
0170 secondary->setSimVertexIndex(simVertexIndex);
0171
0172 if (layer) {
0173 secondary->setOnLayer(layer->isForward(), layer->index());
0174 }
0175
0176 particleBuffer_.push_back(std::move(secondary));
0177 }
0178 }
0179
0180 unsigned fastsim::ParticleManager::addEndVertex(const fastsim::Particle* particle) {
0181 return this->addSimVertex(particle->position(), particle->simTrackIndex());
0182 }
0183
0184 unsigned fastsim::ParticleManager::addSimVertex(const math::XYZTLorentzVector& position, int parentSimTrackIndex) {
0185 int simVertexIndex = simVertices_->size();
0186 simVertices_->emplace_back(position.Vect(), position.T(), parentSimTrackIndex, simVertexIndex);
0187 return simVertexIndex;
0188 }
0189
0190 unsigned fastsim::ParticleManager::addSimTrack(const fastsim::Particle* particle) {
0191 int simTrackIndex = simTracks_->size();
0192 simTracks_->emplace_back(
0193 particle->pdgId(), particle->momentum(), particle->simVertexIndex(), particle->genParticleIndex());
0194 simTracks_->back().setTrackId(simTrackIndex);
0195 return simTrackIndex;
0196 }
0197
0198 std::unique_ptr<fastsim::Particle> fastsim::ParticleManager::nextGenParticle() {
0199
0200
0201
0202
0203
0204
0205
0206 for (; genParticleIterator_ != genParticleEnd_; ++genParticleIterator_, ++genParticleIndex_) {
0207
0208 const HepMC::GenParticle& particle = **genParticleIterator_;
0209 const HepMC::GenVertex* productionVertex = particle.production_vertex();
0210 const HepMC::GenVertex* endVertex = particle.end_vertex();
0211
0212
0213 if (!productionVertex) {
0214 continue;
0215 }
0216 if (std::abs(particle.pdg_id()) < 10 || std::abs(particle.pdg_id()) == 21) {
0217 continue;
0218 }
0219
0220 int exoticRelativeId = 0;
0221 const bool producedWithinBeamPipe =
0222 productionVertex->position().perp2() * lengthUnitConversionFactor2_ < beamPipeRadius2_;
0223 if (!producedWithinBeamPipe && useFastSimsDecayer_) {
0224 exoticRelativesChecker(productionVertex, exoticRelativeId, 0);
0225 if (!isExotic(exoticRelativeId)) {
0226 continue;
0227 }
0228 }
0229
0230
0231 const bool decayedWithinBeamPipe =
0232 endVertex && endVertex->position().perp2() * lengthUnitConversionFactor2_ < beamPipeRadius2_;
0233 if (decayedWithinBeamPipe) {
0234 continue;
0235 }
0236
0237
0238 if (producedWithinBeamPipe && !decayedWithinBeamPipe && useFastSimsDecayer_) {
0239 exoticRelativesChecker(productionVertex, exoticRelativeId, 0);
0240 }
0241
0242
0243 std::unique_ptr<Particle> newParticle(
0244 new Particle(particle.pdg_id(),
0245 math::XYZTLorentzVector(productionVertex->position().x() * lengthUnitConversionFactor_,
0246 productionVertex->position().y() * lengthUnitConversionFactor_,
0247 productionVertex->position().z() * lengthUnitConversionFactor_,
0248 productionVertex->position().t() * timeUnitConversionFactor_),
0249 math::XYZTLorentzVector(particle.momentum().x() * momentumUnitConversionFactor_,
0250 particle.momentum().y() * momentumUnitConversionFactor_,
0251 particle.momentum().z() * momentumUnitConversionFactor_,
0252 particle.momentum().e() * momentumUnitConversionFactor_)));
0253 newParticle->setGenParticleIndex(genParticleIndex_);
0254 if (isExotic(exoticRelativeId)) {
0255 newParticle->setMotherPdgId(exoticRelativeId);
0256 }
0257
0258
0259 if (endVertex) {
0260 double labFrameLifeTime =
0261 (endVertex->position().t() - productionVertex->position().t()) * timeUnitConversionFactor_;
0262 newParticle->setRemainingProperLifeTimeC(labFrameLifeTime / newParticle->gamma() *
0263 fastsim::Constants::speedOfLight);
0264 }
0265
0266
0267
0268 bool foundVtx = false;
0269 for (const auto& simVtx : *simVertices_) {
0270 if (std::abs(simVtx.position().x() - newParticle->position().x()) < 1E-3 &&
0271 std::abs(simVtx.position().y() - newParticle->position().y()) < 1E-3 &&
0272 std::abs(simVtx.position().z() - newParticle->position().z()) < 1E-3) {
0273 newParticle->setSimVertexIndex(simVtx.vertexId());
0274 foundVtx = true;
0275 break;
0276 }
0277 }
0278 if (!foundVtx)
0279 newParticle->setSimVertexIndex(addSimVertex(newParticle->position(), -1));
0280
0281
0282 ++genParticleIterator_;
0283 ++genParticleIndex_;
0284
0285 return newParticle;
0286 }
0287
0288 return std::unique_ptr<Particle>();
0289 }
0290
0291 void fastsim::ParticleManager::exoticRelativesChecker(const HepMC::GenVertex* originVertex,
0292 int& exoticRelativeId_,
0293 int ngendepth = 0) {
0294 if (ngendepth > 99 || exoticRelativeId_ == -1 || isExotic(std::abs(exoticRelativeId_)))
0295 return;
0296 ngendepth += 1;
0297 std::vector<HepMC::GenParticle*>::const_iterator relativesIterator_ = originVertex->particles_in_const_begin();
0298 std::vector<HepMC::GenParticle*>::const_iterator relativesIteratorEnd_ = originVertex->particles_in_const_end();
0299 for (; relativesIterator_ != relativesIteratorEnd_; ++relativesIterator_) {
0300 const HepMC::GenParticle& genRelative = **relativesIterator_;
0301 if (isExotic(std::abs(genRelative.pdg_id()))) {
0302 exoticRelativeId_ = genRelative.pdg_id();
0303 if (ngendepth == 100)
0304 exoticRelativeId_ = -1;
0305 return;
0306 }
0307 const HepMC::GenVertex* vertex_ = genRelative.production_vertex();
0308 if (!vertex_)
0309 return;
0310 exoticRelativesChecker(vertex_, exoticRelativeId_, ngendepth);
0311 }
0312 return;
0313 }