PatternTools/interface/Trajectory.h

0001 #ifndef CommonDet_Trajectory_H
0002 #define CommonDet_Trajectory_H
0003
0004 #include "DataFormats/Common/interface/RefToBase.h"
0005 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0006 #include "DataFormats/TrackCandidate/interface/TrajectoryStopReasons.h"
0007 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0008 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0009 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include <vector>
0013 #include <algorithm>
0014 #include <limits>
0015
0016 /** A class for detailed particle trajectory representation.
0017  *  It is used during trajectory building to "grow" a trajectory.
0018  *  The trajectory is represented as an ordered sequence of
0019  *  TrajectoryMeasurement objects with a stack-like interface.
0020  *  The measurements are added to the Trajectory in the order of
0021  *  increasing precision: each new TrajectoryMeasurement is assumed to improve
0022  *  the precision of the last one, normally by adding a constraint from
0023  *  a new RecHit.
0024  *     However the Trajectory class does not have the means to verify
0025  *  that measurements are added in the correct order, and thus cannot
0026  *  guarantee the order, which is the responsibility of the
0027  *  TrajectoryBuilder. The Trajectory provides some security by
0028  *  allowing to add or remove measurements only on one of it's ends,
0029  *  with push(TM) and pop() methods. The last measurement in a Trajectory
0030  *  can thus be either the innermost (closest to the interaction point)
0031  *  or the outermost, depending on the way the Trajectory was built.
0032  *  The direction of building is represented as a PropagationDirection,
0033  *  which has two possible values: alongMomentum (outwards) and
0034  *  oppositeToMomentum (inwards), and is accessed with the direction()
0035  *  method.
0036  */
0037
0038 class Trajectory {
0039 public:
0040   typedef std::vector<TrajectoryMeasurement> DataContainer;
0041   using ConstRecHitContainer = TrackingRecHit::ConstRecHitContainer;
0042   using RecHitContainer = ConstRecHitContainer;
0043
0044   /** Default constructor of an empty trajectory with undefined seed and
0045    * undefined direction. This constructor is necessary in order to transiently
0046    * copy vector<Trajectory> in the edm::Event
0047    */
0048
0049   Trajectory() = default;
0050
0051   /** Constructor of an empty trajectory with undefined direction.
0052    *  The direction will be defined at the moment of the push of a second
0053    *  measurement, from the relative radii of the first and second
0054    *  measurements.
0055    */
0056
0057   explicit Trajectory(const TrajectorySeed& seed) : theSeed(new TrajectorySeed(seed)), theValid(true) {}
0058
0059   /** Constructor of an empty trajectory with defined direction.
0060    *  No check is made in the push method that measurements are
0061    *  added in the correct direction.
0062    */
0063   Trajectory(const TrajectorySeed& seed, PropagationDirection dir)
0064       : theSeed(new TrajectorySeed(seed)), theDirection(dir), theDirectionValidity(true), theValid(true) {}
0065
0066   /** Constructor of an empty trajectory with defined direction.
0067    *  No check is made in the push method that measurements are
0068    *  added in the correct direction.
0069    */
0070   Trajectory(const std::shared_ptr<const TrajectorySeed>& seed, PropagationDirection dir)
0071       : theSeed(seed), theDirection(dir), theDirectionValidity(true), theValid(true) {}
0072
0073   /** Constructor of an empty trajectory with defined direction.
0074    *  No check is made in the push method that measurements are
0075    *  added in the correct direction.
0076    */
0077   explicit Trajectory(PropagationDirection dir) : theDirection(dir), theDirectionValidity(true), theValid(true) {}
0078
0079   Trajectory(Trajectory const& rh) = default;
0080   Trajectory& operator=(Trajectory const& rh) = default;
0081
0082   Trajectory(Trajectory&& rh)
0083       : theSeed(std::move(rh.theSeed)),
0084         seedRef_(std::move(rh.seedRef_)),
0085         theData(std::move(rh.theData)),
0086         theDPhiCache(rh.theDPhiCache),
0087         theCCCThreshold_(rh.theCCCThreshold_),
0088         theChiSquared(rh.theChiSquared),
0089         theChiSquaredBad(rh.theChiSquaredBad),
0090         theDirection(rh.theDirection),
0091         theDirectionValidity(rh.theDirectionValidity),
0092         theValid(rh.theValid),
0093         theNumberOfFoundHits(rh.theNumberOfFoundHits),
0094         theNumberOfFoundPixelHits(rh.theNumberOfFoundPixelHits),
0095         theNumberOfLostHits(rh.theNumberOfLostHits),
0096         theNumberOfTrailingFoundHits(rh.theNumberOfTrailingFoundHits),
0097         theNumberOfCCCBadHits_(rh.theNumberOfCCCBadHits_),
0098         theNLoops(rh.theNLoops),
0099         stopReason_(rh.stopReason_) {}
0100
0101   Trajectory& operator=(Trajectory&& rh) {
0102     using std::swap;
0103     swap(theData, rh.theData);
0104     theDPhiCache = rh.theDPhiCache;
0105     theCCCThreshold_ = rh.theCCCThreshold_;
0106     theChiSquared = rh.theChiSquared;
0107     theChiSquaredBad = rh.theChiSquaredBad;
0108     theValid = rh.theValid;
0109     theNLoops = rh.theNLoops;
0110     theNumberOfFoundHits = rh.theNumberOfFoundHits;
0111     theNumberOfFoundPixelHits = rh.theNumberOfFoundPixelHits;
0112     theNumberOfLostHits = rh.theNumberOfLostHits;
0113     theNumberOfTrailingFoundHits = rh.theNumberOfTrailingFoundHits;
0114     theNumberOfCCCBadHits_ = rh.theNumberOfCCCBadHits_;
0115     theDirection = rh.theDirection;
0116     theDirectionValidity = rh.theDirectionValidity;
0117     stopReason_ = rh.stopReason_;
0118     swap(theSeed, rh.theSeed);
0119     swap(seedRef_, rh.seedRef_);
0120
0121     return *this;
0122   }
0123
0124   /** Reserves space in the vector to avoid lots of allocations when
0125       push_back-ing measurements */
0126   void reserve(unsigned int n) { theData.reserve(n); }
0127
0128   /** Add a new measurement to a Trajectory.
0129    *  The Chi2 of the trajectory is incremented by the value
0130    *  of tm.estimate() .
0131    */
0132   void push(const TrajectoryMeasurement& tm);
0133   /** same as the one-argument push, but the trajectory Chi2 is incremented
0134    *  by chi2Increment. Useful e.g. in trajectory smoothing.
0135    */
0136   void push(const TrajectoryMeasurement& tm, double chi2Increment);
0137
0138   void push(TrajectoryMeasurement&& tm);
0139   void push(TrajectoryMeasurement&& tm, double chi2Increment);
0140
0141   /** Remove the last measurement from the trajectory.
0142    */
0143   void pop();
0144
0145   /** Access to the last measurement.
0146    *  It's the most precise one in a trajectory before smoothing.
0147    *  It's the outermost measurement if direction() == alongMomentum,
0148    *  the innermost one if direction() == oppositeToMomentum.
0149    */
0150   TrajectoryMeasurement const& lastMeasurement() const {
0151     check();
0152     if (theData.back().recHitR().hit() != nullptr)
0153       return theData.back();
0154     else if (theData.size() > 2)
0155       return *(theData.end() - 2);
0156     else
0157       throw cms::Exception("TrajectoryMeasurement::lastMeasurement - Too few measurements in trajectory");
0158   }
0159
0160   /** Access to the first measurement.
0161    *  It is the least precise one in a trajectory before smoothing.
0162    *  It is precise in a smoothed trajectory.
0163    *  It's the innermost measurement if direction() == alongMomentum,
0164    *  the outermost one if direction() == oppositeToMomentum.
0165    */
0166   TrajectoryMeasurement const& firstMeasurement() const {
0167     check();
0168     if (theData.front().recHitR().hit() != nullptr)
0169       return theData.front();
0170     else if (theData.size() > 2)
0171       return *(theData.begin() + 1);
0172     else
0173       throw cms::Exception("TrajectoryMeasurement::firstMeasurement - Too few measurements in trajectory");
0174   }
0175
0176   /** Return all measurements in a container.
0177    */
0178   DataContainer const& measurements() const { return theData; }
0179   DataContainer& measurements() { return theData; }
0180
0181   /// obsolete name, use measurements() instead.
0182   DataContainer const& data() const { return measurements(); }
0183
0184   /** Return all RecHits in a container.
0185    */
0186   ConstRecHitContainer recHits() const {
0187     ConstRecHitContainer hits;
0188     hits.reserve(theData.size());
0189     for (Trajectory::DataContainer::const_iterator itm = theData.begin(); itm != theData.end(); itm++) {
0190       hits.push_back((*itm).recHit());
0191     }
0192     return hits;
0193   }
0194
0195   /** Just valid hits..
0196    *
0197    */
0198   void validRecHits(ConstRecHitContainer& cont) const;
0199
0200   /** Number of valid RecHits used to determine the trajectory.
0201    *  Can be less than the number of measurements in data() since
0202    *  detector layers crossed without using RecHits from them are also
0203    *  stored as measurements.
0204    */
0205
0206   int foundHits() const { return theNumberOfFoundHits; }
0207
0208   /** Number of valid pixel RecHits used to determine the trajectory.
0209    */
0210   int foundPixelHits() const { return theNumberOfFoundPixelHits; }
0211
0212   /** Number of detector layers crossed without valid RecHits.
0213    *  Used mainly as a criteria for abandoning a trajectory candidate
0214    *  during trajectory building.
0215    */
0216
0217   int lostHits() const { return theNumberOfLostHits; }
0218
0219   /** Number of valid RecHits at the end of the trajectory after last lost hit.
0220    */
0221   int trailingFoundHits() const { return theNumberOfTrailingFoundHits; }
0222
0223   /** Number of hits that are not compatible with the CCC used during
0224    *  patter recognition. Used mainly as a criteria for abandoning a
0225    *  trajectory candidate during trajectory building.
0226    */
0227   int cccBadHits() const { return theNumberOfCCCBadHits_; }
0228
0229   //number of hits in seed
0230   unsigned int seedNHits() const { return seed().nHits(); }
0231
0232   /// True if trajectory has no measurements.
0233   bool empty() const { return theData.empty(); }
0234
0235   /// - Trajectories with at least 1 valid hit:
0236   ///     value of the raw Chi2 of the trajectory, not normalised to the N.D.F.
0237   ///     (evaluated using only valid hits)
0238   /// - Trajectories with only invalid hits:
0239   ///     raw Chi2 (not norm.) of invalid hits w.r.t. the "default" trajectory
0240   ///     (traj. containing only the seed information)
0241   float chiSquared() const { return (theNumberOfFoundHits ? theChiSquared : theChiSquaredBad); }
0242
0243   /// Number of dof of the trajectory. The method accepts a bool in order to properly
0244   /// take into account the presence of magnetic field in the dof computation.
0245   /// By default the MF is considered ON.
0246   int ndof(bool bon = true) const;
0247
0248   /** Direction of "growing" of the trajectory.
0249    *  Possible values are alongMomentum (outwards) and
0250    *  oppositeToMomentum (inwards).
0251    */
0252   PropagationDirection const& direction() const;
0253
0254   /** Returns true if the Trajectory is valid.
0255    *  Trajectories are invalidated e.g. during ambiguity resolution.
0256    */
0257   bool isValid() const { return theValid; }
0258
0259   /// Method to invalidate a trajectory. Useful during ambiguity resolution.
0260   void invalidate() { theValid = false; }
0261
0262   /// Access to the seed used to reconstruct the Trajectory
0263   TrajectorySeed const& seed() const { return *theSeed; }
0264
0265   /** Definition of inactive Det from the Trajectory point of view.
0266    */
0267   static bool inactive(  //const Det& det
0268   ) {
0269     return false;
0270   }  //FIXME
0271
0272   /** Definition of what it means for a hit to be "lost".
0273    *  This definition is also used by the TrajectoryBuilder.
0274    */
0275   static bool lost(const TrackingRecHit& hit);
0276
0277   /** Returns true if the hit type is TrackingRecHit::bad
0278    *  Used in stand-alone trajectory construction
0279    */
0280   static bool isBad(const TrackingRecHit& hit);
0281
0282   /** Returns true if the hit type is TrackingRecHit::bad
0283    *  Used in trajectory filtering
0284    */
0285   static bool pixel(const TrackingRecHit& hit);
0286
0287   /// Redundant method, returns the layer of lastMeasurement() .
0288   const DetLayer* lastLayer() const {
0289     check();
0290     if (theData.back().recHit()->hit() != nullptr)
0291       return theData.back().layer();
0292     else if (theData.size() > 2)
0293       return (theData.end() - 2)->layer();
0294     else
0295       throw cms::Exception("TrajectoryMeasurement::lastMeasurement - Too few measurements in trajectory");
0296   }
0297
0298   /**  return the Reference to the trajectory seed in the original
0299    *   seeds collection. If the collection has been dropped from the
0300    *   Event, the reference may be invalid. Its validity should be tested,
0301    *   before the reference is actually used.
0302    */
0303   edm::RefToBase<TrajectorySeed> seedRef(void) const { return seedRef_; }
0304
0305   void setSeedRef(const edm::RefToBase<TrajectorySeed>& seedRef) { seedRef_ = seedRef; }
0306
0307   TrajectoryStateOnSurface geometricalInnermostState() const;
0308
0309   TrajectoryMeasurement const& closestMeasurement(GlobalPoint) const;
0310
0311   /// Reverse the propagation direction and the order of the trajectory measurements.
0312   /// It doesn't reverse the forward and backward predicted states within each trajectory measurement
0313   void reverse();
0314
0315   const std::shared_ptr<const TrajectorySeed>& sharedSeed() const { return theSeed; }
0316   void setSharedSeed(const std::shared_ptr<const TrajectorySeed>& seed) { theSeed = seed; }
0317
0318   /// accessor to the delta phi angle betweem the directions of the two measurements on the last
0319   /// two layers crossed by the trajectory
0320   float dPhiCacheForLoopersReconstruction() const { return theDPhiCache; }
0321
0322   float cccThreshold() const { return theCCCThreshold_; }
0323
0324   /// method to set the delta phi angle betweem the directions of the two measurements on the last
0325   /// two layers crossed by the trajectory
0326   void setDPhiCacheForLoopersReconstruction(float dphi) { theDPhiCache = dphi; }
0327
0328   bool isLooper() const { return (theNLoops > 0); }
0329   int8_t nLoops() const { return theNLoops; }
0330
0331   void setNLoops(int8_t value) { theNLoops = value; }
0332   void incrementLoops() { theNLoops++; }
0333
0334   void setStopReason(StopReason s) { stopReason_ = s; }
0335   StopReason stopReason() const { return stopReason_; }
0336
0337   int numberOfCCCBadHits(float ccc_threshold);
0338
0339 private:
0340   void pushAux(double chi2Increment);
0341   bool badForCCC(const TrajectoryMeasurement& tm);
0342   void updateBadForCCC(float ccc_threshold);
0343
0344   std::shared_ptr<const TrajectorySeed> theSeed;
0345   edm::RefToBase<TrajectorySeed> seedRef_;
0346
0347   DataContainer theData;
0348
0349   float theDPhiCache = 0;
0350   float theCCCThreshold_ = std::numeric_limits<float>::max();
0351
0352   float theChiSquared = 0;
0353   float theChiSquaredBad = 0;
0354
0355   PropagationDirection theDirection = anyDirection;
0356   bool theDirectionValidity = false;
0357   bool theValid = false;
0358
0359   uint8_t theNumberOfFoundHits = 0;
0360   uint8_t theNumberOfFoundPixelHits = 0;
0361   uint8_t theNumberOfLostHits = 0;
0362   uint8_t theNumberOfTrailingFoundHits = 0;
0363   uint8_t theNumberOfCCCBadHits_ = 0;
0364   int8_t theNLoops = 0;
0365   StopReason stopReason_ = StopReason::UNINITIALIZED;
0366
0367   void check() const;
0368 };
0369
0370 #endif