File indexing completed on 2023-03-17 11:20:54
0001 #include "RecoMuon/TrackerSeedGenerator/plugins/TSGFromPropagation.h"
0002
0003
0004
0005
0006
0007
0008 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0009 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0010 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0011 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0012 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
0013
0014 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0015 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0016 #include "TrackingTools/GeomPropagators/interface/StateOnTrackerBound.h"
0017
0018 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0019 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0020 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0021 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0022
0023 #include "RecoMuon/GlobalTrackingTools/interface/DirectTrackerNavigation.h"
0024 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0025 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0026
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028
0029 TSGFromPropagation::TSGFromPropagation(const edm::ParameterSet& iConfig, edm::ConsumesCollector& iC)
0030 : TSGFromPropagation(iConfig, iC, nullptr) {}
0031
0032 TSGFromPropagation::TSGFromPropagation(const edm::ParameterSet& iConfig,
0033 edm::ConsumesCollector& iC,
0034 const MuonServiceProxy* service)
0035 : theCategory("Muon|RecoMuon|TSGFromPropagation"),
0036 theService(service),
0037 theMaxChi2(iConfig.getParameter<double>("MaxChi2")),
0038 theFixedErrorRescaling(iConfig.getParameter<double>("ErrorRescaling")),
0039 theUseVertexStateFlag(iConfig.getParameter<bool>("UseVertexState")),
0040 theUpdateStateFlag(iConfig.getParameter<bool>("UpdateState")),
0041 theResetMethod([](const edm::ParameterSet& iConfig) {
0042 auto resetMethod = iConfig.getParameter<std::string>("ResetMethod");
0043 if (resetMethod != "discrete" && resetMethod != "fixed" && resetMethod != "matrix") {
0044 edm::LogError("TSGFromPropagation") << "Wrong error rescaling method: " << resetMethod << "\n"
0045 << "Possible choices are: discrete, fixed, matrix.\n"
0046 << "Use discrete method" << std::endl;
0047 resetMethod = "discrete";
0048 }
0049 if ("fixed" == resetMethod) {
0050 return ResetMethod::fixed;
0051 }
0052 if ("matrix" == resetMethod) {
0053 return ResetMethod::matrix;
0054 }
0055 return ResetMethod::discrete;
0056 }(iConfig)),
0057 theSelectStateFlag(iConfig.getParameter<bool>("SelectState")),
0058 thePropagatorName(iConfig.getParameter<std::string>("Propagator")),
0059 theSigmaZ(iConfig.getParameter<double>("SigmaZ")),
0060 theErrorMatrixPset(iConfig.getParameter<edm::ParameterSet>("errorMatrixPset")),
0061 theBeamSpotToken(iC.consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))),
0062 theMeasurementTrackerEventToken(
0063 iC.consumes<MeasurementTrackerEvent>(iConfig.getParameter<edm::InputTag>("MeasurementTrackerEvent"))),
0064 theTrackerToken(iC.esConsumes()) {}
0065
0066 TSGFromPropagation::~TSGFromPropagation() { LogTrace(theCategory) << " TSGFromPropagation dtor called "; }
0067
0068 void TSGFromPropagation::trackerSeeds(const TrackCand& staMuon,
0069 const TrackingRegion& region,
0070 const TrackerTopology* tTopo,
0071 std::vector<TrajectorySeed>& result) {
0072 if (theResetMethod == ResetMethod::discrete)
0073 getRescalingFactor(staMuon);
0074
0075 TrajectoryStateOnSurface staState = outerTkState(staMuon);
0076
0077 if (!staState.isValid()) {
0078 LogTrace(theCategory) << "Error: initial state from L2 muon is invalid.";
0079 return;
0080 }
0081
0082 LogTrace(theCategory) << "begin of trackerSeed:\n staState pos: " << staState.globalPosition()
0083 << " mom: " << staState.globalMomentum() << "pos eta: " << staState.globalPosition().eta()
0084 << "mom eta: " << staState.globalMomentum().eta();
0085
0086 std::vector<const DetLayer*> nls = theNavigation->compatibleLayers(*(staState.freeState()), oppositeToMomentum);
0087
0088 LogTrace(theCategory) << " compatible layers: " << nls.size();
0089
0090 if (nls.empty())
0091 return;
0092
0093 int ndesLayer = 0;
0094
0095 bool usePredictedState = false;
0096
0097 if (theUpdateStateFlag) {
0098 std::vector<TrajectoryMeasurement> alltm;
0099
0100 for (std::vector<const DetLayer*>::const_iterator inl = nls.begin(); inl != nls.end(); inl++, ndesLayer++) {
0101 if ((*inl == nullptr))
0102 break;
0103
0104 alltm = findMeasurements(*inl, staState);
0105 if ((!alltm.empty())) {
0106 LogTrace(theCategory) << "final compatible layer: " << ndesLayer;
0107 break;
0108 }
0109 }
0110
0111 if (alltm.empty()) {
0112 LogTrace(theCategory) << " NO Measurements Found: eta: " << staState.globalPosition().eta() << "pt "
0113 << staState.globalMomentum().perp();
0114 usePredictedState = true;
0115 } else {
0116 LogTrace(theCategory) << " Measurements for seeds: " << alltm.size();
0117 std::stable_sort(alltm.begin(), alltm.end(), increasingEstimate());
0118 if (alltm.size() > 5)
0119 alltm.erase(alltm.begin() + 5, alltm.end());
0120
0121 int i = 0;
0122 for (std::vector<TrajectoryMeasurement>::const_iterator itm = alltm.begin(); itm != alltm.end(); itm++, i++) {
0123 TrajectoryStateOnSurface updatedTSOS = updator()->update(itm->predictedState(), *(itm->recHit()));
0124 if (updatedTSOS.isValid() && passSelection(updatedTSOS)) {
0125 edm::OwnVector<TrackingRecHit> container;
0126 container.push_back(itm->recHit()->hit()->clone());
0127 TrajectorySeed ts = createSeed(updatedTSOS, container, itm->recHit()->geographicalId());
0128 result.push_back(ts);
0129 }
0130 }
0131 LogTrace(theCategory) << "result: " << result.size();
0132 return;
0133 }
0134 }
0135
0136 if (!theUpdateStateFlag || usePredictedState) {
0137 LogTrace(theCategory) << "use predicted state: ";
0138 for (std::vector<const DetLayer*>::const_iterator inl = nls.begin(); inl != nls.end(); inl++) {
0139 if (!result.empty() || *inl == nullptr) {
0140 break;
0141 }
0142 std::vector<DetLayer::DetWithState> compatDets = (*inl)->compatibleDets(staState, *propagator(), *estimator());
0143 LogTrace(theCategory) << " compatDets " << compatDets.size();
0144 if (compatDets.empty())
0145 continue;
0146 TrajectorySeed ts = createSeed(compatDets.front().second, compatDets.front().first->geographicalId());
0147 result.push_back(ts);
0148 }
0149 LogTrace(theCategory) << "result: " << result.size();
0150 return;
0151 }
0152 return;
0153 }
0154
0155 void TSGFromPropagation::init(const MuonServiceProxy* service) {
0156 theFlexErrorRescaling = 1.0;
0157
0158 theEstimator = std::make_unique<Chi2MeasurementEstimator>(theMaxChi2);
0159
0160 theCacheId_TG = 0;
0161
0162 theService = service;
0163
0164 theUpdator = std::make_unique<KFUpdator>();
0165
0166 if (theResetMethod == ResetMethod::matrix && !theErrorMatrixPset.empty()) {
0167 theErrorMatrixAdjuster = std::make_unique<MuonErrorMatrix>(theErrorMatrixPset);
0168 }
0169 }
0170
0171 void TSGFromPropagation::setEvent(const edm::Event& iEvent) {
0172 iEvent.getByToken(theBeamSpotToken, beamSpot);
0173
0174 if (theUpdateStateFlag) {
0175 iEvent.getByToken(theMeasurementTrackerEventToken, theMeasTrackerEvent);
0176 }
0177
0178 unsigned long long newCacheId_TG = theService->eventSetup().get<TrackerRecoGeometryRecord>().cacheIdentifier();
0179
0180 if (newCacheId_TG != theCacheId_TG) {
0181 LogTrace(theCategory) << "Tracker Reco Geometry changed!";
0182 theCacheId_TG = newCacheId_TG;
0183 edm::ESHandle<GeometricSearchTracker> theTracker = theService->eventSetup().getHandle(theTrackerToken);
0184 theNavigation = std::make_unique<DirectTrackerNavigation>(theTracker);
0185 }
0186 }
0187
0188 TrajectoryStateOnSurface TSGFromPropagation::innerState(const TrackCand& staMuon) const {
0189 TrajectoryStateOnSurface innerTS;
0190
0191 if (staMuon.first && staMuon.first->isValid()) {
0192 if (staMuon.first->direction() == alongMomentum) {
0193 innerTS = staMuon.first->firstMeasurement().updatedState();
0194 } else if (staMuon.first->direction() == oppositeToMomentum) {
0195 innerTS = staMuon.first->lastMeasurement().updatedState();
0196 }
0197 } else {
0198 innerTS = trajectoryStateTransform::innerStateOnSurface(
0199 *(staMuon.second), *theService->trackingGeometry(), &*theService->magneticField());
0200 }
0201
0202 adjust(innerTS);
0203
0204 return innerTS;
0205
0206
0207 }
0208
0209 TrajectoryStateOnSurface TSGFromPropagation::outerTkState(const TrackCand& staMuon) const {
0210 TrajectoryStateOnSurface result;
0211
0212 if (theUseVertexStateFlag && staMuon.second->pt() > 1.0) {
0213 FreeTrajectoryState iniState =
0214 trajectoryStateTransform::initialFreeState(*(staMuon.second), &*theService->magneticField());
0215
0216 adjust(iniState);
0217
0218 StateOnTrackerBound fromInside(&*(theService->propagator("PropagatorWithMaterial")));
0219 result = fromInside(iniState);
0220 } else {
0221 StateOnTrackerBound fromOutside(&*propagator());
0222 result = fromOutside(innerState(staMuon));
0223 }
0224 return result;
0225 }
0226
0227 TrajectorySeed TSGFromPropagation::createSeed(const TrajectoryStateOnSurface& tsos, const DetId& id) const {
0228 edm::OwnVector<TrackingRecHit> container;
0229 return createSeed(tsos, container, id);
0230 }
0231
0232 TrajectorySeed TSGFromPropagation::createSeed(const TrajectoryStateOnSurface& tsos,
0233 const edm::OwnVector<TrackingRecHit>& container,
0234 const DetId& id) const {
0235 PTrajectoryStateOnDet const& seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());
0236 return TrajectorySeed(seedTSOS, container, oppositeToMomentum);
0237 }
0238
0239 void TSGFromPropagation::validMeasurements(std::vector<TrajectoryMeasurement>& tms) const {
0240 std::vector<TrajectoryMeasurement>::iterator tmsend = std::remove_if(tms.begin(), tms.end(), isInvalid());
0241 tms.erase(tmsend, tms.end());
0242 return;
0243 }
0244
0245 std::vector<TrajectoryMeasurement> TSGFromPropagation::findMeasurements(
0246 const DetLayer* nl, const TrajectoryStateOnSurface& staState) const {
0247 std::vector<TrajectoryMeasurement> result;
0248
0249 std::vector<DetLayer::DetWithState> compatDets = nl->compatibleDets(staState, *propagator(), *estimator());
0250 if (compatDets.empty())
0251 return result;
0252
0253 for (std::vector<DetLayer::DetWithState>::const_iterator idws = compatDets.begin(); idws != compatDets.end();
0254 ++idws) {
0255 if (idws->second.isValid() && (idws->first)) {
0256 std::vector<TrajectoryMeasurement> tmptm =
0257 theMeasTrackerEvent->idToDet(idws->first->geographicalId())
0258 .fastMeasurements(idws->second, idws->second, *propagator(), *estimator());
0259 validMeasurements(tmptm);
0260
0261
0262
0263
0264 result.insert(result.end(), tmptm.begin(), tmptm.end());
0265
0266 }
0267 }
0268
0269 return result;
0270 }
0271
0272 bool TSGFromPropagation::passSelection(const TrajectoryStateOnSurface& tsos) const {
0273 if (!theSelectStateFlag)
0274 return true;
0275 else {
0276 if (beamSpot.isValid()) {
0277 return ((fabs(zDis(tsos) - beamSpot->z0()) < theSigmaZ));
0278
0279 } else {
0280 return ((fabs(zDis(tsos)) < theSigmaZ));
0281
0282
0283 }
0284 }
0285 }
0286
0287 double TSGFromPropagation::dxyDis(const TrajectoryStateOnSurface& tsos) const {
0288 return fabs(
0289 (-tsos.globalPosition().x() * tsos.globalMomentum().y() + tsos.globalPosition().y() * tsos.globalMomentum().x()) /
0290 tsos.globalMomentum().perp());
0291 }
0292
0293 double TSGFromPropagation::zDis(const TrajectoryStateOnSurface& tsos) const {
0294 return tsos.globalPosition().z() -
0295 tsos.globalPosition().perp() * tsos.globalMomentum().z() / tsos.globalMomentum().perp();
0296 }
0297
0298 void TSGFromPropagation::getRescalingFactor(const TrackCand& staMuon) {
0299 float pt = (staMuon.second)->pt();
0300 if (pt < 13.0)
0301 theFlexErrorRescaling = 3;
0302 else if (pt < 30.0)
0303 theFlexErrorRescaling = 5;
0304 else
0305 theFlexErrorRescaling = 10;
0306 return;
0307 }
0308
0309 void TSGFromPropagation::adjust(FreeTrajectoryState& state) const {
0310
0311 if (theResetMethod == ResetMethod::discrete) {
0312 state.rescaleError(theFlexErrorRescaling);
0313 return;
0314 }
0315
0316
0317 if (theResetMethod == ResetMethod::fixed || !theErrorMatrixAdjuster) {
0318 state.rescaleError(theFixedErrorRescaling);
0319 return;
0320 }
0321
0322 CurvilinearTrajectoryError oMat = state.curvilinearError();
0323 CurvilinearTrajectoryError sfMat = theErrorMatrixAdjuster->get(state.momentum());
0324 MuonErrorMatrix::multiply(oMat, sfMat);
0325
0326 state = FreeTrajectoryState(state.parameters(), oMat);
0327 }
0328
0329 void TSGFromPropagation::adjust(TrajectoryStateOnSurface& state) const {
0330
0331 if (theResetMethod == ResetMethod::discrete) {
0332 state.rescaleError(theFlexErrorRescaling);
0333 return;
0334 }
0335
0336 if (theResetMethod == ResetMethod::fixed || !theErrorMatrixAdjuster) {
0337 state.rescaleError(theFixedErrorRescaling);
0338 return;
0339 }
0340
0341 CurvilinearTrajectoryError oMat = state.curvilinearError();
0342 CurvilinearTrajectoryError sfMat = theErrorMatrixAdjuster->get(state.globalMomentum());
0343 MuonErrorMatrix::multiply(oMat, sfMat);
0344
0345 state =
0346 TrajectoryStateOnSurface(state.weight(), state.globalParameters(), oMat, state.surface(), state.surfaceSide());
0347 }