File indexing completed on 2024-04-06 12:25:28
0001
0002
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
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());
0039 GlobalVector trackMomentum(fTrack.px(), fTrack.py(), fTrack.pz());
0040 if (fTrack.extra().isAvailable()) {
0041 trackPosition = GlobalPoint(fTrack.outerX(), fTrack.outerY(), fTrack.outerZ());
0042 trackMomentum = GlobalVector(fTrack.outerPx(), fTrack.outerPy(), fTrack.outerPz());
0043 }
0044
0045
0046
0047
0048 GlobalTrajectoryParameters trackParams(trackPosition, trackMomentum, fTrack.charge(), &fField);
0049 FreeTrajectoryState trackState(trackParams);
0050
0051
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
0058
0059
0060 return result;
0061 }
0062 }
0063
0064
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
0072
0073
0074 return result;
0075 }
0076 }
0077
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
0085
0086
0087 return result;
0088 }
0089 }
0090
0091
0092 return GlobalPoint(0, 0, 0);
0093 }
0094 }
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
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) {
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 }