ParticlePropagator/src/ParticlePropagator.cc

0001 //CMSSW headers
0002 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0003 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
0004
0005 //FAMOS headers
0006 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
0007 #include "FastSimulation/ParticlePropagator/interface/MagneticFieldMap.h"
0008 #include "FastSimulation/TrackerSetup/interface/TrackerLayer.h"
0009 #include "FastSimulation/Event/interface/FSimTrack.h"
0010 #include "FastSimulation/Event/interface/FSimVertex.h"
0011 #include "FastSimulation/Particle/interface/pdg_functions.h"
0012 #include "FastSimulation/Particle/interface/makeParticle.h"
0013 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0014
0015 ParticlePropagator::ParticlePropagator() : BaseParticlePropagator(), random(nullptr) { ; }
0016
0017 ParticlePropagator::ParticlePropagator(const RawParticle& myPart,
0018                                        double RCyl,
0019                                        double ZCyl,
0020                                        const MagneticFieldMap* aFieldMap,
0021                                        const RandomEngineAndDistribution* engine,
0022                                        const HepPDT::ParticleDataTable* table)
0023     : BaseParticlePropagator(myPart, RCyl, ZCyl, 0.), theFieldMap(aFieldMap), random(engine), theTable(table) {
0024   setMagneticField(fieldMap(particle().X(), particle().Y(), particle().Z()));
0025   initProperDecayTime();
0026 }
0027
0028 ParticlePropagator::ParticlePropagator(const RawParticle& myPart,
0029                                        const MagneticFieldMap* aFieldMap,
0030                                        const RandomEngineAndDistribution* engine,
0031                                        const HepPDT::ParticleDataTable* table)
0032     : BaseParticlePropagator(myPart, 0., 0., 0.),
0033       theFieldMap(aFieldMap),
0034       random(engine),
0035       theTable(table)
0036
0037 {
0038   setMagneticField(fieldMap(particle().X(), particle().Y(), particle().Z()));
0039   initProperDecayTime();
0040 }
0041
0042 ParticlePropagator::ParticlePropagator(const XYZTLorentzVector& mom,
0043                                        const XYZTLorentzVector& vert,
0044                                        float q,
0045                                        const MagneticFieldMap* aFieldMap,
0046                                        const HepPDT::ParticleDataTable* table)
0047     : BaseParticlePropagator(RawParticle(mom, vert, q), 0., 0., 0.),
0048       theFieldMap(aFieldMap),
0049       random(nullptr),
0050       theTable(table) {
0051   setMagneticField(fieldMap(particle().X(), particle().Y(), particle().Z()));
0052 }
0053
0054 ParticlePropagator::ParticlePropagator(const XYZTLorentzVector& mom,
0055                                        const XYZVector& vert,
0056                                        float q,
0057                                        const MagneticFieldMap* aFieldMap,
0058                                        const HepPDT::ParticleDataTable* table)
0059     : BaseParticlePropagator(RawParticle(mom, XYZTLorentzVector(vert.X(), vert.Y(), vert.Z(), 0.0), q), 0., 0., 0.),
0060       theFieldMap(aFieldMap),
0061       random(nullptr),
0062       theTable(table) {
0063   setMagneticField(fieldMap(particle().X(), particle().Y(), particle().Z()));
0064 }
0065
0066 ParticlePropagator::ParticlePropagator(const FSimTrack& simTrack,
0067                                        const MagneticFieldMap* aFieldMap,
0068                                        const RandomEngineAndDistribution* engine,
0069                                        const HepPDT::ParticleDataTable* table)
0070     : BaseParticlePropagator(
0071           makeParticle(table, simTrack.type(), simTrack.momentum(), simTrack.vertex().position()), 0., 0., 0.),
0072       theFieldMap(aFieldMap),
0073       random(engine),
0074       theTable(table) {
0075   setMagneticField(fieldMap(particle().X(), particle().Y(), particle().Z()));
0076   if (simTrack.decayTime() < 0.) {
0077     if (simTrack.nDaughters())
0078       // This particle already decayed, don't decay it twice
0079       this->setProperDecayTime(1E99);
0080     else
0081       // This particle hasn't decayed yet. Decay time according to particle lifetime
0082       initProperDecayTime();
0083   } else {
0084     // Decay time pre-defined at generator level
0085     this->setProperDecayTime(simTrack.decayTime());
0086   }
0087 }
0088
0089 ParticlePropagator::ParticlePropagator(const ParticlePropagator& myPropPart)
0090     : BaseParticlePropagator(myPropPart), theFieldMap(myPropPart.theFieldMap) {
0091   //  setMagneticField(fieldMap(x(),y(),z()));
0092 }
0093
0094 ParticlePropagator::ParticlePropagator(const BaseParticlePropagator& myPropPart,
0095                                        const MagneticFieldMap* aFieldMap,
0096                                        const HepPDT::ParticleDataTable* table)
0097     : BaseParticlePropagator(myPropPart), theFieldMap(aFieldMap), theTable(table) {
0098   setMagneticField(fieldMap(particle().X(), particle().Y(), particle().Z()));
0099 }
0100
0101 void ParticlePropagator::initProperDecayTime() {
0102   // And this is the proper time at which the particle will decay
0103   double properDecayTime = (particle().pid() == 0 || particle().pid() == 22 || abs(particle().pid()) == 11 ||
0104                             abs(particle().pid()) == 2112 || abs(particle().pid()) == 2212 || !random)
0105                                ? 1E99
0106                                : -pdg::cTau(particle().pid(), theTable) * std::log(random->flatShoot());
0107
0108   this->setProperDecayTime(properDecayTime);
0109 }
0110
0111 bool ParticlePropagator::propagateToClosestApproach(double x0, double y0, bool first) {
0112   setMagneticField(fieldMap(0., 0., 0.));
0113   return BaseParticlePropagator::propagateToClosestApproach(x0, y0, first);
0114 }
0115
0116 bool ParticlePropagator::propagateToNominalVertex(const XYZTLorentzVector& v) {
0117   setMagneticField(fieldMap(0., 0., 0.));
0118   return BaseParticlePropagator::propagateToNominalVertex(v);
0119 }
0120
0121 ParticlePropagator ParticlePropagator::propagated() const {
0122   return ParticlePropagator(BaseParticlePropagator::propagated(), theFieldMap, theTable);
0123 }
0124
0125 double ParticlePropagator::fieldMap(double xx, double yy, double zz) {
0126   // Arguments now passed in cm.
0127   //  return MagneticFieldMap::instance()->inTesla(GlobalPoint(xx/10.,yy/10.,zz/10.)).z();
0128   // Return a dummy value for neutral particles!
0129   return particle().charge() == 0.0 || theFieldMap == nullptr ? 4. : theFieldMap->inTeslaZ(GlobalPoint(xx, yy, zz));
0130 }
0131
0132 double ParticlePropagator::fieldMap(const TrackerLayer& layer, double coord, int success) {
0133   // Arguments now passed in cm.
0134   //  return MagneticFieldMap::instance()->inTesla(GlobalPoint(xx/10.,yy/10.,zz/10.)).z();
0135   // Return a dummy value for neutral particles!
0136   return particle().charge() == 0.0 || theFieldMap == nullptr ? 4. : theFieldMap->inTeslaZ(layer, coord, success);
0137 }
0138
0139 bool ParticlePropagator::propagateToBoundSurface(const TrackerLayer& layer) {
0140   fiducial = true;
0141   BoundDisk const* disk = layer.disk();
0142   //  bool disk = layer.forward();
0143   //  double innerradius=-999;
0144   double innerradius = disk ? layer.diskInnerRadius() : -999.;
0145
0146   //  if( disk ) {
0147   //    const Surface& surface = layer.surface();
0148   //    const BoundDisk & myDisk = dynamic_cast<const BoundDisk&>(surface);
0149   //    innerradius=myDisk.innerRadius();
0150   //    innerradius=myDisk->innerRadius();
0151   //  }
0152
0153   bool done = propagate();
0154
0155   // Set the magnetic field at the new location (if succesfully propagated)
0156   if (done && !hasDecayed()) {
0157     if (success == 2)
0158       setMagneticField(fieldMap(layer, particle().r(), success));
0159     else if (success == 1)
0160       setMagneticField(fieldMap(layer, particle().z(), success));
0161   }
0162
0163   // There is some real material here
0164   fiducial = !(!disk && success != 1) && !(disk && (success != 2 || particle().r() < innerradius));
0165
0166   return done;
0167 }
0168
0169 void ParticlePropagator::setPropagationConditions(const TrackerLayer& layer, bool firstLoop) {
0170   // Set the magentic field
0171   // setMagneticField(fieldMap(x(),y(),z()));
0172
0173   // Set R and Z according to the Tracker Layer characteristics.
0174   //  const Surface& surface = layer.surface();
0175
0176   if (layer.forward()) {
0177     //    const BoundDisk & myDisk = dynamic_cast<const BoundDisk&>(surface);
0178     // ParticlePropagator works in mm, whereas the detector geometry is in cm
0179     BaseParticlePropagator::setPropagationConditions(
0180         layer.diskOuterRadius(), fabs(layer.disk()->position().z()), firstLoop);
0181
0182     // ... or if it is a cylinder barrel
0183   } else {
0184     //    const BoundCylinder & myCylinder = dynamic_cast<const BoundCylinder &>(surface);
0185     // ParticlePropagator works now in cm
0186     BaseParticlePropagator::setPropagationConditions(
0187         layer.cylinder()->bounds().width() / 2., layer.cylinder()->bounds().length() / 2., firstLoop);
0188   }
0189 }