EgammaPhotonAlgos/src/ConversionTrackEcalImpactPoint.cc

0001 #include <iostream>
0002 #include <vector>
0003 #include <memory>
0004 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackEcalImpactPoint.h"
0005 // Framework
0006 //
0007
0008 //
0009 #include <vector>
0010 #include <map>
0011
0012 constexpr float epsilon = 0.001;
0013
0014 /** Hard-wired numbers defining the surfaces on which the crystal front faces lie. */
0015 constexpr float barrelRadius = 129.f;       // p81, p50, ECAL TDR
0016 constexpr float barrelHalfLength = 270.9f;  // p81, p50, ECAL TDR
0017 constexpr float endcapRadius = 171.1f;      // fig 3.26, p81, ECAL TDR
0018 constexpr float endcapZ = 320.5f;           // fig 3.26, p81, ECAL TDR
0019
0020 static BoundCylinder* initBarrel() {
0021   Surface::RotationType rot;  // unit rotation matrix
0022
0023   return new Cylinder(
0024       barrelRadius,
0025       Surface::PositionType(0, 0, 0),
0026       rot,
0027       new SimpleCylinderBounds(barrelRadius - epsilon, barrelRadius + epsilon, -barrelHalfLength, barrelHalfLength));
0028 }
0029
0030 static BoundDisk* initNegative() {
0031   Surface::RotationType rot;  // unit rotation matrix
0032   return new BoundDisk(
0033       Surface::PositionType(0, 0, -endcapZ), rot, new SimpleDiskBounds(0, endcapRadius, -epsilon, epsilon));
0034 }
0035
0036 static BoundDisk* initPositive() {
0037   Surface::RotationType rot;  // unit rotation matrix
0038
0039   return new BoundDisk(
0040       Surface::PositionType(0, 0, endcapZ), rot, new SimpleDiskBounds(0, endcapRadius, -epsilon, epsilon));
0041 }
0042
0043 const ReferenceCountingPointer<BoundCylinder> ConversionTrackEcalImpactPoint::theBarrel_ = initBarrel();
0044 const ReferenceCountingPointer<BoundDisk> ConversionTrackEcalImpactPoint::theNegativeEtaEndcap_ = initNegative();
0045 const ReferenceCountingPointer<BoundDisk> ConversionTrackEcalImpactPoint::thePositiveEtaEndcap_ = initPositive();
0046
0047 ConversionTrackEcalImpactPoint::ConversionTrackEcalImpactPoint(const MagneticField* field)
0048     :
0049
0050       theMF_(field) {
0051   forwardPropagator_ = new PropagatorWithMaterial(dir_ = alongMomentum, 0.000511, theMF_);
0052 }
0053
0054 ConversionTrackEcalImpactPoint::~ConversionTrackEcalImpactPoint() { delete forwardPropagator_; }
0055
0056 std::vector<math::XYZPointF> ConversionTrackEcalImpactPoint::find(
0057     const std::vector<reco::TransientTrack>& tracks, const edm::Handle<edm::View<reco::CaloCluster> >& bcHandle) {
0058   std::vector<math::XYZPointF> result;
0059   //
0060   matchingBC_.clear();
0061
0062   std::vector<reco::CaloClusterPtr> matchingBC(2);
0063
0064   //
0065
0066   int iTrk = 0;
0067   for (auto const& track : tracks) {
0068     math::XYZPointF ecalImpactPosition(0., 0., 0.);
0069     const TrajectoryStateOnSurface myTSOS = track.innermostMeasurementState();
0070     if (!(myTSOS.isValid()))
0071       continue;
0072
0073     stateAtECAL_ = forwardPropagator_->propagate(myTSOS, barrel());
0074
0075     if (!stateAtECAL_.isValid() || (stateAtECAL_.isValid() && fabs(stateAtECAL_.globalPosition().eta()) > 1.479)) {
0076       if (track.innermostMeasurementState().globalPosition().eta() > 0.) {
0077         stateAtECAL_ = forwardPropagator_->propagate(myTSOS, positiveEtaEndcap());
0078
0079       } else {
0080         stateAtECAL_ = forwardPropagator_->propagate(myTSOS, negativeEtaEndcap());
0081       }
0082     }
0083
0084     if (stateAtECAL_.isValid())
0085       ecalImpactPosition = stateAtECAL_.globalPosition();
0086
0087     result.push_back(ecalImpactPosition);
0088
0089     if (stateAtECAL_.isValid()) {
0090       int goodBC = 0;
0091       float bcDistanceToTrack = 9999;
0092       reco::CaloClusterPtr matchedBCItr;
0093       int ibc = 0;
0094       goodBC = 0;
0095
0096       for (auto const& bc : *bcHandle) {
0097         float dEta = bc.position().eta() - ecalImpactPosition.eta();
0098         float dPhi = bc.position().phi() - ecalImpactPosition.phi();
0099         if (sqrt(dEta * dEta + dPhi * dPhi) < bcDistanceToTrack) {
0100           goodBC = ibc;
0101           bcDistanceToTrack = sqrt(dEta * dEta + dPhi * dPhi);
0102         }
0103         ibc++;
0104       }
0105
0106       matchingBC[iTrk] = bcHandle->ptrAt(goodBC);
0107     }
0108
0109     iTrk++;
0110   }
0111
0112   matchingBC_ = matchingBC;
0113
0114   return result;
0115 }