File indexing completed on 2024-05-29 23:13:06
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryUpdator.h"
0018 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0019 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0020
0021 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0022 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0023 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0025 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0026 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0027 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0028
0029 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0030 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"
0031
0032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0033 #include "FWCore/Utilities/interface/Exception.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035
0036 #include <algorithm>
0037
0038 using namespace edm;
0039 using namespace std;
0040
0041
0042 MuonTrajectoryUpdator::MuonTrajectoryUpdator(const edm::ParameterSet& par, NavigationDirection fitDirection)
0043 : theFitDirection(fitDirection) {
0044
0045 theMaxChi2 = par.getParameter<double>("MaxChi2");
0046 theEstimator = new Chi2MeasurementEstimator(theMaxChi2);
0047
0048
0049 theUpdator = new KFUpdator();
0050
0051
0052 theGranularity = par.getParameter<int>("Granularity");
0053
0054
0055 theRescaleErrorFlag = par.getParameter<bool>("RescaleError");
0056
0057 if (theRescaleErrorFlag)
0058
0059 theRescaleFactor = par.getParameter<double>("RescaleErrorFactor");
0060
0061
0062 useInvalidHits = par.getParameter<bool>("UseInvalidHits");
0063
0064
0065 theFirstTSOSFlag = true;
0066
0067
0068 theRPCExFlag = par.getParameter<bool>("ExcludeRPCFromFit");
0069 }
0070
0071 MuonTrajectoryUpdator::MuonTrajectoryUpdator(NavigationDirection fitDirection, double chi2, int granularity)
0072 : theMaxChi2(chi2), theGranularity(granularity), theFitDirection(fitDirection) {
0073 theEstimator = new Chi2MeasurementEstimator(theMaxChi2);
0074
0075
0076 theUpdator = new KFUpdator();
0077 }
0078
0079
0080 MuonTrajectoryUpdator::~MuonTrajectoryUpdator() {
0081 delete theEstimator;
0082 delete theUpdator;
0083 }
0084
0085 void MuonTrajectoryUpdator::makeFirstTime() { theFirstTSOSFlag = true; }
0086
0087 pair<bool, TrajectoryStateOnSurface> MuonTrajectoryUpdator::update(const TrajectoryMeasurement* measurement,
0088 Trajectory& trajectory,
0089 const Propagator* propagator) {
0090 const std::string metname = "Muon|RecoMuon|MuonTrajectoryUpdator";
0091
0092 MuonPatternRecoDumper muonDumper;
0093
0094
0095 bool updated = false;
0096
0097 if (!measurement)
0098 return pair<bool, TrajectoryStateOnSurface>(updated, TrajectoryStateOnSurface());
0099
0100
0101 const DetLayer* detLayer = measurement->layer();
0102
0103
0104 TransientTrackingRecHit::ConstRecHitPointer muonRecHit = measurement->recHit();
0105
0106
0107 TransientTrackingRecHit::ConstRecHitContainer recHitsForFit =
0108 MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit, theGranularity);
0109
0110
0111 sort(recHitsForFit, detLayer);
0112
0113 TrajectoryStateOnSurface lastUpdatedTSOS = measurement->predictedState();
0114
0115 LogTrace(metname) << "Number of rechits for the fit: " << recHitsForFit.size() << endl;
0116
0117 TransientTrackingRecHit::ConstRecHitContainer::iterator recHit;
0118 for (recHit = recHitsForFit.begin(); recHit != recHitsForFit.end(); ++recHit) {
0119 if ((*recHit)->isValid()) {
0120
0121 TrajectoryStateOnSurface propagatedTSOS = propagateState(lastUpdatedTSOS, measurement, *recHit, propagator);
0122
0123 if (propagatedTSOS.isValid()) {
0124 pair<bool, double> thisChi2 = estimator()->estimate(propagatedTSOS, *((*recHit).get()));
0125
0126 LogTrace(metname) << "Estimation for Kalman Fit. Chi2: " << thisChi2.second;
0127
0128
0129 bool wantIncludeThisHit = true;
0130 if (theRPCExFlag && (*recHit)->geographicalId().det() == DetId::Muon &&
0131 (*recHit)->geographicalId().subdetId() == MuonSubdetId::RPC) {
0132 wantIncludeThisHit = false;
0133 LogTrace(metname)
0134 << "This is an RPC hit and the present configuration is such that it will be excluded from the fit";
0135 }
0136
0137
0138
0139
0140 if (thisChi2.first) {
0141 updated = true;
0142 if (wantIncludeThisHit) {
0143
0144 LogTrace(metname) << endl
0145 << " Kalman Start"
0146 << "\n"
0147 << "\n";
0148 LogTrace(metname) << " Meas. Position : " << (**recHit).globalPosition() << "\n"
0149 << " Pred. Position : " << propagatedTSOS.globalPosition()
0150 << " Pred Direction : " << propagatedTSOS.globalDirection() << endl;
0151
0152 if (theRescaleErrorFlag && theFirstTSOSFlag) {
0153 propagatedTSOS.rescaleError(theRescaleFactor);
0154 theFirstTSOSFlag = false;
0155 }
0156
0157 lastUpdatedTSOS = measurementUpdator()->update(propagatedTSOS, *((*recHit).get()));
0158
0159 if (!lastUpdatedTSOS.isValid()) {
0160 edm::LogInfo(metname) << "Invalid last TSOS, will skip RecHit ";
0161 lastUpdatedTSOS = propagatedTSOS;
0162 continue;
0163 }
0164
0165 LogTrace(metname) << " Fit Position : " << lastUpdatedTSOS.globalPosition()
0166 << " Fit Direction : " << lastUpdatedTSOS.globalDirection() << "\n"
0167 << " Fit position radius : " << lastUpdatedTSOS.globalPosition().perp()
0168 << "filter updated" << endl;
0169
0170 LogTrace(metname) << muonDumper.dumpTSOS(lastUpdatedTSOS);
0171
0172 LogTrace(metname) << "\n\n Kalman End"
0173 << "\n"
0174 << "\n";
0175
0176 TrajectoryMeasurement&& updatedMeasurement =
0177 updateMeasurement(propagatedTSOS, lastUpdatedTSOS, *recHit, thisChi2.second, detLayer, measurement);
0178
0179 trajectory.push(std::move(updatedMeasurement), thisChi2.second);
0180 } else {
0181 LogTrace(metname) << " Compatible RecHit with good chi2 but made with RPC when it was decided to not "
0182 "include it in the fit"
0183 << " --> trajectory NOT updated, invalid RecHit added." << endl;
0184
0185 MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr =
0186 MuonTransientTrackingRecHit::specificBuild((*recHit)->det(), (*recHit)->hit());
0187 invalidRhPtr->invalidateHit();
0188 TrajectoryMeasurement invalidRhMeasurement(
0189 propagatedTSOS, propagatedTSOS, invalidRhPtr, thisChi2.second, detLayer);
0190 trajectory.push(std::move(invalidRhMeasurement), thisChi2.second);
0191 }
0192 } else {
0193 if (useInvalidHits) {
0194 LogTrace(metname) << " Compatible RecHit with too large chi2"
0195 << " --> trajectory NOT updated, invalid RecHit added." << endl;
0196
0197 MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr =
0198 MuonTransientTrackingRecHit::specificBuild((*recHit)->det(), (*recHit)->hit());
0199 invalidRhPtr->invalidateHit();
0200 TrajectoryMeasurement invalidRhMeasurement(
0201 propagatedTSOS, propagatedTSOS, invalidRhPtr, thisChi2.second, detLayer);
0202 trajectory.push(std::move(invalidRhMeasurement), thisChi2.second);
0203 }
0204 }
0205 }
0206 }
0207 }
0208 recHitsForFit.clear();
0209 return pair<bool, TrajectoryStateOnSurface>(updated, lastUpdatedTSOS);
0210 }
0211
0212 TrajectoryStateOnSurface MuonTrajectoryUpdator::propagateState(
0213 const TrajectoryStateOnSurface& state,
0214 const TrajectoryMeasurement* measurement,
0215 const TransientTrackingRecHit::ConstRecHitPointer& current,
0216 const Propagator* propagator) const {
0217 const TransientTrackingRecHit::ConstRecHitPointer& mother = measurement->recHit();
0218
0219 if (current->geographicalId() == mother->geographicalId())
0220 return measurement->predictedState();
0221
0222 const TrajectoryStateOnSurface tsos = propagator->propagate(state, current->det()->surface());
0223 return tsos;
0224 }
0225
0226
0227 TrajectoryMeasurement MuonTrajectoryUpdator::updateMeasurement(const TrajectoryStateOnSurface& propagatedTSOS,
0228 const TrajectoryStateOnSurface& lastUpdatedTSOS,
0229 const TransientTrackingRecHit::ConstRecHitPointer& recHit,
0230 const double& chi2,
0231 const DetLayer* detLayer,
0232 const TrajectoryMeasurement* initialMeasurement) {
0233 return TrajectoryMeasurement(propagatedTSOS, lastUpdatedTSOS, recHit, chi2, detLayer);
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248 }
0249
0250 void MuonTrajectoryUpdator::sort(TransientTrackingRecHit::ConstRecHitContainer& recHitsForFit,
0251 const DetLayer* detLayer) {
0252 if (detLayer->subDetector() == GeomDetEnumerators::DT) {
0253 if (fitDirection() == insideOut)
0254 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), RadiusComparatorInOut());
0255 else if (fitDirection() == outsideIn)
0256 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), RadiusComparatorOutIn());
0257 else
0258 LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0259 }
0260
0261 else if (detLayer->subDetector() == GeomDetEnumerators::CSC) {
0262 if (fitDirection() == insideOut)
0263 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0264 else if (fitDirection() == outsideIn)
0265 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0266 else
0267 LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0268 }
0269
0270 else if (detLayer->subDetector() == GeomDetEnumerators::GEM) {
0271 if (fitDirection() == insideOut)
0272 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0273 else if (fitDirection() == outsideIn)
0274 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0275 else
0276 LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0277 } else if (detLayer->subDetector() == GeomDetEnumerators::ME0) {
0278 if (fitDirection() == insideOut)
0279 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0280 else if (fitDirection() == outsideIn)
0281 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0282 else
0283 LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0284 }
0285 }