TrackerSimHitProducer.cc

Go to the documentation of this file.

#include <memory>

#include <memory>
#include <vector>

// framework
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/ProducesCollector.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "MagneticField/UniformEngine/interface/UniformMagneticField.h"

// tracking
#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"

// fastsim
#include "FastSimulation/TrajectoryManager/interface/InsideBoundsMeasurementEstimator.h"  //TODO move this
#include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModel.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModelFactory.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"

// data formats
#include "DataFormats/GeometrySurface/interface/Plane.h"
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/LocalVector.h"
#include "DataFormats/GeometryVector/interface/LocalPoint.h"
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"

// other
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "CondFormats/External/interface/DetID.h"

// Author: L. Vanelderen, S. Kurz
// Date: 29 May 2017

typedef std::pair<const GeomDet*, TrajectoryStateOnSurface> DetWithState;

namespace fastsim {

  class TrackerSimHitProducer : public InteractionModel {
  public:
    TrackerSimHitProducer(const std::string& name, const edm::ParameterSet& cfg);

    ~TrackerSimHitProducer() override { ; }


    void interact(Particle& particle,
                  const SimplifiedGeometry& layer,
                  std::vector<std::unique_ptr<Particle>>& secondaries,
                  const RandomEngineAndDistribution& random) override;

    void registerProducts(edm::ProducesCollector) const override;

    void storeProducts(edm::Event& iEvent) override;


    std::pair<double, std::unique_ptr<PSimHit>> createHitOnDetector(const TrajectoryStateOnSurface& particle,
                                                                    int pdgId,
                                                                    double layerThickness,
                                                                    double eLoss,
                                                                    int simTrackId,
                                                                    const GeomDet& detector,
                                                                    GlobalPoint& refPos);

  private:
    const double
        onSurfaceTolerance_;  
    std::unique_ptr<edm::PSimHitContainer> simHitContainer_;  
    double minMomentum_;                                      
    bool doHitsFromInboundParticles_;  
  };
}  // namespace fastsim

fastsim::TrackerSimHitProducer::TrackerSimHitProducer(const std::string& name, const edm::ParameterSet& cfg)
    : fastsim::InteractionModel(name), onSurfaceTolerance_(0.01), simHitContainer_(new edm::PSimHitContainer) {
  // Set the minimal momentum
  minMomentum_ = cfg.getParameter<double>("minMomentumCut");
  // - if not set, particles from outside the beampipe with a negative speed in R direction are propagated but no SimHits
  // - particle with positive (negative) z and negative (positive) speed in z direction: no SimHits
  // -> this is not neccesary since a track reconstruction is not possible in this case anyways
  doHitsFromInboundParticles_ = cfg.getParameter<bool>("doHitsFromInboundParticles");
}

void fastsim::TrackerSimHitProducer::registerProducts(edm::ProducesCollector producesCollector) const {
  producesCollector.produces<edm::PSimHitContainer>("TrackerHits");
}

void fastsim::TrackerSimHitProducer::storeProducts(edm::Event& iEvent) {
  iEvent.put(std::move(simHitContainer_), "TrackerHits");
  simHitContainer_ = std::make_unique<edm::PSimHitContainer>();
}

void fastsim::TrackerSimHitProducer::interact(Particle& particle,
                                              const SimplifiedGeometry& layer,
                                              std::vector<std::unique_ptr<Particle>>& secondaries,
                                              const RandomEngineAndDistribution& random) {
  // the energy deposit in the layer
  double energyDeposit = particle.getEnergyDeposit();
  particle.setEnergyDeposit(0);

  //
  // 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)
  //
  if (!doHitsFromInboundParticles_) {
    if (particle.isLooper()) {
      return;
    }
    if (particle.momentum().X() * particle.position().X() + particle.momentum().Y() * particle.position().Y() < 0) {
      particle.setLooper();
      return;
    }
    if (layer.isForward() && ((particle.position().Z() > 0 && particle.momentum().Z() < 0) ||
                              (particle.position().Z() < 0 && particle.momentum().Z() > 0))) {
      particle.setLooper();
      return;
    }
  }

  //
  // check that layer has tracker modules
  //
  if (!layer.getDetLayer()) {
    return;
  }

  //
  // no material
  //
  if (layer.getThickness(particle.position(), particle.momentum()) < 1E-10) {
    return;
  }

  //
  // only charged particles
  //
  if (particle.charge() == 0) {
    return;
  }

  //
  // save hit only if momentum higher than threshold
  //
  if (particle.momentum().Perp2() < minMomentum_ * minMomentum_) {
    return;
  }

  // create the trajectory of the particle
  UniformMagneticField magneticField(layer.getMagneticFieldZ(particle.position()));
  GlobalPoint position(particle.position().X(), particle.position().Y(), particle.position().Z());
  GlobalVector momentum(particle.momentum().Px(), particle.momentum().Py(), particle.momentum().Pz());
  auto plane = layer.getDetLayer()->surface().tangentPlane(position);
  TrajectoryStateOnSurface trajectory(
      GlobalTrajectoryParameters(position, momentum, TrackCharge(particle.charge()), &magneticField), *plane);

  // find detectors compatible with the particle's trajectory
  AnalyticalPropagator propagator(&magneticField, anyDirection);
  InsideBoundsMeasurementEstimator est;
  std::vector<DetWithState> compatibleDetectors = layer.getDetLayer()->compatibleDets(trajectory, propagator, est);

  // You have to sort the simHits in the order they occur!

  // The old algorithm (sorting by distance to IP) doesn't seem to make sense to me (what if particle moves inwards?)

  // Detector layers have to be sorted by proximity to particle.position
  // Propagate particle backwards a bit to make sure it's outside any components
  std::map<double, std::unique_ptr<PSimHit>> distAndHits;
  // Position relative to which the hits should be sorted
  GlobalPoint positionOutside(particle.position().x() - particle.momentum().x() / particle.momentum().P() * 10.,
                              particle.position().y() - particle.momentum().y() / particle.momentum().P() * 10.,
                              particle.position().z() - particle.momentum().z() / particle.momentum().P() * 10.);

  // FastSim: cheat tracking -> assign hits to closest charged daughter if particle decays
  int pdgId = particle.pdgId();
  int simTrackId =
      particle.getMotherSimTrackIndex() >= 0 ? particle.getMotherSimTrackIndex() : particle.simTrackIndex();

  // loop over the compatible detectors
  for (const auto& detectorWithState : compatibleDetectors) {
    const GeomDet& detector = *detectorWithState.first;
    const TrajectoryStateOnSurface& particleState = detectorWithState.second;

    // if the detector has no components
    if (detector.isLeaf()) {
      std::pair<double, std::unique_ptr<PSimHit>> hitPair = createHitOnDetector(particleState,
                                                                                pdgId,
                                                                                layer.getThickness(particle.position()),
                                                                                energyDeposit,
                                                                                simTrackId,
                                                                                detector,
                                                                                positionOutside);
      if (hitPair.second) {
        distAndHits.insert(distAndHits.end(), std::move(hitPair));
      }
    } else {
      // if the detector has components
      for (const auto component : detector.components()) {
        std::pair<double, std::unique_ptr<PSimHit>> hitPair =
            createHitOnDetector(particleState,
                                pdgId,
                                layer.getThickness(particle.position()),
                                energyDeposit,
                                simTrackId,
                                *component,
                                positionOutside);
        if (hitPair.second) {
          distAndHits.insert(distAndHits.end(), std::move(hitPair));
        }
      }
    }
  }

  // Fill simHitContainer
  for (std::map<double, std::unique_ptr<PSimHit>>::const_iterator it = distAndHits.begin(); it != distAndHits.end();
       it++) {
    simHitContainer_->push_back(*(it->second));
  }
}

// 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)
std::pair<double, std::unique_ptr<PSimHit>> fastsim::TrackerSimHitProducer::createHitOnDetector(
    const TrajectoryStateOnSurface& particle,
    int pdgId,
    double layerThickness,
    double eLoss,
    int simTrackId,
    const GeomDet& detector,
    GlobalPoint& refPos) {
  // determine position and momentum of particle in the coordinate system of the detector
  LocalPoint localPosition;
  LocalVector localMomentum;

  // if the particle is close enough, no further propagation is needed
  if (fabs(detector.toLocal(particle.globalPosition()).z()) < onSurfaceTolerance_) {
    localPosition = particle.localPosition();
    localMomentum = particle.localMomentum();
  }
  // else, propagate
  else {
    // find crossing of particle with detector layer
    HelixArbitraryPlaneCrossing crossing(particle.globalPosition().basicVector(),
                                         particle.globalMomentum().basicVector(),
                                         particle.transverseCurvature(),
                                         anyDirection);
    std::pair<bool, double> path = crossing.pathLength(detector.surface());
    // case propagation succeeds
    if (path.first) {
      localPosition = detector.toLocal(GlobalPoint(crossing.position(path.second)));
      localMomentum = detector.toLocal(GlobalVector(crossing.direction(path.second)));
      localMomentum = localMomentum.unit() * particle.localMomentum().mag();
    }
    // case propagation fails
    else {
      return std::pair<double, std::unique_ptr<PSimHit>>(0, std::unique_ptr<PSimHit>(nullptr));
    }
  }

  // find entry and exit point of particle in detector
  const Plane& detectorPlane = detector.surface();
  double halfThick = 0.5 * detectorPlane.bounds().thickness();
  double pZ = localMomentum.z();
  LocalPoint entry = localPosition + (-halfThick / pZ) * localMomentum;
  LocalPoint exit = localPosition + (halfThick / pZ) * localMomentum;
  double tof = particle.globalPosition().mag() / fastsim::Constants::speedOfLight;  // in nanoseconds

  // make sure the simhit is physically on the module
  double boundX = detectorPlane.bounds().width() / 2.;
  double boundY = detectorPlane.bounds().length() / 2.;
  // Special treatment for TID and TEC trapeziodal modules
  unsigned subdet = DetId(detector.geographicalId()).subdetId();
  if (subdet == 4 || subdet == 6)
    boundX *= 1. - localPosition.y() / detectorPlane.position().perp();
  if (fabs(localPosition.x()) > boundX || fabs(localPosition.y()) > boundY) {
    return std::pair<double, std::unique_ptr<PSimHit>>(0, std::unique_ptr<PSimHit>(nullptr));
  }

  // Create the hit: the energy loss rescaled to the module thickness
  // Total thickness is in radiation lengths, 1 radlen = 9.36 cm
  // Sensitive module thickness is about 30 microns larger than
  // the module thickness itself
  eLoss *= (2. * halfThick - 0.003) / (9.36 * layerThickness);

  // Position of the hit in global coordinates
  GlobalPoint hitPos(detector.surface().toGlobal(localPosition));

  return std::pair<double, std::unique_ptr<PSimHit>>((hitPos - refPos).mag(),
                                                     std::make_unique<PSimHit>(entry,
                                                                               exit,
                                                                               localMomentum.mag(),
                                                                               tof,
                                                                               eLoss,
                                                                               pdgId,
                                                                               detector.geographicalId().rawId(),
                                                                               simTrackId,
                                                                               localMomentum.theta(),
                                                                               localMomentum.phi()));
}

DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::TrackerSimHitProducer, "fastsim::TrackerSimHitProducer");