FinalTrackSelectors/plugins/TrackMerger.cc

0001 #include "TrackMerger.h"
0002
0003 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0004 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0005
0006 #define TRACK_SORT 1  // just use all hits from inner track, then append from outer outside it
0007 #define DET_SORT 0    // sort hits using global position of detector
0008 #define HIT_SORT 0    // sort hits using global position of hit
0009
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 // #define VI_DEBUG
0013
0014 #ifdef VI_DEBUG
0015 #define DPRINT(x) std::cout << x << ": "
0016 #define PRINT std::cout
0017 #else
0018 #define DPRINT(x) LogTrace(x)
0019 #define PRINT LogTrace("")
0020 #endif
0021
0022 TrackMerger::TrackMerger(const edm::ParameterSet &iConfig, edm::ConsumesCollector cc)
0023     : useInnermostState_(iConfig.getParameter<bool>("useInnermostState")),
0024       debug_(false),
0025       theBuilderName(iConfig.getParameter<std::string>("ttrhBuilderName")),
0026       geometryToken_(cc.esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0027       magFieldToken_(cc.esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0028       builderToken_(
0029           cc.esConsumes<TransientTrackingRecHitBuilder, TransientRecHitRecord>(edm::ESInputTag("", theBuilderName))),
0030       trackerTopoToken_(cc.esConsumes<TrackerTopology, TrackerTopologyRcd>()) {}
0031
0032 TrackMerger::~TrackMerger() {}
0033
0034 void TrackMerger::init(const edm::EventSetup &iSetup) {
0035   theGeometry = iSetup.getHandle(geometryToken_);
0036   theMagField = iSetup.getHandle(magFieldToken_);
0037   theBuilder = iSetup.getHandle(builderToken_);
0038   theTrkTopo = iSetup.getHandle(trackerTopoToken_);
0039 }
0040
0041 TrackCandidate TrackMerger::merge(const reco::Track &inner,
0042                                   const reco::Track &outer,
0043                                   DuplicateTrackType duplicateType) const {
0044   DPRINT("TrackMerger") << std::abs(inner.eta()) << " merging " << inner.algo() << '/' << outer.algo() << ' '
0045                         << inner.eta() << '/' << outer.eta() << std::endl;
0046
0047   std::vector<const TrackingRecHit *> hits;
0048   hits.reserve(inner.recHitsSize() + outer.recHitsSize());
0049   DPRINT("TrackMerger") << "Inner track hits: " << std::endl;
0050   for (auto it = inner.recHitsBegin(), ed = inner.recHitsEnd(); it != ed; ++it) {
0051     hits.push_back(&**it);
0052     if (debug_) {
0053       [[clang::suppress]]
0054       DetId id(hits.back()->geographicalId());
0055       PRINT << "   subdet " << id.subdetId() << "  layer " << theTrkTopo->layer(id) << " valid "
0056             << hits.back()->isValid() << "   detid: " << id() << std::endl;
0057     }
0058   }
0059   DPRINT("TrackMerger") << "Outer track hits: " << std::endl;
0060
0061   if (duplicateType == DuplicateTrackType::Disjoint) {
0062 #if TRACK_SORT
0063     DetId lastId(hits.back()->geographicalId());
0064     int lastSubdet = lastId.subdetId();
0065     unsigned int lastLayer = theTrkTopo->layer(lastId);
0066     for (auto it = outer.recHitsBegin(), ed = outer.recHitsEnd(); it != ed; ++it) {
0067       const TrackingRecHit *hit = &**it;
0068       DetId id(hit->geographicalId());
0069       int thisSubdet = id.subdetId();
0070       if (thisSubdet > lastSubdet || (thisSubdet == lastSubdet && theTrkTopo->layer(id) > lastLayer)) {
0071         hits.push_back(hit);
0072         PRINT << "  adding   subdet " << id.subdetId() << "  layer " << theTrkTopo->layer(id) << " valid "
0073               << hit->isValid() << "   detid: " << id() << std::endl;
0074       } else {
0075         PRINT << "  skipping subdet " << thisSubdet << "  layer " << theTrkTopo->layer(id) << " valid "
0076               << hit->isValid() << "   detid: " << id() << std::endl;
0077       }
0078     }
0079 #else
0080     addSecondTrackHits(hits, outer);
0081 #endif
0082   } else if (duplicateType == DuplicateTrackType::Overlapping) {
0083     addSecondTrackHits(hits, outer);
0084   }
0085
0086   math::XYZVector p = (inner.innerMomentum() + outer.outerMomentum());
0087   GlobalVector v(p.x(), p.y(), p.z());
0088   if (!useInnermostState_)
0089     v *= -1;
0090
0091   TrackCandidate::RecHitContainer ownHits;
0092   unsigned int nhits = hits.size();
0093   ownHits.reserve(nhits);
0094
0095   if (duplicateType == DuplicateTrackType::Disjoint) {
0096 #if TRACK_SORT
0097     if (!useInnermostState_)
0098       std::reverse(hits.begin(), hits.end());
0099     for (auto hit : hits)
0100       ownHits.push_back(*hit);
0101 #elif DET_SORT
0102     // OLD METHOD, sometimes fails
0103     std::sort(hits.begin(), hits.end(), MomentumSort(v, &*theGeometry));
0104     for (auto hit : hits)
0105       ownHits.push_back(*hit);
0106 #else
0107     sortByHitPosition(v, hits, ownHits);
0108 #endif
0109   } else if (duplicateType == DuplicateTrackType::Overlapping) {
0110     sortByHitPosition(v, hits, ownHits);
0111   }
0112
0113   PTrajectoryStateOnDet state;
0114   PropagationDirection pdir;
0115   if (useInnermostState_) {
0116     pdir = alongMomentum;
0117     if ((inner.outerPosition() - inner.innerPosition()).Dot(inner.momentum()) >= 0) {
0118       // use inner state
0119       TrajectoryStateOnSurface originalTsosIn(
0120           trajectoryStateTransform::innerStateOnSurface(inner, *theGeometry, &*theMagField));
0121       state = trajectoryStateTransform::persistentState(originalTsosIn, DetId(inner.innerDetId()));
0122     } else {
0123       // use outer state
0124       TrajectoryStateOnSurface originalTsosOut(
0125           trajectoryStateTransform::outerStateOnSurface(inner, *theGeometry, &*theMagField));
0126       state = trajectoryStateTransform::persistentState(originalTsosOut, DetId(inner.outerDetId()));
0127     }
0128   } else {
0129     pdir = oppositeToMomentum;
0130     if ((outer.outerPosition() - inner.innerPosition()).Dot(inner.momentum()) >= 0) {
0131       // use outer state
0132       TrajectoryStateOnSurface originalTsosOut(
0133           trajectoryStateTransform::outerStateOnSurface(outer, *theGeometry, &*theMagField));
0134       state = trajectoryStateTransform::persistentState(originalTsosOut, DetId(outer.outerDetId()));
0135     } else {
0136       // use inner state
0137       TrajectoryStateOnSurface originalTsosIn(
0138           trajectoryStateTransform::innerStateOnSurface(outer, *theGeometry, &*theMagField));
0139       state = trajectoryStateTransform::persistentState(originalTsosIn, DetId(outer.innerDetId()));
0140     }
0141   }
0142   TrajectorySeed seed(state, TrackCandidate::RecHitContainer(), pdir);
0143   TrackCandidate ret(ownHits, seed, state, (useInnermostState_ ? inner : outer).seedRef());
0144   ret.setStopReason((uint8_t)(useInnermostState_ ? inner : outer).stopReason());
0145   return ret;
0146 }
0147
0148 void TrackMerger::addSecondTrackHits(std::vector<const TrackingRecHit *> &hits, const reco::Track &outer) const {
0149   size_t nHitsFirstTrack = hits.size();
0150   for (auto it = outer.recHitsBegin(), ed = outer.recHitsEnd(); it != ed; ++it) {
0151     const TrackingRecHit *hit = &**it;
0152     DetId id(hit->geographicalId());
0153     const auto lay = theTrkTopo->layer(id);
0154     bool shared = false;
0155     bool valid = hit->isValid();
0156     PRINT << "   subdet " << id.subdetId() << "  layer " << theTrkTopo->layer(id) << " valid " << valid
0157           << "   detid: " << id() << std::endl;
0158     size_t iHit = 0;
0159     for (auto hit2 : hits) {
0160       ++iHit;
0161       if (iHit > nHitsFirstTrack)
0162         break;
0163       DetId id2 = hit2->geographicalId();
0164       if (id.subdetId() != id2.subdetId())
0165         continue;
0166       if (theTrkTopo->layer(id2) != lay)
0167         continue;
0168       if (hit->sharesInput(hit2, TrackingRecHit::all)) {
0169         PRINT << "        discared as duplicate of other hit" << id() << std::endl;
0170         shared = true;
0171         break;
0172       }
0173       if (hit2->isValid() && !valid) {
0174         PRINT << "        replacing old invalid hit on detid " << id2() << std::endl;
0175         hit = hit2;
0176         shared = true;
0177         break;
0178       }
0179       PRINT << "        discared as additional hit on layer that already contains hit with detid " << id() << std::endl;
0180       shared = true;
0181       break;
0182     }
0183     if (shared)
0184       continue;
0185     hits.push_back(hit);
0186   }
0187 }
0188
0189 void TrackMerger::sortByHitPosition(const GlobalVector &v,
0190                                     const std::vector<const TrackingRecHit *> &hits,
0191                                     TrackCandidate::RecHitContainer &ownHits) const {
0192   // NEW sort, more accurate
0193   unsigned int nhits = hits.size();
0194   std::vector<TransientTrackingRecHit::RecHitPointer> ttrh(nhits);
0195   for (unsigned int i = 0; i < nhits; ++i)
0196     ttrh[i] = theBuilder->build(hits[i]);
0197   std::sort(ttrh.begin(), ttrh.end(), GlobalMomentumSort(v));
0198   for (const auto &hit : ttrh)
0199     ownHits.push_back(*hit);
0200 }
0201
0202 bool TrackMerger::MomentumSort::operator()(const TrackingRecHit *hit1, const TrackingRecHit *hit2) const {
0203   const GeomDet *det1 = geom_->idToDet(hit1->geographicalId());
0204   const GeomDet *det2 = geom_->idToDet(hit2->geographicalId());
0205   GlobalPoint p1 = det1->toGlobal(LocalPoint(0, 0, 0));
0206   GlobalPoint p2 = det2->toGlobal(LocalPoint(0, 0, 0));
0207   return (p2 - p1).dot(dir_) > 0;
0208 }
0209
0210 bool TrackMerger::GlobalMomentumSort::operator()(const TransientTrackingRecHit::RecHitPointer &hit1,
0211                                                  const TransientTrackingRecHit::RecHitPointer &hit2) const {
0212   GlobalPoint p1 = hit1->isValid() ? hit1->globalPosition() : hit1->det()->position();
0213   GlobalPoint p2 = hit2->isValid() ? hit2->globalPosition() : hit2->det()->position();
0214   return (p2 - p1).dot(dir_) > 0;
0215 }