ConversionTrackEcalImpactPoint.cc

Go to the documentation of this file.

#include <iostream>
#include <vector>
#include <memory>
#include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackEcalImpactPoint.h"
// Framework
//
 
//
#include <vector>
#include <map>
 
constexpr float epsilon = 0.001;
 
constexpr float barrelRadius = 129.f;       // p81, p50, ECAL TDR
constexpr float barrelHalfLength = 270.9f;  // p81, p50, ECAL TDR
constexpr float endcapRadius = 171.1f;      // fig 3.26, p81, ECAL TDR
constexpr float endcapZ = 320.5f;           // fig 3.26, p81, ECAL TDR
 
static BoundCylinder* initBarrel() {
  Surface::RotationType rot;  // unit rotation matrix
 
  return new Cylinder(
      barrelRadius,
      Surface::PositionType(0, 0, 0),
      rot,
      new SimpleCylinderBounds(barrelRadius - epsilon, barrelRadius + epsilon, -barrelHalfLength, barrelHalfLength));
}
 
static BoundDisk* initNegative() {
  Surface::RotationType rot;  // unit rotation matrix
  return new BoundDisk(
      Surface::PositionType(0, 0, -endcapZ), rot, new SimpleDiskBounds(0, endcapRadius, -epsilon, epsilon));
}
 
static BoundDisk* initPositive() {
  Surface::RotationType rot;  // unit rotation matrix
 
  return new BoundDisk(
      Surface::PositionType(0, 0, endcapZ), rot, new SimpleDiskBounds(0, endcapRadius, -epsilon, epsilon));
}
 
const ReferenceCountingPointer<BoundCylinder> ConversionTrackEcalImpactPoint::theBarrel_ = initBarrel();
const ReferenceCountingPointer<BoundDisk> ConversionTrackEcalImpactPoint::theNegativeEtaEndcap_ = initNegative();
const ReferenceCountingPointer<BoundDisk> ConversionTrackEcalImpactPoint::thePositiveEtaEndcap_ = initPositive();
 
ConversionTrackEcalImpactPoint::ConversionTrackEcalImpactPoint(const MagneticField* field)
    :
 
      theMF_(field) {
  forwardPropagator_ = new PropagatorWithMaterial(dir_ = alongMomentum, 0.000511, theMF_);
}
 
ConversionTrackEcalImpactPoint::~ConversionTrackEcalImpactPoint() { delete forwardPropagator_; }
 
std::vector<math::XYZPointF> ConversionTrackEcalImpactPoint::find(
    const std::vector<reco::TransientTrack>& tracks, const edm::Handle<edm::View<reco::CaloCluster> >& bcHandle) {
  std::vector<math::XYZPointF> result;
  //
  matchingBC_.clear();
 
  std::vector<reco::CaloClusterPtr> matchingBC(2);
 
  //
 
  int iTrk = 0;
  for (auto const& track : tracks) {
    math::XYZPointF ecalImpactPosition(0., 0., 0.);
    const TrajectoryStateOnSurface myTSOS = track.innermostMeasurementState();
    if (!(myTSOS.isValid()))
      continue;
 
    stateAtECAL_ = forwardPropagator_->propagate(myTSOS, barrel());
 
    if (!stateAtECAL_.isValid() || (stateAtECAL_.isValid() && fabs(stateAtECAL_.globalPosition().eta()) > 1.479)) {
      if (track.innermostMeasurementState().globalPosition().eta() > 0.) {
        stateAtECAL_ = forwardPropagator_->propagate(myTSOS, positiveEtaEndcap());
 
      } else {
        stateAtECAL_ = forwardPropagator_->propagate(myTSOS, negativeEtaEndcap());
      }
    }
 
    if (stateAtECAL_.isValid())
      ecalImpactPosition = stateAtECAL_.globalPosition();
 
    result.push_back(ecalImpactPosition);
 
    if (stateAtECAL_.isValid()) {
      int goodBC = 0;
      float bcDistanceToTrack = 9999;
      reco::CaloClusterPtr matchedBCItr;
      int ibc = 0;
      goodBC = 0;
 
      for (auto const& bc : *bcHandle) {
        float dEta = bc.position().eta() - ecalImpactPosition.eta();
        float dPhi = bc.position().phi() - ecalImpactPosition.phi();
        if (sqrt(dEta * dEta + dPhi * dPhi) < bcDistanceToTrack) {
          goodBC = ibc;
          bcDistanceToTrack = sqrt(dEta * dEta + dPhi * dPhi);
        }
        ibc++;
      }
 
      matchingBC[iTrk] = bcHandle->ptrAt(goodBC);
    }
 
    iTrk++;
  }
 
  matchingBC_ = matchingBC;
 
  return result;
}