Muons/plugins/FastTSGFromPropagation.cc

0001 #include "FastSimulation/Muons/plugins/FastTSGFromPropagation.h"
0002
0003 #include <memory>
0004
0005 /** \class FastTSGFromPropagation
0006  *
0007  *  Emulate TSGFromPropagation in RecoMuon
0008  *
0009  *  \author Hwidong Yoo - Purdue University
0010  */
0011
0012 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0013 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0014 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0015 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0016 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
0017 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0018 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0019 #include "TrackingTools/GeomPropagators/interface/StateOnTrackerBound.h"
0020
0021 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0022 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0023 // #include "RecoTracker/MeasurementDet/interface/TkStripMeasurementDet.h"
0024 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0025 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0026 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0027
0028 #include "RecoMuon/GlobalTrackingTools/interface/DirectTrackerNavigation.h"
0029 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0030
0031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0032 #include "FWCore/ServiceRegistry/interface/Service.h"
0033 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0034
0035 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0036 #include "FastSimulation/Tracking/interface/TrajectorySeedHitCandidate.h"
0037 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0038 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0039 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0040 #include "TrackingTools/TransientTrackingRecHit/interface/GenericTransientTrackingRecHit.h"
0041 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0042 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0043
0044 #include "FastSimulation/Tracking/interface/FastTrackingUtilities.h"
0045
0046 using namespace std;
0047
0048 FastTSGFromPropagation::FastTSGFromPropagation(const edm::ParameterSet& iConfig, edm::ConsumesCollector& iC)
0049     : FastTSGFromPropagation(iConfig, nullptr, iC) {}
0050
0051 FastTSGFromPropagation::FastTSGFromPropagation(const edm::ParameterSet& iConfig,
0052                                                const MuonServiceProxy* service,
0053                                                edm::ConsumesCollector& iC)
0054     : theCategory("FastSimulation|Muons|FastTSGFromPropagation"),
0055       theNavigation(),
0056       theService(service),
0057       theUpdator(),
0058       theEstimator(),
0059       theSigmaZ(0.0),
0060       theConfig(iConfig),
0061       theSimTrackCollectionToken_(
0062           iC.consumes<edm::SimTrackContainer>(theConfig.getParameter<edm::InputTag>("SimTrackCollectionLabel"))),
0063       recHitCombinationsToken_(
0064           iC.consumes<FastTrackerRecHitCombinationCollection>(theConfig.getParameter<edm::InputTag>("HitProducer"))),
0065       beamSpot_(iC.consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))),
0066       theMeasurementTrackerEventToken_(
0067           iC.consumes<MeasurementTrackerEvent>(iConfig.getParameter<edm::InputTag>("MeasurementTrackerEvent"))),
0068       theGeometryToken(iC.esConsumes<edm::Transition::BeginRun>()),
0069       theTTRHBuilderToken(iC.esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "WithTrackAngle"))),
0070       theTrackerToken(iC.esConsumes()) {}
0071
0072 FastTSGFromPropagation::~FastTSGFromPropagation() { LogTrace(theCategory) << " FastTSGFromPropagation dtor called "; }
0073
0074 void FastTSGFromPropagation::trackerSeeds(const TrackCand& staMuon,
0075                                           const TrackingRegion& region,
0076                                           const TrackerTopology* tTopo,
0077                                           std::vector<TrajectorySeed>& result) {
0078   if (theResetMethod == "discrete")
0079     getRescalingFactor(staMuon);
0080
0081   TrajectoryStateOnSurface staState = outerTkState(staMuon);
0082
0083   if (!staState.isValid()) {
0084     LogTrace(theCategory) << "Error: initial state from L2 muon is invalid.";
0085     return;
0086   }
0087
0088   LogTrace(theCategory) << "begin of trackerSeed:\n staState pos: " << staState.globalPosition()
0089                         << " mom: " << staState.globalMomentum() << "pos eta: " << staState.globalPosition().eta()
0090                         << "mom eta: " << staState.globalMomentum().eta();
0091
0092   std::vector<const DetLayer*> nls = theNavigation->compatibleLayers(*(staState.freeState()), oppositeToMomentum);
0093
0094   LogTrace(theCategory) << " compatible layers: " << nls.size();
0095
0096   if (nls.empty())
0097     return;
0098
0099   int ndesLayer = 0;
0100
0101   bool usePredictedState = false;
0102
0103   if (theUpdateStateFlag) {  //use updated states
0104     std::vector<TrajectoryMeasurement> alltm;
0105
0106     for (std::vector<const DetLayer*>::const_iterator inl = nls.begin(); inl != nls.end(); inl++, ndesLayer++) {
0107       if ((*inl == nullptr))
0108         break;
0109       //         if ( (inl != nls.end()-1 ) && ( (*inl)->subDetector() == GeomDetEnumerators::TEC ) && ( (*(inl+1))->subDetector() == GeomDetEnumerators::TOB ) ) continue;
0110       alltm = findMeasurements(*inl, staState);
0111       if ((!alltm.empty())) {
0112         LogTrace(theCategory) << "final compatible layer: " << ndesLayer;
0113         break;
0114       }
0115     }
0116
0117     if (alltm.empty()) {
0118       LogTrace(theCategory) << " NO Measurements Found: eta: " << staState.globalPosition().eta() << "pt "
0119                             << staState.globalMomentum().perp();
0120       usePredictedState = true;
0121     } else {
0122       LogTrace(theCategory) << " Measurements for seeds: " << alltm.size();
0123       std::stable_sort(alltm.begin(), alltm.end(), increasingEstimate());
0124       if (alltm.size() > 5)
0125         alltm.erase(alltm.begin() + 5, alltm.end());
0126
0127       const edm::SimTrackContainer* simTracks = &(*theSimTracks);
0128       TrajectorySeedHitCandidate theSeedHits;
0129       std::vector<TrajectorySeedHitCandidate> outerHits;
0130
0131       //std::vector<TrajectorySeed>  tmpTS;
0132       bool isMatch = false;
0133       for (std::vector<TrajectoryMeasurement>::const_iterator itm = alltm.begin(); itm != alltm.end(); itm++) {
0134         const TrajectoryStateOnSurface seedState = itm->predictedState();
0135         double preY = seedState.globalPosition().y();
0136
0137         // Check SimTrack
0138         FreeTrajectoryState simtrack_trackerstate;
0139         for (unsigned icomb = 0; icomb < recHitCombinations->size(); ++icomb) {
0140           const auto& recHitCombination = (*recHitCombinations)[icomb];
0141           if (recHitCombination.empty())
0142             continue;
0143           int32_t simTrackId = recHitCombination.back()->simTrackId(0);
0144           const SimTrack& simtrack = (*simTracks)[simTrackId];
0145
0146           GlobalPoint position(simtrack.trackerSurfacePosition().x(),
0147                                simtrack.trackerSurfacePosition().y(),
0148                                simtrack.trackerSurfacePosition().z());
0149           GlobalVector momentum(simtrack.trackerSurfaceMomentum().x(),
0150                                 simtrack.trackerSurfaceMomentum().y(),
0151                                 simtrack.trackerSurfaceMomentum().z());
0152           int charge = (int)simtrack.charge();
0153           GlobalTrajectoryParameters glb_parameters(
0154               position, momentum, charge, &*theService->magneticField().product());
0155           simtrack_trackerstate = FreeTrajectoryState(glb_parameters);
0156
0157           unsigned int outerId = 0;
0158           for (const auto& recHitRef : recHitCombination) {
0159             theSeedHits = TrajectorySeedHitCandidate(recHitRef.get(), tTopo);
0160             unsigned int id = theSeedHits.hit()->geographicalId().rawId();
0161             if (preY < 0) {
0162               if (id > outerId)
0163                 outerId = id;
0164             } else {
0165               if (id > outerId)
0166                 outerId = id;
0167             }
0168           }
0169           for (const auto& recHitRef : recHitCombination) {
0170             theSeedHits = TrajectorySeedHitCandidate(recHitRef.get(), tTopo);
0171             if (itm->recHit()->hit()->geographicalId().rawId() == theSeedHits.hit()->geographicalId().rawId()) {
0172               auto aTrackingRecHit = std::unique_ptr<TrackingRecHit>(theSeedHits.hit()->clone());
0173               TransientTrackingRecHit::ConstRecHitPointer recHit = theTTRHBuilder->build(aTrackingRecHit.get());
0174               if (!recHit)
0175                 continue;
0176               TrajectoryStateOnSurface updatedTSOS = updator()->update(seedState, *(recHit));
0177               if (updatedTSOS.isValid() && passSelection(updatedTSOS)) {
0178                 edm::OwnVector<TrackingRecHit> container;
0179                 container.push_back(recHit->hit()->clone());
0180                 fastTrackingUtilities::setRecHitCombinationIndex(container, icomb);
0181                 TrajectorySeed ts = createSeed(updatedTSOS, container, recHit->geographicalId());
0182                 // check direction
0183                 const TrajectorySeed* aSeed = &ts;
0184                 PTrajectoryStateOnDet PTSOD = aSeed->startingState();
0185
0186                 const GeomDet* g = theGeometry->idToDet(PTSOD.detId());
0187                 TrajectoryStateOnSurface tsos = trajectoryStateTransform::transientState(
0188                     PTSOD, &(g->surface()), &*theService->magneticField().product());
0189                 if (tsos.globalMomentum().basicVector() * seedState.globalMomentum().basicVector() < 0.)
0190                   continue;
0191                 result.push_back(ts);
0192                 isMatch = true;
0193               }
0194             }
0195           }
0196         }
0197       }
0198       if (!isMatch) {
0199         // if there is no hits w.r.t. TM, find outermost hit
0200         for (std::vector<TrajectoryMeasurement>::const_iterator itm = alltm.begin(); itm != alltm.end(); itm++) {
0201           const TrajectoryStateOnSurface seedState = itm->predictedState();
0202           double preY = seedState.globalPosition().y();
0203
0204           // Check SimTrack
0205           TrackingRecHit* aTrackingRecHit;
0206           FreeTrajectoryState simtrack_trackerstate;
0207
0208           for (unsigned icomb = 0; icomb < recHitCombinations->size(); ++icomb) {
0209             const auto& recHitCombination = (*recHitCombinations)[icomb];
0210             if (recHitCombination.empty())
0211               continue;
0212             int32_t simTrackId = recHitCombination.back()->simTrackId(0);
0213             const SimTrack& simtrack = (*simTracks)[simTrackId];
0214
0215             GlobalPoint position(simtrack.trackerSurfacePosition().x(),
0216                                  simtrack.trackerSurfacePosition().y(),
0217                                  simtrack.trackerSurfacePosition().z());
0218             GlobalVector momentum(simtrack.trackerSurfaceMomentum().x(),
0219                                   simtrack.trackerSurfaceMomentum().y(),
0220                                   simtrack.trackerSurfaceMomentum().z());
0221             int charge = (int)simtrack.charge();
0222             GlobalTrajectoryParameters glb_parameters(
0223                 position, momentum, charge, &*theService->magneticField().product());
0224             simtrack_trackerstate = FreeTrajectoryState(glb_parameters);
0225
0226             unsigned int outerId = 0;
0227             for (const auto& recHitRef : recHitCombination) {
0228               theSeedHits = TrajectorySeedHitCandidate(recHitRef.get(), tTopo);
0229               unsigned int id = theSeedHits.hit()->geographicalId().rawId();
0230               if (preY < 0) {
0231                 if (id > outerId)
0232                   outerId = id;
0233               } else {
0234                 if (id > outerId)
0235                   outerId = id;
0236               }
0237             }
0238             for (const auto& recHitRef : recHitCombination) {
0239               theSeedHits = TrajectorySeedHitCandidate(recHitRef.get(), tTopo);
0240               if (outerId == theSeedHits.hit()->geographicalId().rawId()) {
0241                 aTrackingRecHit = theSeedHits.hit()->clone();
0242                 TransientTrackingRecHit::ConstRecHitPointer recHit = theTTRHBuilder->build(aTrackingRecHit);
0243                 if (!recHit)
0244                   continue;
0245                 TrajectoryStateOnSurface updatedTSOS = updator()->update(seedState, *(recHit));
0246                 if (updatedTSOS.isValid() && passSelection(updatedTSOS)) {
0247                   edm::OwnVector<TrackingRecHit> container;
0248                   container.push_back(recHit->hit()->clone());
0249                   fastTrackingUtilities::setRecHitCombinationIndex(container, icomb);
0250                   TrajectorySeed ts = createSeed(updatedTSOS, container, recHit->geographicalId());
0251                   // check direction
0252                   const TrajectorySeed* aSeed = &ts;
0253                   PTrajectoryStateOnDet PTSOD = aSeed->startingState();
0254
0255                   const GeomDet* g = theGeometry->idToDet(PTSOD.detId());
0256                   TrajectoryStateOnSurface tsos = trajectoryStateTransform::transientState(
0257                       PTSOD, &(g->surface()), &*theService->magneticField().product());
0258                   if (tsos.globalMomentum().basicVector() * seedState.globalMomentum().basicVector() < 0.)
0259                     continue;
0260                   result.push_back(ts);
0261                 }
0262               }
0263             }
0264           }
0265         }
0266       }
0267
0268       /*
0269        for( unsigned ir = 0; ir < tmpTS.size(); ir++ ) {
0270      const BasicTrajectorySeed* aSeed = &((tmpTS)[ir]);
0271      PTrajectoryStateOnDet PTSOD = aSeed->startingState();
0272
0273      DetId seedDetId(PTSOD.detId());
0274      const GeomDet * g = theGeometry->idToDet(seedDetId);
0275      TrajectoryStateOnSurface tsos = trajectoryStateTransform::transientState(PTSOD, &(g->surface()),  &*theService->magneticField().product());
0276          cout << "tsos3 = " << tsos.globalMomentum() << endl;
0277      if( _index == ir ) {
0278          cout << "tsos4 = " << tsos.globalMomentum() << endl;
0279          result.push_back(tmpTS[ir]);
0280      }
0281        }
0282        */
0283       LogTrace(theCategory) << "result: " << result.size();
0284       return;
0285     }
0286   }
0287
0288   if (!theUpdateStateFlag || usePredictedState) {  //use predicted states
0289     LogTrace(theCategory) << "use predicted state: ";
0290     for (std::vector<const DetLayer*>::const_iterator inl = nls.begin(); inl != nls.end(); inl++) {
0291       if (!result.empty() || *inl == nullptr) {
0292         break;
0293       }
0294       std::vector<DetLayer::DetWithState> compatDets = (*inl)->compatibleDets(staState, *propagator(), *estimator());
0295       LogTrace(theCategory) << " compatDets " << compatDets.size();
0296       if (compatDets.empty())
0297         continue;
0298       TrajectorySeed ts = createSeed(compatDets.front().second, compatDets.front().first->geographicalId());
0299       result.push_back(ts);
0300     }
0301     LogTrace(theCategory) << "result: " << result.size();
0302     return;
0303   }
0304   return;
0305 }
0306
0307 void FastTSGFromPropagation::init(const MuonServiceProxy* service) {
0308   theMaxChi2 = theConfig.getParameter<double>("MaxChi2");
0309
0310   theFixedErrorRescaling = theConfig.getParameter<double>("ErrorRescaling");
0311
0312   theFlexErrorRescaling = 1.0;
0313
0314   theResetMethod = theConfig.getParameter<std::string>("ResetMethod");
0315
0316   if (theResetMethod != "discrete" && theResetMethod != "fixed" && theResetMethod != "matrix") {
0317     edm::LogError("FastTSGFromPropagation") << "Wrong error rescaling method: " << theResetMethod << "\n"
0318                                             << "Possible choices are: discrete, fixed, matrix.\n"
0319                                             << "Use discrete method" << std::endl;
0320     theResetMethod = "discrete";
0321   }
0322
0323   theEstimator = std::make_unique<Chi2MeasurementEstimator>(theMaxChi2);
0324
0325   theCacheId_TG = 0;
0326
0327   thePropagatorName = theConfig.getParameter<std::string>("Propagator");
0328
0329   theService = service;
0330
0331   theUseVertexStateFlag = theConfig.getParameter<bool>("UseVertexState");
0332
0333   theUpdateStateFlag = theConfig.getParameter<bool>("UpdateState");
0334
0335   theSelectStateFlag = theConfig.getParameter<bool>("SelectState");
0336
0337   theUpdator = std::make_unique<KFUpdator>();
0338
0339   theSigmaZ = theConfig.getParameter<double>("SigmaZ");
0340
0341   edm::ParameterSet errorMatrixPset = theConfig.getParameter<edm::ParameterSet>("errorMatrixPset");
0342   if (theResetMethod == "matrix" && !errorMatrixPset.empty()) {
0343     theAdjustAtIp = errorMatrixPset.getParameter<bool>("atIP");
0344     theErrorMatrixAdjuster = std::make_unique<MuonErrorMatrix>(errorMatrixPset);
0345   } else {
0346     theAdjustAtIp = false;
0347     theErrorMatrixAdjuster.reset();
0348   }
0349
0350   theGeometry = &theService->eventSetup().getData(theGeometryToken);
0351
0352   theTTRHBuilder = theService->eventSetup().getHandle(theTTRHBuilderToken);
0353 }
0354
0355 void FastTSGFromPropagation::setEvent(const edm::Event& iEvent) {
0356   iEvent.getByToken(beamSpot_, theBeamSpot);
0357
0358   // retrieve the MC truth (SimTracks)
0359   iEvent.getByToken(theSimTrackCollectionToken_, theSimTracks);
0360   iEvent.getByToken(recHitCombinationsToken_, recHitCombinations);
0361
0362   if (theUpdateStateFlag) {
0363     iEvent.getByToken(theMeasurementTrackerEventToken_, theMeasTrackerEvent);
0364   }
0365
0366   unsigned long long newCacheId_TG = theService->eventSetup().get<TrackerRecoGeometryRecord>().cacheIdentifier();
0367
0368   if (newCacheId_TG != theCacheId_TG) {
0369     LogTrace(theCategory) << "Tracker Reco Geometry changed!";
0370     theCacheId_TG = newCacheId_TG;
0371     auto theTracker = theService->eventSetup().getHandle(theTrackerToken);
0372     theNavigation = std::make_unique<DirectTrackerNavigation>(theTracker);
0373   }
0374 }
0375
0376 TrajectoryStateOnSurface FastTSGFromPropagation::innerState(const TrackCand& staMuon) const {
0377   TrajectoryStateOnSurface innerTS;
0378
0379   if (staMuon.first && staMuon.first->isValid()) {
0380     if (staMuon.first->direction() == alongMomentum) {
0381       innerTS = staMuon.first->firstMeasurement().updatedState();
0382     } else if (staMuon.first->direction() == oppositeToMomentum) {
0383       innerTS = staMuon.first->lastMeasurement().updatedState();
0384     }
0385   } else {
0386     innerTS = trajectoryStateTransform::innerStateOnSurface(
0387         *(staMuon.second), *theService->trackingGeometry(), &*theService->magneticField());
0388   }
0389   //rescale the error
0390   adjust(innerTS);
0391
0392   return innerTS;
0393 }
0394
0395 TrajectoryStateOnSurface FastTSGFromPropagation::outerTkState(const TrackCand& staMuon) const {
0396   TrajectoryStateOnSurface result;
0397
0398   if (theUseVertexStateFlag && staMuon.second->pt() > 1.0) {
0399     FreeTrajectoryState iniState =
0400         trajectoryStateTransform::initialFreeState(*(staMuon.second), &*theService->magneticField());
0401     //rescale the error at IP
0402     adjust(iniState);
0403
0404     StateOnTrackerBound fromInside(&*(theService->propagator("PropagatorWithMaterial")));
0405     result = fromInside(iniState);
0406   } else {
0407     StateOnTrackerBound fromOutside(&*propagator());
0408     result = fromOutside(innerState(staMuon));
0409   }
0410   return result;
0411 }
0412
0413 TrajectorySeed FastTSGFromPropagation::createSeed(const TrajectoryStateOnSurface& tsos, const DetId& id) const {
0414   edm::OwnVector<TrackingRecHit> container;
0415   return createSeed(tsos, container, id);
0416 }
0417
0418 TrajectorySeed FastTSGFromPropagation::createSeed(const TrajectoryStateOnSurface& tsos,
0419                                                   const edm::OwnVector<TrackingRecHit>& container,
0420                                                   const DetId& id) const {
0421   PTrajectoryStateOnDet seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());
0422   return TrajectorySeed(seedTSOS, container, oppositeToMomentum);
0423 }
0424
0425 void FastTSGFromPropagation::validMeasurements(std::vector<TrajectoryMeasurement>& tms) const {
0426   std::vector<TrajectoryMeasurement>::iterator tmsend = std::remove_if(tms.begin(), tms.end(), isInvalid());
0427   tms.erase(tmsend, tms.end());
0428   return;
0429 }
0430
0431 std::vector<TrajectoryMeasurement> FastTSGFromPropagation::findMeasurements(
0432     const DetLayer* nl, const TrajectoryStateOnSurface& staState) const {
0433   std::vector<TrajectoryMeasurement> result;
0434
0435   std::vector<DetLayer::DetWithState> compatDets = nl->compatibleDets(staState, *propagator(), *estimator());
0436   if (compatDets.empty())
0437     return result;
0438
0439   for (std::vector<DetLayer::DetWithState>::const_iterator idws = compatDets.begin(); idws != compatDets.end();
0440        ++idws) {
0441     if (idws->second.isValid() && (idws->first)) {
0442       std::vector<TrajectoryMeasurement> tmptm =
0443           theMeasTrackerEvent->idToDet(idws->first->geographicalId())
0444               .fastMeasurements(idws->second, idws->second, *propagator(), *estimator());
0445       //validMeasurements(tmptm);
0446       //         if ( tmptm.size() > 2 ) {
0447       //            std::stable_sort(tmptm.begin(),tmptm.end(),increasingEstimate());
0448       //            result.insert(result.end(),tmptm.begin(), tmptm.begin()+2);
0449       //         } else {
0450       result.insert(result.end(), tmptm.begin(), tmptm.end());
0451       //         }
0452     }
0453   }
0454
0455   return result;
0456 }
0457
0458 bool FastTSGFromPropagation::passSelection(const TrajectoryStateOnSurface& tsos) const {
0459   if (!theSelectStateFlag)
0460     return true;
0461   else {
0462     if (theBeamSpot.isValid()) {
0463       return ((fabs(zDis(tsos) - theBeamSpot->z0()) < theSigmaZ));
0464
0465     } else {
0466       return ((fabs(zDis(tsos)) < theSigmaZ));
0467       //      double theDxyCut = 100;
0468       //      return ( (zDis(tsos) < theSigmaZ) && (dxyDis(tsos) < theDxyCut) );
0469     }
0470   }
0471 }
0472
0473 double FastTSGFromPropagation::dxyDis(const TrajectoryStateOnSurface& tsos) const {
0474   return fabs(
0475       (-tsos.globalPosition().x() * tsos.globalMomentum().y() + tsos.globalPosition().y() * tsos.globalMomentum().x()) /
0476       tsos.globalMomentum().perp());
0477 }
0478
0479 double FastTSGFromPropagation::zDis(const TrajectoryStateOnSurface& tsos) const {
0480   return tsos.globalPosition().z() -
0481          tsos.globalPosition().perp() * tsos.globalMomentum().z() / tsos.globalMomentum().perp();
0482 }
0483
0484 void FastTSGFromPropagation::getRescalingFactor(const TrackCand& staMuon) {
0485   float pt = (staMuon.second)->pt();
0486   if (pt < 13.0)
0487     theFlexErrorRescaling = 3;
0488   else if (pt < 30.0)
0489     theFlexErrorRescaling = 5;
0490   else
0491     theFlexErrorRescaling = 10;
0492   return;
0493 }
0494
0495 void FastTSGFromPropagation::adjust(FreeTrajectoryState& state) const {
0496   //rescale the error
0497   if (theResetMethod == "discreate") {
0498     state.rescaleError(theFlexErrorRescaling);
0499     return;
0500   }
0501
0502   //rescale the error
0503   if (theResetMethod == "fixed" || !theErrorMatrixAdjuster) {
0504     state.rescaleError(theFixedErrorRescaling);
0505     return;
0506   }
0507
0508   CurvilinearTrajectoryError oMat = state.curvilinearError();
0509   CurvilinearTrajectoryError sfMat = theErrorMatrixAdjuster->get(state.momentum());  //FIXME with position
0510   MuonErrorMatrix::multiply(oMat, sfMat);
0511
0512   state = FreeTrajectoryState(state.parameters(), oMat);
0513 }
0514
0515 void FastTSGFromPropagation::adjust(TrajectoryStateOnSurface& state) const {
0516   //rescale the error
0517   if (theResetMethod == "discreate") {
0518     state.rescaleError(theFlexErrorRescaling);
0519     return;
0520   }
0521
0522   if (theResetMethod == "fixed" || !theErrorMatrixAdjuster) {
0523     state.rescaleError(theFixedErrorRescaling);
0524     return;
0525   }
0526
0527   CurvilinearTrajectoryError oMat = state.curvilinearError();
0528   CurvilinearTrajectoryError sfMat = theErrorMatrixAdjuster->get(state.globalMomentum());  //FIXME with position
0529   MuonErrorMatrix::multiply(oMat, sfMat);
0530
0531   state =
0532       TrajectoryStateOnSurface(state.weight(), state.globalParameters(), oMat, state.surface(), state.surfaceSide());
0533 }
0534
0535 void FastTSGFromPropagation::stateOnDet(const TrajectoryStateOnSurface& ts,
0536                                         unsigned int detid,
0537                                         PTrajectoryStateOnDet& pts) const {
0538   const AlgebraicSymMatrix55& m = ts.localError().matrix();
0539   int dim = 5;  /// should check if corresponds to m
0540   float localErrors[15];
0541   int k = 0;
0542   for (int i = 0; i < dim; ++i) {
0543     for (int j = 0; j <= i; ++j) {
0544       localErrors[k++] = m(i, j);
0545     }
0546   }
0547   int surfaceSide = static_cast<int>(ts.surfaceSide());
0548   pts = PTrajectoryStateOnDet(ts.localParameters(), ts.globalMomentum().perp(), localErrors, detid, surfaceSide);
0549 }