JetAssociationProducers/src/TrackExtrapolator.cc

0001 #include "RecoJets/JetAssociationProducers/interface/TrackExtrapolator.h"
0002
0003 #include <vector>
0004
0005 //
0006 // constructors and destructor
0007 //
0008 TrackExtrapolator::TrackExtrapolator(const edm::ParameterSet& iConfig)
0009     : tracksSrc_(consumes(iConfig.getParameter<edm::InputTag>("trackSrc"))),
0010       fieldToken_(esConsumes()),
0011       propagatorToken_(esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAlong"))),
0012       ecalDetIdAssociatorToken_(esConsumes(edm::ESInputTag("", "EcalDetIdAssociator"))),
0013       trackQuality_(reco::TrackBase::qualityByName(iConfig.getParameter<std::string>("trackQuality"))) {
0014   if (trackQuality_ == reco::TrackBase::undefQuality) {  // we have a problem
0015     throw cms::Exception("InvalidInput") << "Unknown trackQuality value '"
0016                                          << iConfig.getParameter<std::string>("trackQuality")
0017                                          << "'. See possible values in 'reco::TrackBase::qualityByName'";
0018   }
0019
0020   produces<std::vector<reco::TrackExtrapolation>>();
0021 }
0022
0023 TrackExtrapolator::~TrackExtrapolator() {}
0024
0025 //
0026 // member functions
0027 //
0028
0029 // ------------ method called on each new Event  ------------
0030 void TrackExtrapolator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0031   // get stuff from Event Setup
0032   auto const& field = iSetup.getData(fieldToken_);
0033   auto const& propagator = iSetup.getData(propagatorToken_);
0034   FiducialVolume const& ecalvolume = iSetup.getData(ecalDetIdAssociatorToken_).volume();
0035
0036   // get stuff from Event
0037   edm::Handle<reco::TrackCollection> tracks_h;
0038   iEvent.getByToken(tracksSrc_, tracks_h);
0039
0040   auto extrapolations = std::make_unique<std::vector<reco::TrackExtrapolation>>();
0041
0042   // Get list of tracks we want to extrapolate
0043   std::vector<reco::TrackRef> goodTracks;
0044   for (reco::TrackCollection::const_iterator trkBegin = tracks_h->begin(), trkEnd = tracks_h->end(), itrk = trkBegin;
0045        itrk != trkEnd;
0046        ++itrk) {
0047     reco::TrackBase::TrackQuality trackQuality = reco::TrackBase::TrackQuality(trackQuality_);
0048
0049     // Cut on track quality
0050     if (itrk->quality(trackQuality)) {
0051       goodTracks.push_back(reco::TrackRef(tracks_h, itrk - trkBegin));
0052     }
0053   }
0054   std::vector<reco::TrackBase::Point> vresultPos(1);
0055   std::vector<reco::TrackBase::Vector> vresultMom(1);
0056
0057   // Now loop through the list of tracks and extrapolate them
0058   for (std::vector<reco::TrackRef>::const_iterator trkBegin = goodTracks.begin(),
0059                                                    trkEnd = goodTracks.end(),
0060                                                    itrk = trkBegin;
0061        itrk != trkEnd;
0062        ++itrk) {
0063     if (propagateTrackToVolume(**itrk, field, propagator, ecalvolume, vresultPos[0], vresultMom[0])) {
0064       extrapolations->push_back(reco::TrackExtrapolation(*itrk, vresultPos, vresultMom));
0065     }
0066   }
0067   iEvent.put(std::move(extrapolations));
0068 }
0069
0070 // -----------------------------------------------------------------------------
0071 //
0072 bool TrackExtrapolator::propagateTrackToVolume(const reco::Track& fTrack,
0073                                                const MagneticField& fField,
0074                                                const Propagator& fPropagator,
0075                                                const FiducialVolume& volume,
0076                                                reco::TrackBase::Point& resultPos,
0077                                                reco::TrackBase::Vector& resultMom) {
0078   GlobalPoint trackPosition(fTrack.vx(), fTrack.vy(), fTrack.vz());   // reference point
0079   GlobalVector trackMomentum(fTrack.px(), fTrack.py(), fTrack.pz());  // reference momentum
0080   if (fTrack.extra().isAvailable()) {                                 // use outer point information, if available
0081     trackPosition = GlobalPoint(fTrack.outerX(), fTrack.outerY(), fTrack.outerZ());
0082     trackMomentum = GlobalVector(fTrack.outerPx(), fTrack.outerPy(), fTrack.outerPz());
0083   }
0084
0085   GlobalTrajectoryParameters trackParams(trackPosition, trackMomentum, fTrack.charge(), &fField);
0086   FreeTrajectoryState trackState(trackParams);
0087
0088   TrajectoryStateOnSurface propagatedInfo = fPropagator.propagate(
0089       trackState, *Cylinder::build(volume.minR(), Surface::PositionType(0, 0, 0), Surface::RotationType()));
0090
0091   // if the track went through either side of the endcaps, repropagate the track
0092   double minz = volume.minZ();
0093   if (propagatedInfo.isValid() && propagatedInfo.globalPosition().z() > minz) {
0094     propagatedInfo =
0095         fPropagator.propagate(trackState, *Plane::build(Surface::PositionType(0, 0, minz), Surface::RotationType()));
0096
0097   } else if (propagatedInfo.isValid() && propagatedInfo.globalPosition().z() < -minz) {
0098     propagatedInfo =
0099         fPropagator.propagate(trackState, *Plane::build(Surface::PositionType(0, 0, -minz), Surface::RotationType()));
0100   }
0101
0102   if (propagatedInfo.isValid()) {
0103     resultPos = propagatedInfo.globalPosition();
0104     resultMom = propagatedInfo.globalMomentum();
0105     return true;
0106   } else {
0107     return false;
0108   }
0109 }
0110
0111 //define this as a plug-in
0112 DEFINE_FWK_MODULE(TrackExtrapolator);