JetAssociationAlgorithms/src/JetTracksAssociationDRCalo.cc

0001 // Associate jets with tracks by simple "dR" criteria
0002 // Fedor Ratnikov (UMd), Aug. 28, 2007
0003
0004 #include "RecoJets/JetAssociationAlgorithms/interface/JetTracksAssociationDRCalo.h"
0005
0006 #include "DataFormats/JetReco/interface/Jet.h"
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008
0009 #include "DataFormats/Math/interface/deltaR.h"
0010 #include "DataFormats/Math/interface/Vector3D.h"
0011
0012 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0013 #include "DataFormats/GeometrySurface/interface/Plane.h"
0014
0015 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0016 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0017 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0018 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0019 #include "MagneticField/Engine/interface/MagneticField.h"
0020
0021 namespace {
0022   // basic geometry constants, imported from Geometry/HcalTowerAlgo/src/CaloTowerHardcodeGeometryLoader.cc
0023   const double rBarrel = 129.;
0024   const double zEndcap = 320.;
0025   const double zVF = 1100.;
0026   const double rEndcapMin = zEndcap * tan(2 * atan(exp(-3.)));
0027   const double rVFMin = zEndcap * tan(2 * atan(exp(-5.191)));
0028
0029   struct ImpactPoint {
0030     unsigned index;
0031     double eta;
0032     double phi;
0033   };
0034
0035   GlobalPoint propagateTrackToCalo(const reco::Track& fTrack,
0036                                    const MagneticField& fField,
0037                                    const Propagator& fPropagator) {
0038     GlobalPoint trackPosition(fTrack.vx(), fTrack.vy(), fTrack.vz());   // reference point
0039     GlobalVector trackMomentum(fTrack.px(), fTrack.py(), fTrack.pz());  // reference momentum
0040     if (fTrack.extra().isAvailable()) {                                 // use outer point information, if available
0041       trackPosition = GlobalPoint(fTrack.outerX(), fTrack.outerY(), fTrack.outerZ());
0042       trackMomentum = GlobalVector(fTrack.outerPx(), fTrack.outerPy(), fTrack.outerPz());
0043     }
0044     //     std::cout << "propagateTrackToCalo-> start propagating track"
0045     //        << " x/y/z: " << trackPosition.x() << '/' << trackPosition.y() << '/' << trackPosition.z()
0046     //        << ", pt/eta/phi: " << trackMomentum.perp() << '/' << trackMomentum.eta() << '/' << trackMomentum.barePhi()
0047     //        << std::endl;
0048     GlobalTrajectoryParameters trackParams(trackPosition, trackMomentum, fTrack.charge(), &fField);
0049     FreeTrajectoryState trackState(trackParams);
0050
0051     // first propagate to barrel
0052     TrajectoryStateOnSurface propagatedInfo = fPropagator.propagate(
0053         trackState, *Cylinder::build(rBarrel, Surface::PositionType(0, 0, 0), Surface::RotationType()));
0054     if (propagatedInfo.isValid()) {
0055       GlobalPoint result(propagatedInfo.globalPosition());
0056       if (fabs(result.z()) < zEndcap) {
0057         //  std::cout << "propagateTrackToCalo-> propagated to barrel:"
0058         //        << " x/y/z/r: " << result.x() << '/' << result.y() << '/' << result.z() << '/' << result.perp()
0059         //        << std::endl;
0060         return result;
0061       }
0062     }
0063
0064     // failed with barrel, try endcap
0065     double zTarget = trackMomentum.z() > 0 ? zEndcap : -zEndcap;
0066     propagatedInfo =
0067         fPropagator.propagate(trackState, *Plane::build(Surface::PositionType(0, 0, zTarget), Surface::RotationType()));
0068     if (propagatedInfo.isValid()) {
0069       GlobalPoint result(propagatedInfo.globalPosition());
0070       if (fabs(result.perp()) > rEndcapMin) {
0071         //  std::cout << "propagateTrackToCalo-> propagated to endcap:"
0072         //        << " x/y/z/r: " << result.x() << '/' << result.y() << '/' << result.z() << '/' << result.perp()
0073         //        << std::endl;
0074         return result;
0075       }
0076     }
0077     // failed with endcap, try VF
0078     zTarget = trackMomentum.z() > 0 ? zVF : -zVF;
0079     propagatedInfo =
0080         fPropagator.propagate(trackState, *Plane::build(Surface::PositionType(0, 0, zTarget), Surface::RotationType()));
0081     if (propagatedInfo.isValid()) {
0082       GlobalPoint result(propagatedInfo.globalPosition());
0083       if (fabs(result.perp()) > rVFMin) {
0084         //  std::cout << "propagateTrackToCalo-> propagated to VF:"
0085         //        << " x/y/z/r: " << result.x() << '/' << result.y() << '/' << result.z() << '/' << result.perp()
0086         //        << std::endl;
0087         return result;
0088       }
0089     }
0090     // no luck
0091     //     std::cout << "propagateTrackToCalo-> failed to propagate track to calorimeter" << std::endl;
0092     return GlobalPoint(0, 0, 0);
0093   }
0094 }  // namespace
0095
0096 JetTracksAssociationDRCalo::JetTracksAssociationDRCalo(double fDr) : mDeltaR2Threshold(fDr * fDr) {}
0097
0098 void JetTracksAssociationDRCalo::produce(reco::JetTracksAssociation::Container* fAssociation,
0099                                          const std::vector<edm::RefToBase<reco::Jet> >& fJets,
0100                                          const std::vector<reco::TrackRef>& fTracks,
0101                                          const MagneticField& fField,
0102                                          const Propagator& fPropagator) const {
0103   // cache track parameters
0104   std::vector<ImpactPoint> impacts;
0105   for (unsigned t = 0; t < fTracks.size(); ++t) {
0106     GlobalPoint impact = propagateTrackToCalo(*(fTracks[t]), fField, fPropagator);
0107     if (impact.mag() > 0) {  // successful extrapolation
0108       ImpactPoint goodTrack;
0109       goodTrack.index = t;
0110       goodTrack.eta = impact.eta();
0111       goodTrack.phi = impact.barePhi();
0112       impacts.push_back(goodTrack);
0113     }
0114   }
0115
0116   for (unsigned j = 0; j < fJets.size(); ++j) {
0117     reco::TrackRefVector assoTracks;
0118     const reco::Jet* jet = &*(fJets[j]);
0119     double jetEta = jet->eta();
0120     double jetPhi = jet->phi();
0121     for (unsigned t = 0; t < impacts.size(); ++t) {
0122       double dR2 = deltaR2(jetEta, jetPhi, impacts[t].eta, impacts[t].phi);
0123       if (dR2 < mDeltaR2Threshold)
0124         assoTracks.push_back(fTracks[impacts[t].index]);
0125     }
0126     reco::JetTracksAssociation::setValue(fAssociation, fJets[j], assoTracks);
0127   }
0128 }
0129
0130 math::XYZPoint JetTracksAssociationDRCalo::propagateTrackToCalorimeter(const reco::Track& fTrack,
0131                                                                        const MagneticField& fField,
0132                                                                        const Propagator& fPropagator) {
0133   GlobalPoint result(propagateTrackToCalo(fTrack, fField, fPropagator));
0134   return math::XYZPoint(result.x(), result.y(), result.z());
0135 }