TrackCandidateProducer.cc

Go to the documentation of this file.

// system
#include <memory>
#include <vector>
#include <map>
 
// framework
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
 
// data format
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Common/interface/OwnVector.h"
#include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
#include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHitCollection.h"
#include "SimDataFormats/Track/interface/SimTrack.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
 
// geometry / magnetic field / propagation
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
 
// fastsim
#include "FastSimulation/Tracking/interface/TrackingLayer.h"
#include "FastSimulation/Tracking/interface/FastTrackingUtilities.h"
#include "FastSimulation/Tracking/interface/FastTrackerRecHitSplitter.h"
#include "FastSimulation/Tracking/interface/SeedMatcher.h"
 
class TrackCandidateProducer : public edm::stream::EDProducer<> {
public:
  explicit TrackCandidateProducer(const edm::ParameterSet& conf);
 
  void produce(edm::Event& e, const edm::EventSetup& es) override;
 
private:
  // tokens & labels
  edm::EDGetTokenT<edm::View<TrajectorySeed> > seedToken;
  edm::EDGetTokenT<FastTrackerRecHitCombinationCollection> recHitCombinationsToken;
  edm::EDGetTokenT<std::vector<bool> > hitMasksToken;
  edm::EDGetTokenT<edm::SimTrackContainer> simTrackToken;
  std::string propagatorLabel;
 
  // other data
  bool rejectOverlaps;
  bool splitHits;
  FastTrackerRecHitSplitter hitSplitter;
  double maxSeedMatchEstimator;
};
 
TrackCandidateProducer::TrackCandidateProducer(const edm::ParameterSet& conf) : hitSplitter() {
  // produces
  produces<TrackCandidateCollection>();
 
  // consumes
  if (conf.exists("hitMasks")) {
    hitMasksToken = consumes<std::vector<bool> >(conf.getParameter<edm::InputTag>("hitMasks"));
  }
  seedToken = consumes<edm::View<TrajectorySeed> >(conf.getParameter<edm::InputTag>("src"));
  recHitCombinationsToken =
      consumes<FastTrackerRecHitCombinationCollection>(conf.getParameter<edm::InputTag>("recHitCombinations"));
  simTrackToken = consumes<edm::SimTrackContainer>(conf.getParameter<edm::InputTag>("simTracks"));
 
  // other parameters
  maxSeedMatchEstimator = conf.getUntrackedParameter<double>("maxSeedMatchEstimator", 0);
  rejectOverlaps = conf.getParameter<bool>("OverlapCleaning");
  splitHits = conf.getParameter<bool>("SplitHits");
  propagatorLabel = conf.getParameter<std::string>("propagator");
}
 
void TrackCandidateProducer::produce(edm::Event& e, const edm::EventSetup& es) {
  // get records
  edm::ESHandle<MagneticField> magneticField;
  es.get<IdealMagneticFieldRecord>().get(magneticField);
 
  edm::ESHandle<TrackerGeometry> trackerGeometry;
  es.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
 
  edm::ESHandle<TrackerTopology> trackerTopology;
  es.get<TrackerTopologyRcd>().get(trackerTopology);
 
  edm::ESHandle<Propagator> propagator;
  es.get<TrackingComponentsRecord>().get(propagatorLabel, propagator);
 
  // get products
  edm::Handle<edm::View<TrajectorySeed> > seeds;
  e.getByToken(seedToken, seeds);
 
  edm::Handle<FastTrackerRecHitCombinationCollection> recHitCombinations;
  e.getByToken(recHitCombinationsToken, recHitCombinations);
 
  edm::Handle<edm::SimTrackContainer> simTracks;
  e.getByToken(simTrackToken, simTracks);
 
  edm::Handle<std::vector<bool> > hitMasks;
  if (!hitMasksToken.isUninitialized()) {
    e.getByToken(hitMasksToken, hitMasks);
  }
 
  // output collection
  std::unique_ptr<TrackCandidateCollection> output(new TrackCandidateCollection);
 
  // loop over the seeds
  for (unsigned seedIndex = 0; seedIndex < seeds->size(); ++seedIndex) {
    const TrajectorySeed seed = (*seeds)[seedIndex];
    std::vector<int32_t> recHitCombinationIndices;
 
    // match hitless seeds to simTracks and find corresponding recHitCombination
    if (seed.nHits() == 0) {
      recHitCombinationIndices = SeedMatcher::matchRecHitCombinations(
          seed, *recHitCombinations, *simTracks, maxSeedMatchEstimator, *propagator, *magneticField, *trackerGeometry);
    }
    // for normal seeds, retrieve the corresponding recHitCombination from the seed hits
    else {
      int32_t icomb = fastTrackingUtilities::getRecHitCombinationIndex(seed);
      recHitCombinationIndices.push_back(icomb);
    }
 
    // loop over the matched recHitCombinations
    for (auto icomb : recHitCombinationIndices) {
      if (icomb < 0 || unsigned(icomb) >= recHitCombinations->size()) {
        throw cms::Exception("TrackCandidateProducer")
            << " found seed with recHitCombination out or range: " << icomb << std::endl;
      }
      const FastTrackerRecHitCombination& recHitCombination = (*recHitCombinations)[icomb];
 
      // container for select hits
      std::vector<const FastTrackerRecHit*> selectedRecHits;
 
      // add the seed hits
      TrajectorySeed::range seedHitRange = seed.recHits();  //Hits in a seed
      for (TrajectorySeed::const_iterator ihit = seedHitRange.first; ihit != seedHitRange.second; ++ihit) {
        selectedRecHits.push_back(static_cast<const FastTrackerRecHit*>(&*ihit));
      }
 
      // prepare to skip seed hits
      const FastTrackerRecHit* lastHitToSkip = nullptr;
      if (!selectedRecHits.empty()) {
        lastHitToSkip = selectedRecHits.back();
      }
 
      // inOut or outIn tracking ?
      bool hitsAlongMomentum = (seed.direction() == alongMomentum);
 
      // add hits from combination to hit selection
      for (unsigned hitIndex = hitsAlongMomentum ? 0 : recHitCombination.size() - 1;
           hitIndex < recHitCombination.size();
           hitsAlongMomentum ? ++hitIndex : --hitIndex) {
        const FastTrackerRecHit* selectedRecHit = recHitCombination[hitIndex].get();
 
        // skip seed hits
        if (lastHitToSkip) {
          if (lastHitToSkip->sameId(selectedRecHit)) {
            lastHitToSkip = nullptr;
          }
          continue;
        }
 
        // apply hit masking
        if (hitMasks.isValid() && fastTrackingUtilities::hitIsMasked(selectedRecHit, *hitMasks)) {
          continue;
        }
 
        //  if overlap rejection is not switched on, accept all hits
        //  always accept the first hit
        //  also accept a hit if it is not on the layer of the previous hit
        if (!rejectOverlaps || selectedRecHits.empty() ||
            (TrackingLayer::createFromDetId(selectedRecHits.back()->geographicalId(), *trackerTopology.product()) !=
             TrackingLayer::createFromDetId(selectedRecHit->geographicalId(), *trackerTopology.product()))) {
          selectedRecHits.push_back(selectedRecHit);
        }
        //  else:
        //    overlap rejection is switched on
        //    the hit is on the same layer as the previous hit
        //  accept the one with smallest error
        else if (fastTrackingUtilities::hitLocalError(selectedRecHit) <
                 fastTrackingUtilities::hitLocalError(selectedRecHits.back())) {
          selectedRecHits.back() = selectedRecHit;
        }
      }
 
      // split hits / store copies for the track candidate
      edm::OwnVector<TrackingRecHit> hitsForTrackCandidate;
      for (unsigned index = 0; index < selectedRecHits.size(); ++index) {
        if (splitHits) {
          // add split hits to splitSelectedRecHits
          hitSplitter.split(*selectedRecHits[index], hitsForTrackCandidate, hitsAlongMomentum);
        } else {
          hitsForTrackCandidate.push_back(selectedRecHits[index]->clone());
        }
      }
 
      // set the recHitCombinationIndex
      fastTrackingUtilities::setRecHitCombinationIndex(hitsForTrackCandidate, icomb);
 
      // create track candidate state
      //   1. get seed state (defined on the surface of the most outer hit)
      DetId seedDetId(seed.startingState().detId());
      const GeomDet* gdet = trackerGeometry->idToDet(seedDetId);
      TrajectoryStateOnSurface seedTSOS =
          trajectoryStateTransform::transientState(seed.startingState(), &(gdet->surface()), magneticField.product());
      //   2. backPropagate the seedState to the surfuce of the most inner hit
      const GeomDet* initialLayer = trackerGeometry->idToDet(hitsForTrackCandidate.front().geographicalId());
      const TrajectoryStateOnSurface initialTSOS = propagator->propagate(seedTSOS, initialLayer->surface());
      //   3. check validity and transform
      if (!initialTSOS.isValid())
        continue;
      PTrajectoryStateOnDet PTSOD = trajectoryStateTransform::persistentState(
          initialTSOS, hitsForTrackCandidate.front().geographicalId().rawId());
 
      // add track candidate to output collection
      output->push_back(
          TrackCandidate(hitsForTrackCandidate, seed, PTSOD, edm::RefToBase<TrajectorySeed>(seeds, seedIndex)));
    }
  }
 
  // Save the track candidates
  e.put(std::move(output));
}
 
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(TrackCandidateProducer);