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;
}