SimplifiedGeometryPropagator/src/ParticleManager.cc

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       // prepare unit convsersions
0037       //  --------------------------------------------
0038       // |          |      hepmc               |  cms |
0039       //  --------------------------------------------
0040       // | length   | genEvent_->length_unit   |  cm  |
0041       // | momentum | genEvent_->momentum_unit |  GeV |
0042       // | time     | length unit (t*c)        |  ns  |
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   // add the main vertex from the signal event to the simvertex collection
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   // retrieve particle from buffer
0068   if (!particleBuffer_.empty()) {
0069     particle = std::move(particleBuffer_.back());
0070     particleBuffer_.pop_back();
0071   }
0072   // or from genParticle list
0073   else {
0074     particle = nextGenParticle();
0075     if (!particle)
0076       return nullptr;
0077   }
0078
0079   // if filter does not accept, skip particle
0080   if (!particleFilter_->accepts(*particle)) {
0081     return nextParticle(random);
0082   }
0083
0084   // lifetime or charge of particle are not yet set
0085   if (!particle->remainingProperLifeTimeIsSet() || !particle->chargeIsSet()) {
0086     // retrieve the particle data
0087     const HepPDT::ParticleData* particleData = particleDataTable_->particle(HepPDT::ParticleID(particle->pdgId()));
0088     if (!particleData) {
0089       // in very few events the Decayer (pythia) produces high mass resonances that are for some reason not present in the table (even though they should technically be)
0090       // they have short lifetimes, so decay them right away (charge and lifetime cannot be taken from table)
0091       particle->setRemainingProperLifeTimeC(0.);
0092       particle->setCharge(0.);
0093     }
0094
0095     // set lifetime
0096     if (!particle->remainingProperLifeTimeIsSet()) {
0097       // The lifetime is 0. in the Pythia Particle Data Table! Calculate from width instead (ct=hbar/width).
0098       // ct=particleData->lifetime().value();
0099       double width = particleData->totalWidth().value();
0100       if (width > 1.0e-35) {
0101         particle->setRemainingProperLifeTimeC(-log(random.flatShoot()) * 6.582119e-25 / width / 10.);  // ct in cm
0102       } else {
0103         particle->setStable();
0104       }
0105     }
0106
0107     // set charge
0108     if (!particle->chargeIsSet()) {
0109       particle->setCharge(particleData->charge());
0110     }
0111   }
0112
0113   // add corresponding simTrack to simTrack collection
0114   unsigned simTrackIndex = addSimTrack(particle.get());
0115   particle->setSimTrackIndex(simTrackIndex);
0116
0117   // and return
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   // vertex must be within the accepted volume
0126   if (!particleFilter_->acceptsVtx(vertexPosition)) {
0127     return;
0128   }
0129
0130   // no need to create vertex in case no particles are produced
0131   if (secondaries.empty()) {
0132     return;
0133   }
0134
0135   // add simVertex
0136   unsigned simVertexIndex = addSimVertex(vertexPosition, parentSimTrackIndex);
0137
0138   // closest charged daughter continues the track of the mother particle
0139   // simplified tracking algorithm for fastSim
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         // larger than max requirement on deltaR
0148         secondary->resetMother();
0149       } else {
0150         if (secondary->getMotherDeltaR() < distMin) {
0151           distMin = secondary->getMotherDeltaR();
0152           idxMin = idx;
0153         }
0154       }
0155     }
0156   }
0157
0158   // add secondaries to buffer
0159   idx = -1;
0160   for (auto& secondary : secondaries) {
0161     idx++;
0162     if (idxMin != -1) {
0163       // reset all but the particle with the lowest deltaR (which is at idxMin)
0164       if (secondary->getMotherDeltaR() != -1 && idx != idxMin) {
0165         secondary->resetMother();
0166       }
0167     }
0168
0169     // set origin vertex
0170     secondary->setSimVertexIndex(simVertexIndex);
0171     //
0172     if (layer) {
0173       secondary->setOnLayer(layer->isForward(), layer->index());
0174     }
0175     // ...and add particle to buffer
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   // only consider particles that start in the beam pipe and end outside the beam pipe
0200   // try to get the decay time from pythia
0201   // use hepmc units
0202   // make the link simtrack to simvertex
0203   // try not to change the simvertex structure
0204
0205   // loop over gen particles
0206   for (; genParticleIterator_ != genParticleEnd_; ++genParticleIterator_, ++genParticleIndex_) {
0207     // some handy pointers and references
0208     const HepMC::GenParticle& particle = **genParticleIterator_;
0209     const HepMC::GenVertex* productionVertex = particle.production_vertex();
0210     const HepMC::GenVertex* endVertex = particle.end_vertex();
0211
0212     // skip incoming particles
0213     if (!productionVertex) {
0214       continue;
0215     }
0216     if (std::abs(particle.pdg_id()) < 10 || std::abs(particle.pdg_id()) == 21) {
0217       continue;
0218     }
0219     // particles which do not descend from exotics must be produced within the beampipe
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     // FastSim will not make hits out of particles that decay before reaching the beam pipe
0231     const bool decayedWithinBeamPipe =
0232         endVertex && endVertex->position().perp2() * lengthUnitConversionFactor2_ < beamPipeRadius2_;
0233     if (decayedWithinBeamPipe) {
0234       continue;
0235     }
0236
0237     // SM particles that descend from exotics and cross the beam pipe radius should make hits but not be decayed
0238     if (producedWithinBeamPipe && !decayedWithinBeamPipe && useFastSimsDecayer_) {
0239       exoticRelativesChecker(productionVertex, exoticRelativeId, 0);
0240     }
0241
0242     // make the particle
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     // try to get the life time of the particle from the genEvent
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     // Find production vertex if it already exists. Otherwise create new vertex
0267     // Possible to recreate the whole GenEvent using SimTracks/SimVertices (see FBaseSimEvent::fill(..))
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     // iterator/index has to be increased in case of return (is not done by the loop then)
0282     ++genParticleIterator_;
0283     ++genParticleIndex_;
0284     // and return
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 }