PatternTools/src/Trajectory.cc

0001 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0002 #include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h"
0003 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0004 #include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h"
0005 #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"
0006
0007 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include <algorithm>
0012
0013 namespace {
0014   template <typename DataContainer>
0015   unsigned short countTrailingValidHits(DataContainer const& meas) {
0016     unsigned short n = 0;
0017     for (auto it = meas.rbegin(); it != meas.rend(); ++it) {
0018       if (Trajectory::lost(*(*it).recHit()))
0019         break;
0020       if ((*it).recHit()->isValid())
0021         ++n;
0022     }
0023     return n;
0024   }
0025
0026 }  // namespace
0027
0028 using namespace std;
0029
0030 void Trajectory::pop() {
0031   if (!empty()) {
0032     if (theData.back().recHit()->isValid()) {
0033       theNumberOfFoundHits--;
0034       theChiSquared -= theData.back().estimate();
0035       if (badForCCC(theData.back()))
0036         theNumberOfCCCBadHits_--;
0037       if (pixel(*(theData.back().recHit())))
0038         theNumberOfFoundPixelHits--;
0039     } else if (lost(*(theData.back().recHit()))) {
0040       theNumberOfLostHits--;
0041     } else if (isBad(*(theData.back().recHit())) && theData.back().recHit()->geographicalId().det() == DetId::Muon) {
0042       theChiSquaredBad -= theData.back().estimate();
0043     }
0044
0045     theData.pop_back();
0046     theNumberOfTrailingFoundHits = countTrailingValidHits(theData);
0047   }
0048 }
0049
0050 void Trajectory::push(const TrajectoryMeasurement& tm) { push(tm, tm.estimate()); }
0051
0052 void Trajectory::push(TrajectoryMeasurement&& tm) { push(tm, tm.estimate()); }
0053
0054 void Trajectory::push(const TrajectoryMeasurement& tm, double chi2Increment) {
0055   theData.push_back(tm);
0056   pushAux(chi2Increment);
0057 }
0058
0059 void Trajectory::push(TrajectoryMeasurement&& tm, double chi2Increment) {
0060   theData.push_back(tm);
0061   pushAux(chi2Increment);
0062 }
0063
0064 void Trajectory::pushAux(double chi2Increment) {
0065   const TrajectoryMeasurement& tm = theData.back();
0066   if (tm.recHit()->isValid()) {
0067     theChiSquared += chi2Increment;
0068     theNumberOfFoundHits++;
0069     theNumberOfTrailingFoundHits++;
0070     if (badForCCC(tm))
0071       theNumberOfCCCBadHits_++;
0072     if (pixel(*(tm.recHit())))
0073       theNumberOfFoundPixelHits++;
0074   }
0075   // else if (lost( tm.recHit()) && !inactive(tm.recHit().det())) theNumberOfLostHits++;
0076   else if (lost(*(tm.recHit()))) {
0077     theNumberOfLostHits++;
0078     theNumberOfTrailingFoundHits = 0;
0079   }
0080
0081   else if (isBad(*(tm.recHit())) && tm.recHit()->geographicalId().det() == DetId::Muon) {
0082     theChiSquaredBad += chi2Increment;
0083   }
0084
0085   // in case of a Trajectory constructed without direction,
0086   // determine direction from the radii of the first two measurements
0087
0088   if (!theDirectionValidity && theData.size() >= 2) {
0089     if (theData[0].updatedState().globalPosition().perp2() < theData.back().updatedState().globalPosition().perp2())
0090       theDirection = alongMomentum;
0091     else
0092       theDirection = oppositeToMomentum;
0093     theDirectionValidity = true;
0094   }
0095 }
0096
0097 int Trajectory::ndof(bool bon) const {
0098   Trajectory::RecHitContainer&& transRecHits = recHits();
0099
0100   int dof = 0;
0101   int dofBad = 0;
0102
0103   for (auto& rechit : transRecHits) {
0104     if ((rechit)->isValid())
0105       dof += (rechit)->dimension();
0106     else if (isBad(*rechit) && (rechit)->geographicalId().det() == DetId::Muon)
0107       dofBad += (rechit)->dimension();
0108   }
0109
0110   // If dof!=0 (there is at least 1 valid hit),
0111   //    return ndof=ndof(fit)
0112   // If dof=0 (all rec hits are invalid, only for STA trajectories),
0113   //    return ndof=ndof(invalid hits)
0114   if (dof) {
0115     int constr = bon ? 5 : 4;
0116     return std::max(dof - constr, 0);
0117   } else {
0118     // A STA can have < 5 (invalid) hits
0119     // if this is the case ==> ndof = 1
0120     // (to avoid divisions by 0)
0121     int constr = bon ? 5 : 4;
0122     return std::max(dofBad - constr, 1);
0123   }
0124 }
0125
0126 void Trajectory::validRecHits(ConstRecHitContainer& hits) const {
0127   hits.reserve(foundHits());
0128   for (auto const& tm : theData)
0129     if (tm.recHit()->isValid())
0130       hits.push_back(tm.recHit());
0131 }
0132
0133 PropagationDirection const& Trajectory::direction() const {
0134   if (theDirectionValidity)
0135     return theDirection;
0136   else
0137     throw cms::Exception("TrackingTools/PatternTools", "Trajectory::direction() requested but not set");
0138 }
0139
0140 void Trajectory::check() const {
0141   if (theData.empty())
0142     throw cms::Exception("TrackingTools/PatternTools",
0143                          "Trajectory::check() - information requested from empty Trajectory");
0144 }
0145
0146 bool Trajectory::lost(const TrackingRecHit& hit) {
0147   if (hit.isValid())
0148     return false;
0149   else {
0150     //     // A DetLayer is always inactive in this logic.
0151     //     // The DetLayer is the Det of an invalid RecHit only if no DetUnit
0152     //     // is compatible with the predicted state, so we don't really expect
0153     //     // a hit in this case.
0154
0155     if (hit.geographicalId().rawId() == 0) {
0156       return false;
0157     } else {
0158       return hit.getType() == TrackingRecHit::missing;
0159     }
0160   }
0161 }
0162
0163 bool Trajectory::isBad(const TrackingRecHit& hit) {
0164   if (hit.isValid())
0165     return false;
0166   else {
0167     if (hit.geographicalId().rawId() == 0) {
0168       return false;
0169     } else {
0170       return hit.getType() == TrackingRecHit::bad;
0171     }
0172   }
0173 }
0174
0175 bool Trajectory::pixel(const TrackingRecHit& hit) {
0176   if (!trackerHitRTTI::isFromDetOrFast(hit))
0177     return false;
0178   auto const* thit = static_cast<const BaseTrackerRecHit*>(hit.hit());
0179   return thit->isPixel();
0180 }
0181
0182 bool Trajectory::badForCCC(const TrajectoryMeasurement& tm) {
0183   if (!trackerHitRTTI::isFromDet(*tm.recHit()))
0184     return false;
0185   auto const* thit = static_cast<const BaseTrackerRecHit*>(tm.recHit()->hit());
0186   if (!thit)
0187     return false;
0188   if (thit->isPixel() || thit->isPhase2())
0189     return false;
0190   if (!tm.updatedState().isValid())
0191     return false;
0192   return siStripClusterTools::chargePerCM(thit->rawId(),
0193                                           thit->firstClusterRef().stripCluster(),
0194                                           tm.updatedState().localParameters()) < theCCCThreshold_;
0195 }
0196
0197 void Trajectory::updateBadForCCC(float ccc_threshold) {
0198   // If the supplied threshold is the same as the currently cached
0199   // one, then return the current number of bad hits for CCC,
0200   // otherwise do a new full rescan.
0201   if (ccc_threshold == theCCCThreshold_)
0202     return;
0203
0204   theCCCThreshold_ = ccc_threshold;
0205   theNumberOfCCCBadHits_ = 0;
0206   for (auto const& h : theData) {
0207     if (badForCCC(h))
0208       theNumberOfCCCBadHits_++;
0209   }
0210 }
0211
0212 int Trajectory::numberOfCCCBadHits(float ccc_threshold) {
0213   updateBadForCCC(ccc_threshold);
0214   return theNumberOfCCCBadHits_;
0215 }
0216
0217 TrajectoryStateOnSurface Trajectory::geometricalInnermostState() const {
0218   check();
0219
0220   //if trajectory is in one half, return the end closer to origin point
0221   if (firstMeasurement().updatedState().globalMomentum().perp() > 1.0 &&
0222       (firstMeasurement().updatedState().globalPosition().basicVector().dot(
0223            firstMeasurement().updatedState().globalMomentum().basicVector()) *
0224            lastMeasurement().updatedState().globalPosition().basicVector().dot(
0225                lastMeasurement().updatedState().globalMomentum().basicVector()) >
0226        0)) {
0227     return (firstMeasurement().updatedState().globalPosition().mag() <
0228             lastMeasurement().updatedState().globalPosition().mag())
0229                ? firstMeasurement().updatedState()
0230                : lastMeasurement().updatedState();
0231   }
0232
0233   //more complicated in case of traversing and low-pt trajectories with loops
0234   return closestMeasurement(GlobalPoint(0.0, 0.0, 0.0)).updatedState();
0235 }
0236
0237 namespace {
0238   /// used to determine closest measurement to given point
0239   struct LessMag {
0240     LessMag(GlobalPoint point) : thePoint(point) {}
0241     bool operator()(const TrajectoryMeasurement& lhs, const TrajectoryMeasurement& rhs) const {
0242       if (lhs.updatedState().isValid() && rhs.updatedState().isValid())
0243         return (lhs.updatedState().globalPosition() - thePoint).mag2() <
0244                (rhs.updatedState().globalPosition() - thePoint).mag2();
0245       else {
0246         edm::LogError("InvalidStateOnMeasurement")
0247             << "an updated state is not valid. result of LessMag comparator will be wrong.";
0248         return false;
0249       }
0250     }
0251     GlobalPoint thePoint;
0252   };
0253
0254 }  // namespace
0255
0256 TrajectoryMeasurement const& Trajectory::closestMeasurement(GlobalPoint point) const {
0257   check();
0258   vector<TrajectoryMeasurement>::const_iterator iter =
0259       std::min_element(measurements().begin(), measurements().end(), LessMag(point));
0260
0261   return (*iter);
0262 }
0263
0264 void Trajectory::reverse() {
0265   // reverse the direction (without changing it if it's not along or opposite)
0266   if (theDirection == alongMomentum)
0267     theDirection = oppositeToMomentum;
0268   else if (theDirection == oppositeToMomentum)
0269     theDirection = alongMomentum;
0270   // reverse the order of the hits
0271   std::reverse(theData.begin(), theData.end());
0272 }