PatternTools/src/TempTrajectory.cc

0001 #include "TrackingTools/PatternTools/interface/TempTrajectory.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 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0007 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0010 #include "FWCore/Utilities/interface/Likely.h"
0011
0012 namespace {
0013   template <typename DataContainer>
0014   unsigned short countTrailingValidHits(DataContainer const& meas) {
0015     unsigned short n = 0;
0016     for (auto it = meas.rbegin(); it != meas.rend(); --it) {  // it is not consistent with std...
0017       if (TempTrajectory::lost(*(*it).recHit()))
0018         break;
0019       if ((*it).recHit()->isValid())
0020         ++n;
0021     }
0022     return n;
0023   }
0024 }  // namespace
0025
0026 TempTrajectory::TempTrajectory(Trajectory&& traj) : thePayload(std::make_unique<Payload>()) {
0027   assert(traj.isValid());
0028   thePayload->theDirection = traj.direction();
0029   thePayload->theNHseed = traj.seedNHits();
0030   thePayload->theNLoops = traj.nLoops();
0031   thePayload->theDPhiCache = traj.dPhiCacheForLoopersReconstruction();
0032   thePayload->theCCCThreshold_ = traj.cccThreshold();
0033   thePayload->stopReason_ = traj.stopReason();
0034   for (auto& it : traj.measurements()) {
0035     push(it);
0036   }
0037   traj.measurements().clear();
0038 }
0039
0040 void TempTrajectory::pop() {
0041   if (!empty()) {
0042     if (theData.back().recHit()->isValid()) {
0043       thePayload->theNumberOfFoundHits--;
0044       if (badForCCC(theData.back()))
0045         thePayload->theNumberOfCCCBadHits_--;
0046       if (Trajectory::pixel(*(theData.back().recHit())))
0047         thePayload->theNumberOfFoundPixelHits--;
0048     } else if (lost(*(theData.back().recHit()))) {
0049       thePayload->theNumberOfLostHits--;
0050     }
0051     theData.pop_back();
0052     thePayload->theNumberOfTrailingFoundHits = countTrailingValidHits(theData);
0053   }
0054 }
0055
0056 void TempTrajectory::pushAux(double chi2Increment) {
0057   const TrajectoryMeasurement& tm = theData.back();
0058   if (tm.recHit()->isValid()) {
0059     thePayload->theNumberOfFoundHits++;
0060     thePayload->theNumberOfTrailingFoundHits++;
0061     if (badForCCC(tm))
0062       thePayload->theNumberOfCCCBadHits_++;
0063     if (Trajectory::pixel(*(tm.recHit())))
0064       thePayload->theNumberOfFoundPixelHits++;
0065   }
0066   //else if (lost( tm.recHit()) && !inactive(tm.recHit().det())) theNumberOfLostHits++;
0067   else if (lost(*(tm.recHit()))) {
0068     thePayload->theNumberOfLostHits++;
0069     thePayload->theNumberOfTrailingFoundHits = 0;
0070   }
0071
0072   thePayload->theChiSquared += chi2Increment;
0073 }
0074
0075 void TempTrajectory::push(const TempTrajectory& segment) {
0076   assert(segment.thePayload->theDirection == thePayload->theDirection);
0077   assert(segment.thePayload->theCCCThreshold_ == thePayload->theCCCThreshold_);
0078
0079   const int N = segment.measurements().size();
0080   TrajectoryMeasurement const* tmp[N];
0081   int i = 0;
0082   //for (DataContainer::const_iterator it = segment.measurements().rbegin(), ed = segment.measurements().rend(); it != ed; --it)
0083   for (auto const& tm : segment.measurements())
0084     tmp[i++] = &tm;
0085   while (i != 0)
0086     theData.push_back(*tmp[--i]);
0087   thePayload->theNumberOfFoundHits += segment.thePayload->theNumberOfFoundHits;
0088   thePayload->theNumberOfFoundPixelHits += segment.thePayload->theNumberOfFoundPixelHits;
0089   thePayload->theNumberOfLostHits += segment.thePayload->theNumberOfLostHits;
0090   thePayload->theNumberOfCCCBadHits_ += segment.thePayload->theNumberOfCCCBadHits_;
0091   thePayload->theNumberOfTrailingFoundHits = countTrailingValidHits(theData);
0092   thePayload->theChiSquared += segment.thePayload->theChiSquared;
0093 }
0094
0095 void TempTrajectory::join(TempTrajectory& segment) {
0096   assert(segment.thePayload->theDirection == thePayload->theDirection);
0097
0098   if (thePayload->theCCCThreshold_ != segment.thePayload->theCCCThreshold_)
0099     segment.updateBadForCCC(thePayload->theCCCThreshold_);
0100   if (segment.theData.shared()) {
0101     push(segment);
0102     segment.theData.clear();  // obey the contract, and increase the chances it will be not shared one day
0103   } else {
0104     theData.join(segment.theData);
0105     thePayload->theNumberOfFoundHits += segment.thePayload->theNumberOfFoundHits;
0106     thePayload->theNumberOfFoundPixelHits += segment.thePayload->theNumberOfFoundPixelHits;
0107     thePayload->theNumberOfLostHits += segment.thePayload->theNumberOfLostHits;
0108     thePayload->theNumberOfCCCBadHits_ += segment.thePayload->theNumberOfCCCBadHits_;
0109     thePayload->theNumberOfTrailingFoundHits = countTrailingValidHits(theData);
0110     thePayload->theChiSquared += segment.thePayload->theChiSquared;
0111   }
0112 }
0113
0114 PropagationDirection TempTrajectory::direction() const { return PropagationDirection(thePayload->theDirection); }
0115
0116 void TempTrajectory::check() const {
0117   if (theData.empty())
0118     throw cms::Exception("TrackingTools/PatternTools",
0119                          "Trajectory::check() - information requested from empty Trajectory");
0120 }
0121
0122 bool TempTrajectory::lost(const TrackingRecHit& hit) {
0123   if LIKELY (hit.isValid())
0124     return false;
0125
0126   //     // A DetLayer is always inactive in this logic.
0127   //     // The DetLayer is the Det of an invalid RecHit only if no DetUnit
0128   //     // is compatible with the predicted state, so we don't really expect
0129   //     // a hit in this case.
0130
0131   if (hit.geographicalId().rawId() == 0) {
0132     return false;
0133   }
0134   return hit.getType() == TrackingRecHit::missing;
0135 }
0136
0137 bool TempTrajectory::badForCCC(const TrajectoryMeasurement& tm) {
0138   if (!trackerHitRTTI::isFromDet(*tm.recHit()))
0139     return false;
0140   auto const* thit = static_cast<const BaseTrackerRecHit*>(tm.recHit()->hit());
0141   if (!thit)
0142     return false;
0143   if (thit->isPixel() || thit->isPhase2())
0144     return false;
0145   if (!tm.updatedState().isValid())
0146     return false;
0147   return siStripClusterTools::chargePerCM(thit->rawId(),
0148                                           thit->firstClusterRef().stripCluster(),
0149                                           tm.updatedState().localParameters()) < thePayload->theCCCThreshold_;
0150 }
0151
0152 void TempTrajectory::updateBadForCCC(float ccc_threshold) {
0153   // If the supplied threshold is the same as the currently cached
0154   // one, then return the current number of bad hits for CCC,
0155   // otherwise do a new full rescan.
0156   if (ccc_threshold == thePayload->theCCCThreshold_)
0157     return;
0158
0159   thePayload->theCCCThreshold_ = ccc_threshold;
0160   thePayload->theNumberOfCCCBadHits_ = 0;
0161   for (auto const& h : theData) {
0162     if (badForCCC(h))
0163       thePayload->theNumberOfCCCBadHits_++;
0164   }
0165 }
0166
0167 int TempTrajectory::numberOfCCCBadHits(float ccc_threshold) {
0168   updateBadForCCC(ccc_threshold);
0169   return thePayload->theNumberOfCCCBadHits_;
0170 }
0171
0172 Trajectory TempTrajectory::toTrajectory() const {
0173   assert(isValid());
0174   PropagationDirection p = PropagationDirection(thePayload->theDirection);
0175   Trajectory traj(p);
0176   traj.setNLoops(thePayload->theNLoops);
0177   traj.setStopReason(thePayload->stopReason_);
0178   traj.numberOfCCCBadHits(thePayload->theCCCThreshold_);
0179
0180   traj.reserve(theData.size());
0181   const TrajectoryMeasurement* tmp[theData.size()];
0182   int i = 0;
0183   for (DataContainer::const_iterator it = theData.rbegin(), ed = theData.rend(); it != ed; --it)
0184     tmp[i++] = &(*it);
0185   while (i != 0)
0186     traj.push(*tmp[--i]);
0187   return traj;
0188 }