File indexing completed on 2024-04-06 12:31:34
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 }
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
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
0086
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
0111
0112
0113
0114 if (dof) {
0115 int constr = bon ? 5 : 4;
0116 return std::max(dof - constr, 0);
0117 } else {
0118
0119
0120
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
0151
0152
0153
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
0199
0200
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
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
0234 return closestMeasurement(GlobalPoint(0.0, 0.0, 0.0)).updatedState();
0235 }
0236
0237 namespace {
0238
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 }
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
0266 if (theDirection == alongMomentum)
0267 theDirection = oppositeToMomentum;
0268 else if (theDirection == oppositeToMomentum)
0269 theDirection = alongMomentum;
0270
0271 std::reverse(theData.begin(), theData.end());
0272 }