Tracking/plugins/TrackCandidateProducer.cc

0001 // system
0002 #include <memory>
0003 #include <vector>
0004 #include <map>
0005
0006 // framework
0007 #include "FWCore/Framework/interface/stream/EDProducer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 // data format
0016 #include "DataFormats/Common/interface/Handle.h"
0017 #include "DataFormats/Common/interface/OwnVector.h"
0018 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0019 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0020 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
0021 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHitCollection.h"
0022 #include "SimDataFormats/Track/interface/SimTrack.h"
0023 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0024
0025 // geometry / magnetic field / propagation
0026 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0027 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0028 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0029 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0030 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0031 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0032 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0033 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0034 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0035
0036 // fastsim
0037 #include "FastSimulation/Tracking/interface/TrackingLayer.h"
0038 #include "FastSimulation/Tracking/interface/FastTrackingUtilities.h"
0039 #include "FastSimulation/Tracking/interface/FastTrackerRecHitSplitter.h"
0040 #include "FastSimulation/Tracking/interface/SeedMatcher.h"
0041
0042 class TrackCandidateProducer : public edm::stream::EDProducer<> {
0043 public:
0044   explicit TrackCandidateProducer(const edm::ParameterSet& conf);
0045
0046   void produce(edm::Event&, const edm::EventSetup&) override;
0047
0048 private:
0049   // tokens & labels
0050   edm::EDGetTokenT<edm::View<TrajectorySeed> > seedToken;
0051   edm::EDGetTokenT<FastTrackerRecHitCombinationCollection> recHitCombinationsToken;
0052   edm::EDGetTokenT<std::vector<bool> > hitMasksToken;
0053   edm::EDGetTokenT<edm::SimTrackContainer> simTrackToken;
0054   std::string propagatorLabel;
0055   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldESToken_;
0056   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeometryESToken_;
0057   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopologyESToken_;
0058   const edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorESToken_;
0059
0060   // other data
0061   bool rejectOverlaps;
0062   bool splitHits;
0063   FastTrackerRecHitSplitter hitSplitter;
0064   double maxSeedMatchEstimator;
0065 };
0066
0067 TrackCandidateProducer::TrackCandidateProducer(const edm::ParameterSet& conf)
0068     : propagatorLabel(conf.getParameter<std::string>("propagator")),
0069       magneticFieldESToken_(esConsumes()),
0070       trackerGeometryESToken_(esConsumes()),
0071       trackerTopologyESToken_(esConsumes()),
0072       propagatorESToken_(esConsumes(edm::ESInputTag("", propagatorLabel))),
0073       hitSplitter() {
0074   // produces
0075   produces<TrackCandidateCollection>();
0076
0077   // consumes
0078   if (conf.exists("hitMasks")) {
0079     hitMasksToken = consumes<std::vector<bool> >(conf.getParameter<edm::InputTag>("hitMasks"));
0080   }
0081   seedToken = consumes<edm::View<TrajectorySeed> >(conf.getParameter<edm::InputTag>("src"));
0082   recHitCombinationsToken =
0083       consumes<FastTrackerRecHitCombinationCollection>(conf.getParameter<edm::InputTag>("recHitCombinations"));
0084   simTrackToken = consumes<edm::SimTrackContainer>(conf.getParameter<edm::InputTag>("simTracks"));
0085
0086   // other parameters
0087   maxSeedMatchEstimator = conf.getUntrackedParameter<double>("maxSeedMatchEstimator", 0);
0088   rejectOverlaps = conf.getParameter<bool>("OverlapCleaning");
0089   splitHits = conf.getParameter<bool>("SplitHits");
0090 }
0091
0092 void TrackCandidateProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0093   // get records
0094   auto const& magneticField = es.getData(magneticFieldESToken_);
0095   auto const& trackerGeometry = es.getData(trackerGeometryESToken_);
0096   auto const& trackerTopology = es.getData(trackerTopologyESToken_);
0097   auto const& propagator = es.getData(propagatorESToken_);
0098
0099   // get products
0100   edm::Handle<edm::View<TrajectorySeed> > seeds;
0101   e.getByToken(seedToken, seeds);
0102
0103   edm::Handle<FastTrackerRecHitCombinationCollection> recHitCombinations;
0104   e.getByToken(recHitCombinationsToken, recHitCombinations);
0105
0106   edm::Handle<edm::SimTrackContainer> simTracks;
0107   e.getByToken(simTrackToken, simTracks);
0108
0109   edm::Handle<std::vector<bool> > hitMasks;
0110   if (!hitMasksToken.isUninitialized()) {
0111     e.getByToken(hitMasksToken, hitMasks);
0112   }
0113
0114   // output collection
0115   std::unique_ptr<TrackCandidateCollection> output(new TrackCandidateCollection);
0116
0117   // loop over the seeds
0118   for (unsigned seedIndex = 0; seedIndex < seeds->size(); ++seedIndex) {
0119     const TrajectorySeed seed = (*seeds)[seedIndex];
0120     std::vector<int32_t> recHitCombinationIndices;
0121
0122     // match hitless seeds to simTracks and find corresponding recHitCombination
0123     if (seed.nHits() == 0) {
0124       recHitCombinationIndices = SeedMatcher::matchRecHitCombinations(
0125           seed, *recHitCombinations, *simTracks, maxSeedMatchEstimator, propagator, magneticField, trackerGeometry);
0126     }
0127     // for normal seeds, retrieve the corresponding recHitCombination from the seed hits
0128     else {
0129       int32_t icomb = fastTrackingUtilities::getRecHitCombinationIndex(seed);
0130       recHitCombinationIndices.push_back(icomb);
0131     }
0132
0133     // loop over the matched recHitCombinations
0134     for (auto icomb : recHitCombinationIndices) {
0135       if (icomb < 0 || unsigned(icomb) >= recHitCombinations->size()) {
0136         throw cms::Exception("TrackCandidateProducer")
0137             << " found seed with recHitCombination out or range: " << icomb << std::endl;
0138       }
0139       const FastTrackerRecHitCombination& recHitCombination = (*recHitCombinations)[icomb];
0140
0141       // container for select hits
0142       std::vector<const FastTrackerRecHit*> selectedRecHits;
0143
0144       // add the seed hits
0145       for (auto const& recHit : seed.recHits()) {
0146         selectedRecHits.push_back(static_cast<const FastTrackerRecHit*>(&recHit));
0147       }
0148
0149       // prepare to skip seed hits
0150       const FastTrackerRecHit* lastHitToSkip = nullptr;
0151       if (!selectedRecHits.empty()) {
0152         lastHitToSkip = selectedRecHits.back();
0153       }
0154
0155       // inOut or outIn tracking ?
0156       bool hitsAlongMomentum = (seed.direction() == alongMomentum);
0157
0158       // add hits from combination to hit selection
0159       for (unsigned hitIndex = hitsAlongMomentum ? 0 : recHitCombination.size() - 1;
0160            hitIndex < recHitCombination.size();
0161            hitsAlongMomentum ? ++hitIndex : --hitIndex) {
0162         const FastTrackerRecHit* selectedRecHit = recHitCombination[hitIndex].get();
0163
0164         // skip seed hits
0165         if (lastHitToSkip) {
0166           if (lastHitToSkip->sameId(selectedRecHit)) {
0167             lastHitToSkip = nullptr;
0168           }
0169           continue;
0170         }
0171
0172         // apply hit masking
0173         if (hitMasks.isValid() && fastTrackingUtilities::hitIsMasked(selectedRecHit, *hitMasks)) {
0174           continue;
0175         }
0176
0177         //  if overlap rejection is not switched on, accept all hits
0178         //  always accept the first hit
0179         //  also accept a hit if it is not on the layer of the previous hit
0180         if (!rejectOverlaps || selectedRecHits.empty() ||
0181             (TrackingLayer::createFromDetId(selectedRecHits.back()->geographicalId(), trackerTopology) !=
0182              TrackingLayer::createFromDetId(selectedRecHit->geographicalId(), trackerTopology))) {
0183           selectedRecHits.push_back(selectedRecHit);
0184         }
0185         //  else:
0186         //    overlap rejection is switched on
0187         //    the hit is on the same layer as the previous hit
0188         //  accept the one with smallest error
0189         else if (fastTrackingUtilities::hitLocalError(selectedRecHit) <
0190                  fastTrackingUtilities::hitLocalError(selectedRecHits.back())) {
0191           selectedRecHits.back() = selectedRecHit;
0192         }
0193       }
0194
0195       // split hits / store copies for the track candidate
0196       edm::OwnVector<TrackingRecHit> hitsForTrackCandidate;
0197       for (unsigned index = 0; index < selectedRecHits.size(); ++index) {
0198         if (splitHits) {
0199           // add split hits to splitSelectedRecHits
0200           hitSplitter.split(*selectedRecHits[index], hitsForTrackCandidate, hitsAlongMomentum);
0201         } else {
0202           hitsForTrackCandidate.push_back(selectedRecHits[index]->clone());
0203         }
0204       }
0205
0206       // set the recHitCombinationIndex
0207       fastTrackingUtilities::setRecHitCombinationIndex(hitsForTrackCandidate, icomb);
0208
0209       // create track candidate state
0210       //   1. get seed state (defined on the surface of the most outer hit)
0211       DetId seedDetId(seed.startingState().detId());
0212       const GeomDet* gdet = trackerGeometry.idToDet(seedDetId);
0213       TrajectoryStateOnSurface seedTSOS =
0214           trajectoryStateTransform::transientState(seed.startingState(), &(gdet->surface()), &magneticField);
0215       //   2. backPropagate the seedState to the surfuce of the most inner hit
0216       const GeomDet* initialLayer = trackerGeometry.idToDet(hitsForTrackCandidate.front().geographicalId());
0217       const TrajectoryStateOnSurface initialTSOS = propagator.propagate(seedTSOS, initialLayer->surface());
0218       //   3. check validity and transform
0219       if (!initialTSOS.isValid())
0220         continue;
0221       PTrajectoryStateOnDet PTSOD = trajectoryStateTransform::persistentState(
0222           initialTSOS, hitsForTrackCandidate.front().geographicalId().rawId());
0223
0224       // add track candidate to output collection
0225       output->push_back(
0226           TrackCandidate(hitsForTrackCandidate, seed, PTSOD, edm::RefToBase<TrajectorySeed>(seeds, seedIndex)));
0227     }
0228   }
0229
0230   // Save the track candidates
0231   e.put(std::move(output));
0232 }
0233
0234 #include "FWCore/Framework/interface/MakerMacros.h"
0235 DEFINE_FWK_MODULE(TrackCandidateProducer);