SimplifiedGeometryPropagator/plugins/TrackerSimHitProducer.cc

0001 #include <memory>
0002
0003 #include <memory>
0004 #include <vector>
0005
0006 // framework
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/ProducesCollector.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "MagneticField/UniformEngine/interface/UniformMagneticField.h"
0011
0012 // tracking
0013 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0014 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0015 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0016 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0017 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0018
0019 // fastsim
0020 #include "FastSimulation/TrajectoryManager/interface/InsideBoundsMeasurementEstimator.h"  //TODO move this
0021 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
0022 #include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
0023 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModel.h"
0024 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModelFactory.h"
0025 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
0026
0027 // data formats
0028 #include "DataFormats/GeometrySurface/interface/Plane.h"
0029 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0030 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0031 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0032 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0033 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0034 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0035
0036 // other
0037 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0038 #include "CondFormats/External/interface/DetID.h"
0039
0040 ///////////////////////////////////////////////
0041 // Author: L. Vanelderen, S. Kurz
0042 // Date: 29 May 2017
0043 //////////////////////////////////////////////////////////
0044
0045 typedef std::pair<const GeomDet*, TrajectoryStateOnSurface> DetWithState;
0046
0047 namespace fastsim {
0048   //! Produces SimHits in the tracker layers.
0049   /*!
0050         Also assigns the energy loss of the particle (ionization) with the SimHit.
0051         Furthermore, SimHits from different SubModules have to be sorted by their occurance!
0052         \sa EnergyLoss
0053     */
0054   class TrackerSimHitProducer : public InteractionModel {
0055   public:
0056     //! Constructor.
0057     TrackerSimHitProducer(const std::string& name, const edm::ParameterSet& cfg);
0058
0059     //! Default destructor.
0060     ~TrackerSimHitProducer() override { ; }
0061
0062     //! Perform the interaction.
0063     /*!
0064             \param particle The particle that interacts with the matter.
0065             \param layer The detector layer that interacts with the particle.
0066             \param secondaries Particles that are produced in the interaction (if any).
0067             \param random The Random Engine.
0068         */
0069     void interact(Particle& particle,
0070                   const SimplifiedGeometry& layer,
0071                   std::vector<std::unique_ptr<Particle>>& secondaries,
0072                   const RandomEngineAndDistribution& random) override;
0073
0074     //! Register the SimHit collection.
0075     void registerProducts(edm::ProducesCollector) const override;
0076
0077     //! Store the SimHit collection.
0078     void storeProducts(edm::Event& iEvent) override;
0079
0080     //! Helper funtion to create the actual SimHit on a detector (sub-) module.
0081     /*!
0082             \param particle Representation of the particle's trajectory
0083             \param pdgId The pdgId of the particle
0084             \param layerThickness The thickness of the layer (in radLengths)
0085             \param eLoss The energy that particle deposited in the detector (lost via ionisation)
0086             \param simTrackId The SimTrack this hit should be assigned to
0087             \param detector The detector element that is hit
0088             \param refPos Reference position that is used to sort the hits if layer consists of sub-detectors
0089             \ return Returns the SimHit and the distance relative to refPos since hits have to be ordered (in time) afterwards.
0090         */
0091     std::pair<double, std::unique_ptr<PSimHit>> createHitOnDetector(const TrajectoryStateOnSurface& particle,
0092                                                                     int pdgId,
0093                                                                     double layerThickness,
0094                                                                     double eLoss,
0095                                                                     int simTrackId,
0096                                                                     const GeomDet& detector,
0097                                                                     GlobalPoint& refPos);
0098
0099   private:
0100     const double
0101         onSurfaceTolerance_;  //!< Max distance between particle and active (sub-) module. Otherwise particle has to be propagated.
0102     std::unique_ptr<edm::PSimHitContainer> simHitContainer_;  //!< The SimHit.
0103     double minMomentum_;                                      //!< Set the minimal momentum of incoming particle
0104     bool doHitsFromInboundParticles_;  //!< If not set, incoming particles (negative speed relative to center of detector) don't create a SimHits since reconstruction anyways not possible
0105   };
0106 }  // namespace fastsim
0107
0108 fastsim::TrackerSimHitProducer::TrackerSimHitProducer(const std::string& name, const edm::ParameterSet& cfg)
0109     : fastsim::InteractionModel(name), onSurfaceTolerance_(0.01), simHitContainer_(new edm::PSimHitContainer) {
0110   // Set the minimal momentum
0111   minMomentum_ = cfg.getParameter<double>("minMomentumCut");
0112   // - if not set, particles from outside the beampipe with a negative speed in R direction are propagated but no SimHits
0113   // - particle with positive (negative) z and negative (positive) speed in z direction: no SimHits
0114   // -> this is not neccesary since a track reconstruction is not possible in this case anyways
0115   doHitsFromInboundParticles_ = cfg.getParameter<bool>("doHitsFromInboundParticles");
0116 }
0117
0118 void fastsim::TrackerSimHitProducer::registerProducts(edm::ProducesCollector producesCollector) const {
0119   producesCollector.produces<edm::PSimHitContainer>("TrackerHits");
0120 }
0121
0122 void fastsim::TrackerSimHitProducer::storeProducts(edm::Event& iEvent) {
0123   iEvent.put(std::move(simHitContainer_), "TrackerHits");
0124   simHitContainer_ = std::make_unique<edm::PSimHitContainer>();
0125 }
0126
0127 void fastsim::TrackerSimHitProducer::interact(Particle& particle,
0128                                               const SimplifiedGeometry& layer,
0129                                               std::vector<std::unique_ptr<Particle>>& secondaries,
0130                                               const RandomEngineAndDistribution& random) {
0131   // the energy deposit in the layer
0132   double energyDeposit = particle.getEnergyDeposit();
0133   particle.setEnergyDeposit(0);
0134
0135   //
0136   // don't save hits from particle that did a loop or are inbound (coming from the outer part of the tracker, going towards the center)
0137   //
0138   if (!doHitsFromInboundParticles_) {
0139     if (particle.isLooper()) {
0140       return;
0141     }
0142     if (particle.momentum().X() * particle.position().X() + particle.momentum().Y() * particle.position().Y() < 0) {
0143       particle.setLooper();
0144       return;
0145     }
0146     if (layer.isForward() && ((particle.position().Z() > 0 && particle.momentum().Z() < 0) ||
0147                               (particle.position().Z() < 0 && particle.momentum().Z() > 0))) {
0148       particle.setLooper();
0149       return;
0150     }
0151   }
0152
0153   //
0154   // check that layer has tracker modules
0155   //
0156   if (!layer.getDetLayer()) {
0157     return;
0158   }
0159
0160   //
0161   // no material
0162   //
0163   if (layer.getThickness(particle.position(), particle.momentum()) < 1E-10) {
0164     return;
0165   }
0166
0167   //
0168   // only charged particles
0169   //
0170   if (particle.charge() == 0) {
0171     return;
0172   }
0173
0174   //
0175   // save hit only if momentum higher than threshold
0176   //
0177   if (particle.momentum().Perp2() < minMomentum_ * minMomentum_) {
0178     return;
0179   }
0180
0181   // create the trajectory of the particle
0182   UniformMagneticField magneticField(layer.getMagneticFieldZ(particle.position()));
0183   GlobalPoint position(particle.position().X(), particle.position().Y(), particle.position().Z());
0184   GlobalVector momentum(particle.momentum().Px(), particle.momentum().Py(), particle.momentum().Pz());
0185   auto plane = layer.getDetLayer()->surface().tangentPlane(position);
0186   TrajectoryStateOnSurface trajectory(
0187       GlobalTrajectoryParameters(position, momentum, TrackCharge(particle.charge()), &magneticField), *plane);
0188
0189   // find detectors compatible with the particle's trajectory
0190   AnalyticalPropagator propagator(&magneticField, anyDirection);
0191   InsideBoundsMeasurementEstimator est;
0192   std::vector<DetWithState> compatibleDetectors = layer.getDetLayer()->compatibleDets(trajectory, propagator, est);
0193
0194   ////////
0195   // You have to sort the simHits in the order they occur!
0196   ////////
0197
0198   // The old algorithm (sorting by distance to IP) doesn't seem to make sense to me (what if particle moves inwards?)
0199
0200   // Detector layers have to be sorted by proximity to particle.position
0201   // Propagate particle backwards a bit to make sure it's outside any components
0202   std::map<double, std::unique_ptr<PSimHit>> distAndHits;
0203   // Position relative to which the hits should be sorted
0204   GlobalPoint positionOutside(particle.position().x() - particle.momentum().x() / particle.momentum().P() * 10.,
0205                               particle.position().y() - particle.momentum().y() / particle.momentum().P() * 10.,
0206                               particle.position().z() - particle.momentum().z() / particle.momentum().P() * 10.);
0207
0208   // FastSim: cheat tracking -> assign hits to closest charged daughter if particle decays
0209   int pdgId = particle.pdgId();
0210   int simTrackId =
0211       particle.getMotherSimTrackIndex() >= 0 ? particle.getMotherSimTrackIndex() : particle.simTrackIndex();
0212
0213   // loop over the compatible detectors
0214   for (const auto& detectorWithState : compatibleDetectors) {
0215     const GeomDet& detector = *detectorWithState.first;
0216     const TrajectoryStateOnSurface& particleState = detectorWithState.second;
0217
0218     // if the detector has no components
0219     if (detector.isLeaf()) {
0220       std::pair<double, std::unique_ptr<PSimHit>> hitPair = createHitOnDetector(particleState,
0221                                                                                 pdgId,
0222                                                                                 layer.getThickness(particle.position()),
0223                                                                                 energyDeposit,
0224                                                                                 simTrackId,
0225                                                                                 detector,
0226                                                                                 positionOutside);
0227       if (hitPair.second) {
0228         distAndHits.insert(distAndHits.end(), std::move(hitPair));
0229       }
0230     } else {
0231       // if the detector has components
0232       for (const auto component : detector.components()) {
0233         std::pair<double, std::unique_ptr<PSimHit>> hitPair =
0234             createHitOnDetector(particleState,
0235                                 pdgId,
0236                                 layer.getThickness(particle.position()),
0237                                 energyDeposit,
0238                                 simTrackId,
0239                                 *component,
0240                                 positionOutside);
0241         if (hitPair.second) {
0242           distAndHits.insert(distAndHits.end(), std::move(hitPair));
0243         }
0244       }
0245     }
0246   }
0247
0248   // Fill simHitContainer
0249   for (std::map<double, std::unique_ptr<PSimHit>>::const_iterator it = distAndHits.begin(); it != distAndHits.end();
0250        it++) {
0251     simHitContainer_->push_back(*(it->second));
0252   }
0253 }
0254
0255 // Also returns distance to simHit since hits have to be ordered (in time) afterwards. Necessary to do explicit copy of TrajectoryStateOnSurface particle (not call by reference)
0256 std::pair<double, std::unique_ptr<PSimHit>> fastsim::TrackerSimHitProducer::createHitOnDetector(
0257     const TrajectoryStateOnSurface& particle,
0258     int pdgId,
0259     double layerThickness,
0260     double eLoss,
0261     int simTrackId,
0262     const GeomDet& detector,
0263     GlobalPoint& refPos) {
0264   // determine position and momentum of particle in the coordinate system of the detector
0265   LocalPoint localPosition;
0266   LocalVector localMomentum;
0267
0268   // if the particle is close enough, no further propagation is needed
0269   if (fabs(detector.toLocal(particle.globalPosition()).z()) < onSurfaceTolerance_) {
0270     localPosition = particle.localPosition();
0271     localMomentum = particle.localMomentum();
0272   }
0273   // else, propagate
0274   else {
0275     // find crossing of particle with detector layer
0276     HelixArbitraryPlaneCrossing crossing(particle.globalPosition().basicVector(),
0277                                          particle.globalMomentum().basicVector(),
0278                                          particle.transverseCurvature(),
0279                                          anyDirection);
0280     std::pair<bool, double> path = crossing.pathLength(detector.surface());
0281     // case propagation succeeds
0282     if (path.first) {
0283       localPosition = detector.toLocal(GlobalPoint(crossing.position(path.second)));
0284       localMomentum = detector.toLocal(GlobalVector(crossing.direction(path.second)));
0285       localMomentum = localMomentum.unit() * particle.localMomentum().mag();
0286     }
0287     // case propagation fails
0288     else {
0289       return std::pair<double, std::unique_ptr<PSimHit>>(0, std::unique_ptr<PSimHit>(nullptr));
0290     }
0291   }
0292
0293   // find entry and exit point of particle in detector
0294   const Plane& detectorPlane = detector.surface();
0295   double halfThick = 0.5 * detectorPlane.bounds().thickness();
0296   double pZ = localMomentum.z();
0297   LocalPoint entry = localPosition + (-halfThick / pZ) * localMomentum;
0298   LocalPoint exit = localPosition + (halfThick / pZ) * localMomentum;
0299   double tof = particle.globalPosition().mag() / fastsim::Constants::speedOfLight;  // in nanoseconds
0300
0301   // make sure the simhit is physically on the module
0302   double boundX = detectorPlane.bounds().width() / 2.;
0303   double boundY = detectorPlane.bounds().length() / 2.;
0304   // Special treatment for TID and TEC trapeziodal modules
0305   unsigned subdet = DetId(detector.geographicalId()).subdetId();
0306   if (subdet == 4 || subdet == 6)
0307     boundX *= 1. - localPosition.y() / detectorPlane.position().perp();
0308   if (fabs(localPosition.x()) > boundX || fabs(localPosition.y()) > boundY) {
0309     return std::pair<double, std::unique_ptr<PSimHit>>(0, std::unique_ptr<PSimHit>(nullptr));
0310   }
0311
0312   // Create the hit: the energy loss rescaled to the module thickness
0313   // Total thickness is in radiation lengths, 1 radlen = 9.36 cm
0314   // Sensitive module thickness is about 30 microns larger than
0315   // the module thickness itself
0316   eLoss *= (2. * halfThick - 0.003) / (9.36 * layerThickness);
0317
0318   // Position of the hit in global coordinates
0319   GlobalPoint hitPos(detector.surface().toGlobal(localPosition));
0320
0321   return std::pair<double, std::unique_ptr<PSimHit>>((hitPos - refPos).mag(),
0322                                                      std::make_unique<PSimHit>(entry,
0323                                                                                exit,
0324                                                                                localMomentum.mag(),
0325                                                                                tof,
0326                                                                                eLoss,
0327                                                                                pdgId,
0328                                                                                detector.geographicalId().rawId(),
0329                                                                                simTrackId,
0330                                                                                localMomentum.theta(),
0331                                                                                localMomentum.phi()));
0332 }
0333
0334 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::TrackerSimHitProducer, "fastsim::TrackerSimHitProducer");