Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:33:10

0001 /** \class MuonTrackAnalyzer
0002  *  Analyzer of the Muon tracks
0003  *
0004  *  \author R. Bellan - INFN Torino <riccardo.bellan@cern.ch>
0005  */
0006 
0007 #include "Validation/RecoMuon/src/MuonTrackAnalyzer.h"
0008 
0009 // Collaborating Class Header
0010 #include "FWCore/Framework/interface/ConsumesCollector.h"
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 
0016 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0017 
0018 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0019 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0020 #include "DataFormats/Math/interface/deltaR.h"
0021 
0022 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0023 
0024 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0025 
0026 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0027 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0028 #include "RecoMuon/TrackingTools/interface/MuonUpdatorAtVertex.h"
0029 
0030 #include "Validation/RecoMuon/src/Histograms.h"
0031 #include "Validation/RecoMuon/src/HTrack.h"
0032 
0033 #include "TFile.h"
0034 #include "TH1F.h"
0035 #include "TH2F.h"
0036 
0037 using namespace std;
0038 using namespace edm;
0039 
0040 /// Constructor
0041 MuonTrackAnalyzer::MuonTrackAnalyzer(const ParameterSet &ps) {
0042   // service parameters
0043   pset = ps;
0044   ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
0045   // the services
0046   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0047 
0048   theSimTracksLabel = edm::InputTag("g4SimHits");
0049   theSimTracksToken = consumes<edm::SimTrackContainer>(theSimTracksLabel);
0050 
0051   theTracksLabel = pset.getParameter<InputTag>("Tracks");
0052   theTracksToken = consumes<reco::TrackCollection>(theTracksLabel);
0053   doTracksAnalysis = pset.getUntrackedParameter<bool>("DoTracksAnalysis", true);
0054 
0055   doSeedsAnalysis = pset.getUntrackedParameter<bool>("DoSeedsAnalysis", false);
0056   if (doSeedsAnalysis) {
0057     theSeedsLabel = pset.getParameter<InputTag>("MuonSeed");
0058     theSeedsToken = consumes<TrajectorySeedCollection>(theSeedsLabel);
0059     ParameterSet updatorPar = pset.getParameter<ParameterSet>("MuonUpdatorAtVertexParameters");
0060     theSeedPropagatorName = updatorPar.getParameter<string>("Propagator");
0061 
0062     theUpdator = new MuonUpdatorAtVertex(updatorPar, theService);
0063   }
0064 
0065   theCSCSimHitLabel = pset.getParameter<InputTag>("CSCSimHit");
0066   theDTSimHitLabel = pset.getParameter<InputTag>("DTSimHit");
0067   theRPCSimHitLabel = pset.getParameter<InputTag>("RPCSimHit");
0068   theCSCSimHitToken = consumes<std::vector<PSimHit> >(theCSCSimHitLabel);
0069   theDTSimHitToken = consumes<std::vector<PSimHit> >(theDTSimHitLabel);
0070   theRPCSimHitToken = consumes<std::vector<PSimHit> >(theRPCSimHitLabel);
0071 
0072   theEtaRange = (EtaRange)pset.getParameter<int>("EtaRange");
0073 
0074   // number of sim tracks
0075   numberOfSimTracks = 0;
0076   // number of reco tracks
0077   numberOfRecTracks = 0;
0078 
0079   dbe_ = edm::Service<DQMStore>().operator->();
0080   out = pset.getUntrackedParameter<string>("rootFileName");
0081   dirName_ = pset.getUntrackedParameter<std::string>("dirName");
0082   subsystemname_ = pset.getUntrackedParameter<std::string>("subSystemFolder", "YourSubsystem");
0083 }
0084 
0085 /// Destructor
0086 MuonTrackAnalyzer::~MuonTrackAnalyzer() {
0087   if (theService)
0088     delete theService;
0089 }
0090 
0091 void MuonTrackAnalyzer::bookHistograms(DQMStore::IBooker &ibooker,
0092                                        edm::Run const &iRun,
0093                                        edm::EventSetup const & /* iSetup */) {
0094   ibooker.cd();
0095 
0096   InputTag algo = theTracksLabel;
0097   string dirName = dirName_;
0098   if (!algo.process().empty())
0099     dirName += algo.process() + "_";
0100   if (!algo.label().empty())
0101     dirName += algo.label() + "_";
0102   if (!algo.instance().empty())
0103     dirName += algo.instance() + "";
0104   if (dirName.find("Tracks") < dirName.length()) {
0105     dirName.replace(dirName.find("Tracks"), 6, "");
0106   }
0107   std::replace(dirName.begin(), dirName.end(), ':', '_');
0108   ibooker.setCurrentFolder(dirName);
0109 
0110   //ibooker.goUp();
0111   std::string simName = dirName;
0112   simName += "/SimTracks";
0113   hSimTracks = new HTrackVariables(ibooker, simName, "SimTracks");
0114 
0115   ibooker.cd();
0116   ibooker.setCurrentFolder(dirName);
0117 
0118   // Create the root file
0119   //theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0120 
0121   if (doSeedsAnalysis) {
0122     ibooker.cd();
0123     ibooker.setCurrentFolder(dirName);
0124     hRecoSeedInner = new HTrack(ibooker, dirName, "RecoSeed", "Inner");
0125     hRecoSeedPCA = new HTrack(ibooker, dirName, "RecoSeed", "PCA");
0126   }
0127 
0128   if (doTracksAnalysis) {
0129     ibooker.cd();
0130     ibooker.setCurrentFolder(dirName);
0131     hRecoTracksPCA = new HTrack(ibooker, dirName, "RecoTracks", "PCA");
0132     hRecoTracksInner = new HTrack(ibooker, dirName, "RecoTracks", "Inner");
0133     hRecoTracksOuter = new HTrack(ibooker, dirName, "RecoTracks", "Outer");
0134 
0135     ibooker.cd();
0136     ibooker.setCurrentFolder(dirName);
0137 
0138     // General Histos
0139 
0140     hChi2 = ibooker.book1D("chi2", "#chi^2", 200, 0, 200);
0141     hChi2VsEta = ibooker.book2D("chi2VsEta", "#chi^2 VS #eta", 120, -3., 3., 200, 0, 200);
0142 
0143     hChi2Norm = ibooker.book1D("chi2Norm", "Normalized #chi^2", 400, 0, 100);
0144     hChi2NormVsEta = ibooker.book2D("chi2NormVsEta", "Normalized #chi^2 VS #eta", 120, -3., 3., 400, 0, 100);
0145 
0146     hHitsPerTrack = ibooker.book1D("HitsPerTrack", "Number of hits per track", 55, 0, 55);
0147     hHitsPerTrackVsEta =
0148         ibooker.book2D("HitsPerTrackVsEta", "Number of hits per track VS #eta", 120, -3., 3., 55, 0, 55);
0149 
0150     hDof = ibooker.book1D("dof", "Number of Degree of Freedom", 55, 0, 55);
0151     hDofVsEta = ibooker.book2D("dofVsEta", "Number of Degree of Freedom VS #eta", 120, -3., 3., 55, 0, 55);
0152 
0153     hChi2Prob = ibooker.book1D("chi2Prob", "#chi^2 probability", 200, 0, 1);
0154     hChi2ProbVsEta = ibooker.book2D("chi2ProbVsEta", "#chi^2 probability VS #eta", 120, -3., 3., 200, 0, 1);
0155 
0156     hNumberOfTracks = ibooker.book1D("NumberOfTracks", "Number of reconstructed tracks per event", 200, 0, 200);
0157     hNumberOfTracksVsEta = ibooker.book2D(
0158         "NumberOfTracksVsEta", "Number of reconstructed tracks per event VS #eta", 120, -3., 3., 10, 0, 10);
0159 
0160     hChargeVsEta = ibooker.book2D("ChargeVsEta", "Charge vs #eta gen", 120, -3., 3., 4, -2., 2.);
0161     hChargeVsPt = ibooker.book2D("ChargeVsPt", "Charge vs P_{T} gen", 250, 0, 200, 4, -2., 2.);
0162     hPtRecVsPtGen = ibooker.book2D("PtRecVsPtGen", "P_{T} rec vs P_{T} gen", 250, 0, 200, 250, 0, 200);
0163 
0164     hDeltaPtVsEta = ibooker.book2D("DeltaPtVsEta", "#Delta P_{t} vs #eta gen", 120, -3., 3., 500, -250., 250.);
0165     hDeltaPt_In_Out_VsEta =
0166         ibooker.book2D("DeltaPt_In_Out_VsEta_", "P^{in}_{t} - P^{out}_{t} vs #eta gen", 120, -3., 3., 500, -250., 250.);
0167   }
0168 }
0169 
0170 void MuonTrackAnalyzer::analyze(const Event &event, const EventSetup &eventSetup) {
0171   LogDebug("MuonTrackAnalyzer") << "Run: " << event.id().run() << " Event: " << event.id().event();
0172 
0173   // Update the services
0174   theService->update(eventSetup);
0175 
0176   Handle<SimTrackContainer> simTracks;
0177   event.getByToken(theSimTracksToken, simTracks);
0178   fillPlots(event, simTracks);
0179 
0180   if (doTracksAnalysis)
0181     tracksAnalysis(event, eventSetup, simTracks);
0182 
0183   if (doSeedsAnalysis)
0184     seedsAnalysis(event, eventSetup, simTracks);
0185 }
0186 
0187 void MuonTrackAnalyzer::seedsAnalysis(const Event &event,
0188                                       const EventSetup &eventSetup,
0189                                       Handle<SimTrackContainer> simTracks) {
0190   MuonPatternRecoDumper debug;
0191 
0192   // Get the RecTrack collection from the event
0193   Handle<TrajectorySeedCollection> seeds;
0194   event.getByToken(theSeedsToken, seeds);
0195 
0196   LogTrace("MuonTrackAnalyzer") << "Number of reconstructed seeds: " << seeds->size() << endl;
0197 
0198   for (TrajectorySeedCollection::const_iterator seed = seeds->begin(); seed != seeds->end(); ++seed) {
0199     TrajectoryStateOnSurface seedTSOS = getSeedTSOS(*seed);
0200     pair<SimTrack, double> sim = getSimTrack(seedTSOS, simTracks);
0201     fillPlots(seedTSOS, sim.first, hRecoSeedInner, debug);
0202 
0203     std::pair<bool, FreeTrajectoryState> propSeed = theUpdator->propagateToNominalLine(seedTSOS);
0204     if (propSeed.first)
0205       fillPlots(propSeed.second, sim.first, hRecoSeedPCA, debug);
0206     else
0207       LogTrace("MuonTrackAnalyzer") << "Error in seed propagation" << endl;
0208   }
0209 }
0210 
0211 void MuonTrackAnalyzer::tracksAnalysis(const Event &event,
0212                                        const EventSetup &eventSetup,
0213                                        Handle<SimTrackContainer> simTracks) {
0214   MuonPatternRecoDumper debug;
0215 
0216   // Get the RecTrack collection from the event
0217   Handle<reco::TrackCollection> tracks;
0218   event.getByToken(theTracksToken, tracks);
0219 
0220   LogTrace("MuonTrackAnalyzer") << "Reconstructed tracks: " << tracks->size() << endl;
0221   hNumberOfTracks->Fill(tracks->size());
0222 
0223   if (!tracks->empty())
0224     numberOfRecTracks++;
0225 
0226   // Loop over the Rec tracks
0227   for (reco::TrackCollection::const_iterator t = tracks->begin(); t != tracks->end(); ++t) {
0228     reco::TransientTrack track(*t, &*theService->magneticField(), theService->trackingGeometry());
0229 
0230     TrajectoryStateOnSurface outerTSOS = track.outermostMeasurementState();
0231     TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
0232     TrajectoryStateOnSurface pcaTSOS = track.impactPointState();
0233 
0234     pair<SimTrack, double> sim = getSimTrack(pcaTSOS, simTracks);
0235     SimTrack simTrack = sim.first;
0236     hNumberOfTracksVsEta->Fill(simTrack.momentum().eta(), tracks->size());
0237     fillPlots(track, simTrack);
0238 
0239     LogTrace("MuonTrackAnalyzer") << "State at the outer surface: " << endl;
0240     fillPlots(outerTSOS, simTrack, hRecoTracksOuter, debug);
0241 
0242     LogTrace("MuonTrackAnalyzer") << "State at the inner surface: " << endl;
0243     fillPlots(innerTSOS, simTrack, hRecoTracksInner, debug);
0244 
0245     LogTrace("MuonTrackAnalyzer") << "State at PCA: " << endl;
0246     fillPlots(pcaTSOS, simTrack, hRecoTracksPCA, debug);
0247 
0248     double deltaPt_in_out = innerTSOS.globalMomentum().perp() - outerTSOS.globalMomentum().perp();
0249     hDeltaPt_In_Out_VsEta->Fill(simTrack.momentum().eta(), deltaPt_in_out);
0250 
0251     double deltaPt_pca_sim = pcaTSOS.globalMomentum().perp() - sqrt(simTrack.momentum().Perp2());
0252     hDeltaPtVsEta->Fill(simTrack.momentum().eta(), deltaPt_pca_sim);
0253 
0254     hChargeVsEta->Fill(simTrack.momentum().eta(), pcaTSOS.charge());
0255 
0256     hChargeVsPt->Fill(sqrt(simTrack.momentum().perp2()), pcaTSOS.charge());
0257 
0258     hPtRecVsPtGen->Fill(sqrt(simTrack.momentum().perp2()), pcaTSOS.globalMomentum().perp());
0259   }
0260   LogTrace("MuonTrackAnalyzer") << "--------------------------------------------" << endl;
0261 }
0262 
0263 void MuonTrackAnalyzer::fillPlots(const Event &event, edm::Handle<edm::SimTrackContainer> &simTracks) {
0264   if (!checkMuonSimHitPresence(event, simTracks))
0265     return;
0266 
0267   // Loop over the Sim tracks
0268   SimTrackContainer::const_iterator simTrack;
0269   LogTrace("MuonTrackAnalyzer") << "Simulated tracks: " << simTracks->size() << endl;
0270 
0271   for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack)
0272     if (abs((*simTrack).type()) == 13) {
0273       if (!isInTheAcceptance((*simTrack).momentum().eta()))
0274         continue;  // FIXME!!
0275 
0276       numberOfSimTracks++;
0277 
0278       LogTrace("MuonTrackAnalyzer") << "Simualted muon:" << endl;
0279       LogTrace("MuonTrackAnalyzer") << "Sim pT: " << sqrt((*simTrack).momentum().perp2()) << endl;
0280       LogTrace("MuonTrackAnalyzer") << "Sim Eta: " << (*simTrack).momentum().eta() << endl;  // FIXME
0281 
0282       hSimTracks->Fill((*simTrack).momentum().mag(),
0283                        sqrt((*simTrack).momentum().perp2()),
0284                        (*simTrack).momentum().eta(),
0285                        (*simTrack).momentum().phi(),
0286                        -(*simTrack).type() / abs((*simTrack).type()));  // Double FIXME
0287       LogTrace("MuonTrackAnalyzer") << "hSimTracks filled" << endl;
0288     }
0289 
0290   LogTrace("MuonTrackAnalyzer") << endl;
0291 }
0292 
0293 void MuonTrackAnalyzer::fillPlots(reco::TransientTrack &track, SimTrack &simTrack) {
0294   LogTrace("MuonTrackAnalyzer") << "Analizer: New track, chi2: " << track.chi2() << " dof: " << track.ndof() << endl;
0295   hChi2->Fill(track.chi2());
0296   hDof->Fill(track.ndof());
0297   hChi2Norm->Fill(track.normalizedChi2());
0298   hHitsPerTrack->Fill(track.recHitsSize());
0299 
0300   hChi2Prob->Fill(ChiSquaredProbability(track.chi2(), track.ndof()));
0301 
0302   hChi2VsEta->Fill(simTrack.momentum().eta(), track.chi2());
0303   hChi2NormVsEta->Fill(simTrack.momentum().eta(), track.normalizedChi2());
0304   hChi2ProbVsEta->Fill(simTrack.momentum().eta(), ChiSquaredProbability(track.chi2(), track.ndof()));
0305   hHitsPerTrackVsEta->Fill(simTrack.momentum().eta(), track.recHitsSize());
0306   hDofVsEta->Fill(simTrack.momentum().eta(), track.ndof());
0307 }
0308 
0309 void MuonTrackAnalyzer::fillPlots(TrajectoryStateOnSurface &recoTSOS,
0310                                   SimTrack &simTrack,
0311                                   HTrack *histo,
0312                                   MuonPatternRecoDumper &debug) {
0313   LogTrace("MuonTrackAnalyzer") << debug.dumpTSOS(recoTSOS) << endl;
0314   histo->Fill(recoTSOS);
0315 
0316   GlobalVector tsosVect = recoTSOS.globalMomentum();
0317   math::XYZVectorD reco(tsosVect.x(), tsosVect.y(), tsosVect.z());
0318   double deltaRVal = deltaR<double>(reco.eta(), reco.phi(), simTrack.momentum().eta(), simTrack.momentum().phi());
0319   histo->FillDeltaR(deltaRVal);
0320 
0321   histo->computeResolutionAndPull(recoTSOS, simTrack);
0322 }
0323 
0324 void MuonTrackAnalyzer::fillPlots(FreeTrajectoryState &recoFTS,
0325                                   SimTrack &simTrack,
0326                                   HTrack *histo,
0327                                   MuonPatternRecoDumper &debug) {
0328   LogTrace("MuonTrackAnalyzer") << debug.dumpFTS(recoFTS) << endl;
0329   histo->Fill(recoFTS);
0330 
0331   GlobalVector ftsVect = recoFTS.momentum();
0332   math::XYZVectorD reco(ftsVect.x(), ftsVect.y(), ftsVect.z());
0333   double deltaRVal = deltaR<double>(reco.eta(), reco.phi(), simTrack.momentum().eta(), simTrack.momentum().phi());
0334   histo->FillDeltaR(deltaRVal);
0335 
0336   histo->computeResolutionAndPull(recoFTS, simTrack);
0337 }
0338 
0339 pair<SimTrack, double> MuonTrackAnalyzer::getSimTrack(TrajectoryStateOnSurface &tsos,
0340                                                       Handle<SimTrackContainer> simTracks) {
0341   //   // Loop over the Sim tracks
0342   //   SimTrackContainer::const_iterator simTrack;
0343 
0344   //   SimTrack result;
0345   //   int mu=0;
0346   //   for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack)
0347   //     if (abs((*simTrack).type()) == 13) {
0348   //       result = *simTrack;
0349   //       ++mu;
0350   //     }
0351 
0352   //   if(mu != 1) LogTrace("MuonTrackAnalyzer") << "WARNING!! more than 1 simulated muon!!" <<endl;
0353   //   return result;
0354 
0355   // Loop over the Sim tracks
0356   SimTrackContainer::const_iterator simTrack;
0357 
0358   SimTrack result;
0359 
0360   double bestDeltaR = 10e5;
0361   for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0362     if (abs((*simTrack).type()) != 13)
0363       continue;
0364 
0365     //    double newDeltaR = tsos.globalMomentum().basicVector().deltaR(simTrack->momentum().vect());
0366     GlobalVector tsosVect = tsos.globalMomentum();
0367     math::XYZVectorD vect(tsosVect.x(), tsosVect.y(), tsosVect.z());
0368     double newDeltaR = deltaR<double>(vect.eta(), vect.phi(), simTrack->momentum().eta(), simTrack->momentum().phi());
0369 
0370     if (newDeltaR < bestDeltaR) {
0371       LogTrace("MuonTrackAnalyzer") << "Matching Track with DeltaR = " << newDeltaR << endl;
0372       bestDeltaR = newDeltaR;
0373       result = *simTrack;
0374     }
0375   }
0376   return pair<SimTrack, double>(result, bestDeltaR);
0377 }
0378 
0379 bool MuonTrackAnalyzer::isInTheAcceptance(double eta) {
0380   switch (theEtaRange) {
0381     case all:
0382       return (abs(eta) <= 2.4) ? true : false;
0383     case barrel:
0384       return (abs(eta) < 1.1) ? true : false;
0385     case endcap:
0386       return (abs(eta) >= 1.1 && abs(eta) <= 2.4) ? true : false;
0387     default: {
0388       LogTrace("MuonTrackAnalyzer") << "No correct Eta range selected!! " << endl;
0389       return false;
0390     }
0391   }
0392 }
0393 
0394 bool MuonTrackAnalyzer::checkMuonSimHitPresence(const Event &event, edm::Handle<edm::SimTrackContainer> simTracks) {
0395   // Get the SimHit collection from the event
0396   Handle<PSimHitContainer> dtSimHits;
0397   event.getByToken(theDTSimHitToken, dtSimHits);
0398 
0399   Handle<PSimHitContainer> cscSimHits;
0400   event.getByToken(theCSCSimHitToken, cscSimHits);
0401 
0402   Handle<PSimHitContainer> rpcSimHits;
0403   event.getByToken(theRPCSimHitToken, rpcSimHits);
0404 
0405   map<unsigned int, vector<const PSimHit *> > mapOfMuonSimHits;
0406 
0407   for (PSimHitContainer::const_iterator simhit = dtSimHits->begin(); simhit != dtSimHits->end(); ++simhit) {
0408     if (abs(simhit->particleType()) != 13)
0409       continue;
0410     mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
0411   }
0412 
0413   for (PSimHitContainer::const_iterator simhit = cscSimHits->begin(); simhit != cscSimHits->end(); ++simhit) {
0414     if (abs(simhit->particleType()) != 13)
0415       continue;
0416     mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
0417   }
0418 
0419   for (PSimHitContainer::const_iterator simhit = rpcSimHits->begin(); simhit != rpcSimHits->end(); ++simhit) {
0420     if (abs(simhit->particleType()) != 13)
0421       continue;
0422     mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
0423   }
0424 
0425   bool presence = false;
0426 
0427   for (SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0428     if (abs(simTrack->type()) != 13)
0429       continue;
0430 
0431     map<unsigned int, vector<const PSimHit *> >::const_iterator mapIterator =
0432         mapOfMuonSimHits.find(simTrack->trackId());
0433 
0434     if (mapIterator != mapOfMuonSimHits.end())
0435       presence = true;
0436   }
0437 
0438   return presence;
0439 }
0440 
0441 TrajectoryStateOnSurface MuonTrackAnalyzer::getSeedTSOS(const TrajectorySeed &seed) {
0442   // Get the Trajectory State on Det (persistent version of a TSOS) from the seed
0443   PTrajectoryStateOnDet pTSOD = seed.startingState();
0444 
0445   // Transform it in a TrajectoryStateOnSurface
0446 
0447   DetId seedDetId(pTSOD.detId());
0448 
0449   const GeomDet *gdet = theService->trackingGeometry()->idToDet(seedDetId);
0450 
0451   TrajectoryStateOnSurface initialState =
0452       trajectoryStateTransform::transientState(pTSOD, &(gdet->surface()), &*theService->magneticField());
0453 
0454   // Get the layer on which the seed relies
0455   const DetLayer *initialLayer = theService->detLayerGeometry()->idToLayer(seedDetId);
0456 
0457   PropagationDirection detLayerOrder = oppositeToMomentum;
0458 
0459   // ask for compatible layers
0460   vector<const DetLayer *> detLayers;
0461   detLayers =
0462       theService->muonNavigationSchool()->compatibleLayers(*initialLayer, *initialState.freeState(), detLayerOrder);
0463 
0464   TrajectoryStateOnSurface result = initialState;
0465   if (!detLayers.empty()) {
0466     const DetLayer *finalLayer = detLayers.back();
0467     const TrajectoryStateOnSurface propagatedState =
0468         theService->propagator(theSeedPropagatorName)->propagate(initialState, finalLayer->surface());
0469     if (propagatedState.isValid())
0470       result = propagatedState;
0471   }
0472 
0473   return result;
0474 }