File indexing completed on 2023-03-17 11:20:57
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 LogTrace(metname) << " Fit Position : " << lastUpdatedTSOS.globalPosition()
0160 << " Fit Direction : " << lastUpdatedTSOS.globalDirection() << "\n"
0161 << " Fit position radius : " << lastUpdatedTSOS.globalPosition().perp()
0162 << "filter updated" << endl;
0163
0164 LogTrace(metname) << muonDumper.dumpTSOS(lastUpdatedTSOS);
0165
0166 LogTrace(metname) << "\n\n Kalman End"
0167 << "\n"
0168 << "\n";
0169
0170 TrajectoryMeasurement&& updatedMeasurement =
0171 updateMeasurement(propagatedTSOS, lastUpdatedTSOS, *recHit, thisChi2.second, detLayer, measurement);
0172
0173 trajectory.push(std::move(updatedMeasurement), thisChi2.second);
0174 } else {
0175 LogTrace(metname) << " Compatible RecHit with good chi2 but made with RPC when it was decided to not "
0176 "include it in the fit"
0177 << " --> trajectory NOT updated, invalid RecHit added." << endl;
0178
0179 MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr =
0180 MuonTransientTrackingRecHit::specificBuild((*recHit)->det(), (*recHit)->hit());
0181 invalidRhPtr->invalidateHit();
0182 TrajectoryMeasurement invalidRhMeasurement(
0183 propagatedTSOS, propagatedTSOS, invalidRhPtr, thisChi2.second, detLayer);
0184 trajectory.push(std::move(invalidRhMeasurement), thisChi2.second);
0185 }
0186 } else {
0187 if (useInvalidHits) {
0188 LogTrace(metname) << " Compatible RecHit with too large chi2"
0189 << " --> trajectory NOT updated, invalid RecHit added." << endl;
0190
0191 MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr =
0192 MuonTransientTrackingRecHit::specificBuild((*recHit)->det(), (*recHit)->hit());
0193 invalidRhPtr->invalidateHit();
0194 TrajectoryMeasurement invalidRhMeasurement(
0195 propagatedTSOS, propagatedTSOS, invalidRhPtr, thisChi2.second, detLayer);
0196 trajectory.push(std::move(invalidRhMeasurement), thisChi2.second);
0197 }
0198 }
0199 }
0200 }
0201 }
0202 recHitsForFit.clear();
0203 return pair<bool, TrajectoryStateOnSurface>(updated, lastUpdatedTSOS);
0204 }
0205
0206 TrajectoryStateOnSurface MuonTrajectoryUpdator::propagateState(
0207 const TrajectoryStateOnSurface& state,
0208 const TrajectoryMeasurement* measurement,
0209 const TransientTrackingRecHit::ConstRecHitPointer& current,
0210 const Propagator* propagator) const {
0211 const TransientTrackingRecHit::ConstRecHitPointer& mother = measurement->recHit();
0212
0213 if (current->geographicalId() == mother->geographicalId())
0214 return measurement->predictedState();
0215
0216 const TrajectoryStateOnSurface tsos = propagator->propagate(state, current->det()->surface());
0217 return tsos;
0218 }
0219
0220
0221 TrajectoryMeasurement MuonTrajectoryUpdator::updateMeasurement(const TrajectoryStateOnSurface& propagatedTSOS,
0222 const TrajectoryStateOnSurface& lastUpdatedTSOS,
0223 const TransientTrackingRecHit::ConstRecHitPointer& recHit,
0224 const double& chi2,
0225 const DetLayer* detLayer,
0226 const TrajectoryMeasurement* initialMeasurement) {
0227 return TrajectoryMeasurement(propagatedTSOS, lastUpdatedTSOS, recHit, chi2, detLayer);
0228
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242 }
0243
0244 void MuonTrajectoryUpdator::sort(TransientTrackingRecHit::ConstRecHitContainer& recHitsForFit,
0245 const DetLayer* detLayer) {
0246 if (detLayer->subDetector() == GeomDetEnumerators::DT) {
0247 if (fitDirection() == insideOut)
0248 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), RadiusComparatorInOut());
0249 else if (fitDirection() == outsideIn)
0250 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), RadiusComparatorOutIn());
0251 else
0252 LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0253 }
0254
0255 else if (detLayer->subDetector() == GeomDetEnumerators::CSC) {
0256 if (fitDirection() == insideOut)
0257 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0258 else if (fitDirection() == outsideIn)
0259 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0260 else
0261 LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0262 }
0263
0264 else if (detLayer->subDetector() == GeomDetEnumerators::GEM) {
0265 if (fitDirection() == insideOut)
0266 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0267 else if (fitDirection() == outsideIn)
0268 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0269 else
0270 LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0271 } else if (detLayer->subDetector() == GeomDetEnumerators::ME0) {
0272 if (fitDirection() == insideOut)
0273 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0274 else if (fitDirection() == outsideIn)
0275 stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0276 else
0277 LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0278 }
0279 }