SimplifiedGeometryPropagator/src/ParticleFilter.cc

0001 #include "FastSimulation/SimplifiedGeometryPropagator/interface/ParticleFilter.h"
0002 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "vdt/vdtMath.h"
0005
0006 fastsim::ParticleFilter::ParticleFilter(const edm::ParameterSet& cfg) {
0007   // Charged particles must have pt greater than chargedPtMin [GeV]
0008   double chargedPtMin = cfg.getParameter<double>("chargedPtMin");
0009   chargedPtMin2_ = chargedPtMin * chargedPtMin;
0010
0011   // (All) particles must have energy greater than EMin [GeV]
0012   EMin_ = cfg.getParameter<double>("EMin");
0013
0014   // Allow *ALL* protons with energy > protonEMin
0015   protonEMin_ = cfg.getParameter<double>("protonEMin");
0016
0017   // List of invisible particles if extension needed
0018   // Predefined: Neutrinos, Neutralino_1
0019   skipParticles_ = cfg.getParameter<std::vector<int>>("invisibleParticles");
0020
0021   // Particles must have abs(eta) < etaMax (if close enough to 0,0,0)
0022   double etaMax = cfg.getParameter<double>("etaMax");
0023   cos2ThetaMax_ = (vdt::fast_exp(2. * etaMax) - 1.) / (vdt::fast_exp(2. * etaMax) + 1.);
0024   cos2ThetaMax_ *= cos2ThetaMax_;
0025
0026   // Particles must have vertex inside the tracker
0027   vertexRMax2_ = 129.0 * 129.0;
0028   vertexZMax_ = 303.353;
0029 }
0030
0031 bool fastsim::ParticleFilter::accepts(const fastsim::Particle& particle) const {
0032   int absPdgId = abs(particle.pdgId());
0033
0034   // skip invisible particles
0035   if (absPdgId == 12 || absPdgId == 14 || absPdgId == 16 || absPdgId == 1000022) {
0036     return false;
0037   }
0038   // keep all high-energy protons
0039   else if (absPdgId == 2212 && particle.momentum().E() >= protonEMin_) {
0040     return true;
0041   }
0042
0043   // cut on eta if the origin vertex is close to the beam
0044   else if (particle.position().Perp2() < 25. &&
0045            particle.momentum().Pz() * particle.momentum().Pz() / particle.momentum().P2() > cos2ThetaMax_) {
0046     return false;
0047   }
0048
0049   // possible to extend list of invisible particles
0050   for (unsigned InvIdx = 0; InvIdx < skipParticles_.size(); InvIdx++) {
0051     if (absPdgId == abs(skipParticles_.at(InvIdx))) {
0052       return false;
0053     }
0054   }
0055
0056   // particles must have vertex in volume of tracker
0057   return acceptsVtx(particle.position()) && acceptsEn(particle);
0058   //return (acceptsVtx(particle.position()) || particle.momentum().Pz()*particle.momentum().Pz()/particle.momentum().P2() > (vdt::fast_exp(2.*3.0)-1.) / (vdt::fast_exp(2.*3.0)+1.)*(vdt::fast_exp(2.*3.0)-1.) / (vdt::fast_exp(2.*3.0)+1.)) && acceptsEn(particle);
0059 }
0060
0061 bool fastsim::ParticleFilter::acceptsEn(const fastsim::Particle& particle) const {
0062   int absPdgId = abs(particle.pdgId());
0063
0064   // keep all high-energy protons
0065   if (absPdgId == 2212 && particle.momentum().E() >= protonEMin_) {
0066     return true;
0067   }
0068
0069   // cut on the energy
0070   else if (particle.momentum().E() < EMin_) {
0071     return false;
0072   }
0073
0074   // cut on pt of charged particles
0075   else if (particle.charge() != 0 && particle.momentum().Perp2() < chargedPtMin2_) {
0076     return false;
0077   }
0078
0079   return true;
0080 }
0081
0082 bool fastsim::ParticleFilter::acceptsVtx(const math::XYZTLorentzVector& originVertex) const {
0083   // origin vertex is within the tracker volume
0084   return (originVertex.Perp2() < vertexRMax2_ && fabs(originVertex.Z()) < vertexZMax_);
0085 }