TrackerSeedGenerator/plugins/TSGFromPropagation.cc

0001 #include "RecoMuon/TrackerSeedGenerator/plugins/TSGFromPropagation.h"
0002
0003 /** \class TSGFromPropagation
0004  *
0005  *  \author Chang Liu - Purdue University
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) {  //use updated states
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       //         if ( (inl != nls.end()-1 ) && ( (*inl)->subDetector() == GeomDetEnumerators::TEC ) && ( (*(inl+1))->subDetector() == GeomDetEnumerators::TOB ) ) continue;
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) {  //use predicted states
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   //rescale the error
0202   adjust(innerTS);
0203
0204   return innerTS;
0205
0206   //    return trajectoryStateTransform::innerStateOnSurface(*(staMuon.second),*theService->trackingGeometry(), &*theService->magneticField());
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     //rescale the error at IP
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       //         if ( tmptm.size() > 2 ) {
0261       //            std::stable_sort(tmptm.begin(),tmptm.end(),increasingEstimate());
0262       //            result.insert(result.end(),tmptm.begin(), tmptm.begin()+2);
0263       //         } else {
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       //      double theDxyCut = 100;
0282       //      return ( (zDis(tsos) < theSigmaZ) && (dxyDis(tsos) < theDxyCut) );
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   //rescale the error
0311   if (theResetMethod == ResetMethod::discrete) {
0312     state.rescaleError(theFlexErrorRescaling);
0313     return;
0314   }
0315
0316   //rescale the error
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());  //FIXME with position
0324   MuonErrorMatrix::multiply(oMat, sfMat);
0325
0326   state = FreeTrajectoryState(state.parameters(), oMat);
0327 }
0328
0329 void TSGFromPropagation::adjust(TrajectoryStateOnSurface& state) const {
0330   //rescale the error
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());  //FIXME with position
0343   MuonErrorMatrix::multiply(oMat, sfMat);
0344
0345   state =
0346       TrajectoryStateOnSurface(state.weight(), state.globalParameters(), oMat, state.surface(), state.surfaceSide());
0347 }