TrackEfficiencyMonitor.cc

Go to the documentation of this file.

/*
 *  See header file for a description of this class.
 *
 *  \author Jeremy Andrea
 */

#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
//#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "MagneticField/Engine/interface/MagneticField.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "DQMServices/Core/interface/DQMStore.h"
#include "DQM/TrackingMonitor/interface/TrackEfficiencyMonitor.h"
#include <string>

// needed to compute the efficiency

#include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
#include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"

#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
#include "TrackingTools/GeomPropagators/interface/SmartPropagator.h"
#include "DataFormats/Math/interface/Vector3D.h"
#include "DataFormats/GeometrySurface/interface/Plane.h"
#include "DataFormats/GeometrySurface/interface/BoundPlane.h"
#include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"

#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "DataFormats/GeometrySurface/interface/Cylinder.h"

#include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
#include <RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h>

#include "RecoTracker/Record/interface/CkfComponentsRecord.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "RecoTracker/Record/interface/NavigationSchoolRecord.h"



//-----------------------------------------------------------------------------------
TrackEfficiencyMonitor::TrackEfficiencyMonitor(const edm::ParameterSet& iConfig) 
//-----------------------------------------------------------------------------------
{
  dqmStore_            = edm::Service<DQMStore>().operator->(); 
 
  theRadius_           = iConfig.getParameter<double>("theRadius");
  theMaxZ_             = iConfig.getParameter<double>("theMaxZ");
  isBFieldOff_         = iConfig.getParameter<bool>("isBFieldOff");
  trackEfficiency_     = iConfig.getParameter<bool>("trackEfficiency");
  theTKTracksLabel_    = iConfig.getParameter<edm::InputTag>("TKTrackCollection");
  theSTATracksLabel_   = iConfig.getParameter<edm::InputTag>("STATrackCollection");
  
  theTKTracksToken_  = consumes<reco::TrackCollection>(theTKTracksLabel_);
  theSTATracksToken_ = consumes<reco::TrackCollection>(theSTATracksLabel_);

  conf_ = iConfig;
}



//-----------------------------------------------------------------------------------
TrackEfficiencyMonitor::~TrackEfficiencyMonitor() 
//-----------------------------------------------------------------------------------
{}

//-----------------------------------------------------------------------------------
void TrackEfficiencyMonitor::bookHistograms(DQMStore::IBooker & ibooker,
                        edm::Run const & /* iRun */,
                        edm::EventSetup const & /* iSetup */)
//-----------------------------------------------------------------------------------
{
  std::string MEFolderName = conf_.getParameter<std::string>("FolderName"); 
  std::string AlgoName     = conf_.getParameter<std::string>("AlgoName");
  
  
  ibooker.setCurrentFolder(MEFolderName);
  
  //
  int    muonXBin = conf_.getParameter<int>   ("muonXBin");
  double muonXMin = conf_.getParameter<double>("muonXMin");
  double muonXMax = conf_.getParameter<double>("muonXMax");
 
  histname = "muonX_";
  muonX = ibooker.book1D(histname+AlgoName, histname+AlgoName, muonXBin, muonXMin, muonXMax);
  muonX->setAxisTitle("");
  
  //
  int    muonYBin = conf_.getParameter<int>   ("muonYBin");
  double muonYMin = conf_.getParameter<double>("muonYMin");
  double muonYMax = conf_.getParameter<double>("muonYMax");
 
  histname = "muonY_";
  muonY = ibooker.book1D(histname+AlgoName, histname+AlgoName, muonYBin, muonYMin, muonYMax);
  muonY->setAxisTitle("");
  
  //
  int    muonZBin = conf_.getParameter<int>   ("muonZBin");
  double muonZMin = conf_.getParameter<double>("muonZMin");
  double muonZMax = conf_.getParameter<double>("muonZMax");
 
  histname = "muonZ_";
  muonZ = ibooker.book1D(histname+AlgoName, histname+AlgoName, muonZBin, muonZMin, muonZMax);
  muonZ->setAxisTitle("");
  
  //
  int    muonEtaBin = conf_.getParameter<int>   ("muonEtaBin");
  double muonEtaMin = conf_.getParameter<double>("muonEtaMin");
  double muonEtaMax = conf_.getParameter<double>("muonEtaMax");
 
  histname = "muonEta_";
  muonEta = ibooker.book1D(histname+AlgoName, histname+AlgoName, muonEtaBin, muonEtaMin, muonEtaMax);
  muonEta->setAxisTitle("");
  
  //
  int    muonPhiBin = conf_.getParameter<int>   ("muonPhiBin");
  double muonPhiMin = conf_.getParameter<double>("muonPhiMin");
  double muonPhiMax = conf_.getParameter<double>("muonPhiMax");
 
  histname = "muonPhi_";
  muonPhi = ibooker.book1D(histname+AlgoName, histname+AlgoName, muonPhiBin, muonPhiMin, muonPhiMax);
  muonPhi->setAxisTitle("");
  
  //
  int    muonD0Bin = conf_.getParameter<int>   ("muonD0Bin");
  double muonD0Min = conf_.getParameter<double>("muonD0Min");
  double muonD0Max = conf_.getParameter<double>("muonD0Max");
 
  histname = "muonD0_";
  muonD0 = ibooker.book1D(histname+AlgoName, histname+AlgoName, muonD0Bin, muonD0Min, muonD0Max);
  muonD0->setAxisTitle("");
  
  //
  int    muonCompatibleLayersBin = conf_.getParameter<int>   ("muonCompatibleLayersBin");
  double muonCompatibleLayersMin = conf_.getParameter<double>("muonCompatibleLayersMin");
  double muonCompatibleLayersMax = conf_.getParameter<double>("muonCompatibleLayersMax");
 
  histname = "muonCompatibleLayers_";
  muonCompatibleLayers = ibooker.book1D(histname+AlgoName, histname+AlgoName, muonCompatibleLayersBin, muonCompatibleLayersMin, muonCompatibleLayersMax);
  muonCompatibleLayers->setAxisTitle("");

  //------------------------------------------------------------------------------------
  
  //
  int    trackXBin = conf_.getParameter<int>   ("trackXBin");
  double trackXMin = conf_.getParameter<double>("trackXMin");
  double trackXMax = conf_.getParameter<double>("trackXMax");
 
  histname = "trackX_";
  trackX = ibooker.book1D(histname+AlgoName, histname+AlgoName, trackXBin, trackXMin, trackXMax);
  trackX->setAxisTitle("");
  
  //
  int    trackYBin = conf_.getParameter<int>   ("trackYBin");
  double trackYMin = conf_.getParameter<double>("trackYMin");
  double trackYMax = conf_.getParameter<double>("trackYMax");
 
  histname = "trackY_";
  trackY = ibooker.book1D(histname+AlgoName, histname+AlgoName, trackYBin, trackYMin, trackYMax);
  trackY->setAxisTitle("");
  
  //
  int    trackZBin = conf_.getParameter<int>   ("trackZBin");
  double trackZMin = conf_.getParameter<double>("trackZMin");
  double trackZMax = conf_.getParameter<double>("trackZMax");
 
  histname = "trackZ_";
  trackZ = ibooker.book1D(histname+AlgoName, histname+AlgoName, trackZBin, trackZMin, trackZMax);
  trackZ->setAxisTitle("");
  
  //
  int    trackEtaBin = conf_.getParameter<int>   ("trackEtaBin");
  double trackEtaMin = conf_.getParameter<double>("trackEtaMin");
  double trackEtaMax = conf_.getParameter<double>("trackEtaMax");
 
  histname = "trackEta_";
  trackEta = ibooker.book1D(histname+AlgoName, histname+AlgoName, trackEtaBin, trackEtaMin, trackEtaMax);
  trackEta->setAxisTitle("");
  
  //
  int    trackPhiBin = conf_.getParameter<int>   ("trackPhiBin");
  double trackPhiMin = conf_.getParameter<double>("trackPhiMin");
  double trackPhiMax = conf_.getParameter<double>("trackPhiMax");
 
  histname = "trackPhi_";
  trackPhi = ibooker.book1D(histname+AlgoName, histname+AlgoName, trackPhiBin, trackPhiMin, trackPhiMax);
  trackPhi->setAxisTitle("");
  
  //
  int    trackD0Bin = conf_.getParameter<int>   ("trackD0Bin");
  double trackD0Min = conf_.getParameter<double>("trackD0Min");
  double trackD0Max = conf_.getParameter<double>("trackD0Max");
 
  histname = "trackD0_";
  trackD0 = ibooker.book1D(histname+AlgoName, histname+AlgoName, trackD0Bin, trackD0Min, trackD0Max);
  trackD0->setAxisTitle("");
  
  //
  int    trackCompatibleLayersBin = conf_.getParameter<int>   ("trackCompatibleLayersBin");
  double trackCompatibleLayersMin = conf_.getParameter<double>("trackCompatibleLayersMin");
  double trackCompatibleLayersMax = conf_.getParameter<double>("trackCompatibleLayersMax");
 
  histname = "trackCompatibleLayers_";
  trackCompatibleLayers = ibooker.book1D(histname+AlgoName, histname+AlgoName, trackCompatibleLayersBin, trackCompatibleLayersMin, trackCompatibleLayersMax);
  trackCompatibleLayers->setAxisTitle("");

  //------------------------------------------------------------------------------------
 
  //
  int    deltaXBin = conf_.getParameter<int>   ("deltaXBin");
  double deltaXMin = conf_.getParameter<double>("deltaXMin");
  double deltaXMax = conf_.getParameter<double>("deltaXMax");
 
  histname = "deltaX_";
  deltaX = ibooker.book1D(histname+AlgoName, histname+AlgoName, deltaXBin, deltaXMin, deltaXMax);
  deltaX->setAxisTitle("");
  
  //
  int    deltaYBin = conf_.getParameter<int>   ("deltaYBin");
  double deltaYMin = conf_.getParameter<double>("deltaYMin");
  double deltaYMax = conf_.getParameter<double>("deltaYMax");
 
  histname = "deltaY_";
  deltaY = ibooker.book1D(histname+AlgoName, histname+AlgoName, deltaYBin, deltaYMin, deltaYMax);
  deltaY->setAxisTitle("");
  
  //
  int    signDeltaXBin = conf_.getParameter<int>   ("signDeltaXBin");
  double signDeltaXMin = conf_.getParameter<double>("signDeltaXMin");
  double signDeltaXMax = conf_.getParameter<double>("signDeltaXMax");
 
  histname = "signDeltaX_";
  signDeltaX = ibooker.book1D(histname+AlgoName, histname+AlgoName, signDeltaXBin, signDeltaXMin, signDeltaXMax);
  signDeltaX->setAxisTitle("");
  
  //
  int    signDeltaYBin = conf_.getParameter<int>   ("signDeltaYBin");
  double signDeltaYMin = conf_.getParameter<double>("signDeltaYMin");
  double signDeltaYMax = conf_.getParameter<double>("signDeltaYMax");
 
  histname = "signDeltaY_";
  signDeltaY = ibooker.book1D(histname+AlgoName, histname+AlgoName, signDeltaYBin, signDeltaYMin, signDeltaYMax);
  signDeltaY->setAxisTitle("");

}


//-----------------------------------------------------------------------------------
void TrackEfficiencyMonitor::beginJob(void) 
//-----------------------------------------------------------------------------------
{
  
}



//-----------------------------------------------------------------------------------
void TrackEfficiencyMonitor::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) 
//-----------------------------------------------------------------------------------
{
   

  edm::Handle<reco::TrackCollection> tkTracks;
  iEvent.getByToken(theTKTracksToken_, tkTracks);
  edm::Handle<reco::TrackCollection> staTracks;
  iEvent.getByToken(theSTATracksToken_, staTracks);
  edm::ESHandle<NavigationSchool> nav;
  iSetup.get<NavigationSchoolRecord>().get("CosmicNavigationSchool", nav); 
  iSetup.get<CkfComponentsRecord>().get(measurementTrackerHandle);
 
  nCompatibleLayers = 0; 
  
  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theTTrackBuilder);
  iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny",thePropagator);
  iSetup.get<IdealMagneticFieldRecord>().get(bField); 
  iSetup.get<TrackerRecoGeometryRecord>().get(theTracker);
  theNavigation = new DirectTrackerNavigation(theTracker);
  
  
  if(trackEfficiency_){
  //---------------------------------------------------
  // Select muons with good quality
  // If B field is on, no up-down matching between the muons
  //---------------------------------------------------  
  bool isGoodMuon = false;
  double mudd0 = 0., mudphi = 0., muddsz = 0., mudeta = 0.;
  if(isBFieldOff_)
  {
    if ( staTracks->size() == 2 )
    {
      for ( unsigned int bindex = 0; bindex < staTracks->size(); ++bindex ) 
      {
    
    if (0 == bindex){
       mudd0+=(*staTracks)[bindex].d0(); 
       mudphi+=(*staTracks)[bindex].phi();
       muddsz+=(*staTracks)[bindex].dsz(); 
       mudeta+=(*staTracks)[bindex].eta();}
    if (1 == bindex){
       mudd0-=(*staTracks)[bindex].d0(); 
       mudphi-=(*staTracks)[bindex].phi();
       muddsz-=(*staTracks)[bindex].dsz(); 
       mudeta-=(*staTracks)[bindex].eta();}
       }
      if ((fabs(mudd0)<15.0)&&(fabs(mudphi)<0.045)&&(fabs(muddsz)<20.0)&&(fabs(mudeta)<0.060))  isGoodMuon = true;
     }
     
    if(isGoodMuon) testTrackerTracks(tkTracks,staTracks, *nav.product());
  
  }
  else if ( staTracks->size() == 1 || staTracks->size() == 2) testTrackerTracks(tkTracks,staTracks, *nav.product());
  }
  
  
  if(!trackEfficiency_ && tkTracks->size() == 1 ){
    if( (tkTracks->front()).normalizedChi2() < 5 &&
    (tkTracks->front()).hitPattern().numberOfValidHits() > 8)
     testSTATracks(tkTracks,staTracks);  
  }
  
  
  delete theNavigation;
   
}



//-----------------------------------------------------------------------------------
void TrackEfficiencyMonitor::endJob(void) 
//-----------------------------------------------------------------------------------
{ 
  bool outputMEsInRootFile = conf_.getParameter<bool>("OutputMEsInRootFile");
  std::string outputFileName = conf_.getParameter<std::string>("OutputFileName");
  if(outputMEsInRootFile){
    dqmStore_->showDirStructure(); 
    dqmStore_->save(outputFileName); 
    }

  //if ( theNavigation ) delete theNavigation; 

}



//-----------------------------------------------------------------------------------
void TrackEfficiencyMonitor::testTrackerTracks(edm::Handle<reco::TrackCollection> tkTracks, edm::Handle<reco::TrackCollection> staTracks, const NavigationSchool& navigationSchool)
//-----------------------------------------------------------------------------------
{
  
  const std::string metname = "testStandAloneMuonTracks";

  //---------------------------------------------------
  // get the index of the "up" muon
  // histograms will only be computed for the "up" muon
  //---------------------------------------------------
 
  int nUpMuon = 0;
  int idxUpMuon = -1;
  for(unsigned int i =0; i< staTracks->size(); i++){
     if( checkSemiCylinder((*staTracks)[i]) == TrackEfficiencyMonitor::Up ){
        nUpMuon++;
    idxUpMuon = i;
     }
  }
    
   
  if( nUpMuon == 1)
  {
    
    //---------------------------------------------------
    //get the muon informations
    //---------------------------------------------------
    
    reco::TransientTrack staTT = theTTrackBuilder->build((*staTracks)[idxUpMuon]);
    double ipX   = staTT.impactPointState().globalPosition().x();
    double ipY   = staTT.impactPointState().globalPosition().y();
    double ipZ   = staTT.impactPointState().globalPosition().z();
    double eta   = staTT.impactPointState().globalDirection().eta();
    double phi   = staTT.impactPointState().globalDirection().phi();
    double d0    = (*staTracks)[idxUpMuon].d0();
    
    TrajectoryStateOnSurface theTSOS = staTT.outermostMeasurementState(); 
    TrajectoryStateOnSurface theTSOSCompLayers = staTT.outermostMeasurementState(); 
    
    
    //---------------------------------------------------
    //check if the muon is in the tracker acceptance
    //---------------------------------------------------    
    bool isInTrackerAcceptance = false;   
    isInTrackerAcceptance = trackerAcceptance( theTSOS, theRadius_ , theMaxZ_ );
    
    //---------------------------------------------------st
    //count the number of compatible layers
    //---------------------------------------------------       
    nCompatibleLayers = compatibleLayers(navigationSchool, theTSOSCompLayers);
    
    if(isInTrackerAcceptance && (*staTracks)[idxUpMuon].hitPattern().numberOfValidHits() > 28)
    {
      //---------------------------------------------------
      //count the number of good muon candidates
      //---------------------------------------------------   
        
      TrajectoryStateOnSurface staState;      
      LocalVector diffLocal;   
     
     
      bool isTrack = false;
      if(!tkTracks->empty())
      {
       //---------------------------------------------------   
       //look for associated tracks
       //---------------------------------------------------   
        float DR2min = 1000;
        reco::TrackCollection::const_iterator closestTrk = tkTracks->end();

       for(reco::TrackCollection::const_iterator tkTrack = tkTracks->begin(); tkTrack != tkTracks->end(); ++tkTrack)
       {
          reco::TransientTrack tkTT = theTTrackBuilder->build(*tkTrack);
          TrajectoryStateOnSurface tkInner = tkTT.innermostMeasurementState();
          staState = thePropagator->propagate(staTT.outermostMeasurementState(),tkInner.surface());  
      failedToPropagate = 1;
            
      if(staState.isValid())
      {
            failedToPropagate = 0;
            diffLocal = tkInner.localPosition() - staState.localPosition();
        double DR2 = diffLocal.x()*diffLocal.x()+diffLocal.y()*diffLocal.y();
        if (DR2<DR2min) { DR2min = DR2; closestTrk = tkTrack;}
        if ( pow(DR2,0.5) < 100. ) isTrack = true;
            }
          }
     
      
      if (DR2min != 1000) 
      { 
        
        reco::TransientTrack tkTT = theTTrackBuilder->build(*closestTrk);
            TrajectoryStateOnSurface tkInner = tkTT.innermostMeasurementState();
        staState = thePropagator->propagate(staTT.outermostMeasurementState(),tkInner.surface());
        deltaX->Fill(diffLocal.x());
            deltaY->Fill(diffLocal.y());
        signDeltaX->Fill(diffLocal.x()/(tkInner.localError().positionError().xx() + staState.localError().positionError().xx()));
        signDeltaY->Fill(diffLocal.y()/(tkInner.localError().positionError().yy() + staState.localError().positionError().yy()));
               
       } 
        }
      
      
      if(failedToPropagate == 0)
      {       
     muonX->Fill(ipX);
         muonY->Fill(ipY);
         muonZ->Fill(ipZ);
         muonEta->Fill(eta);
         muonPhi->Fill(phi);
     muonD0->Fill(d0);
         muonCompatibleLayers->Fill(nCompatibleLayers);

         if(isTrack)
         { 
            trackX->Fill(ipX);
            trackY->Fill(ipY);
            trackZ->Fill(ipZ);
            trackEta->Fill(eta);
            trackPhi->Fill(phi);
        trackD0->Fill(d0);
        trackCompatibleLayers->Fill(nCompatibleLayers);
          } 
        }
     }
   }
    
}











//-----------------------------------------------------------------------------------
void TrackEfficiencyMonitor::testSTATracks(edm::Handle<reco::TrackCollection> tkTracks, edm::Handle<reco::TrackCollection> staTracks)
//-----------------------------------------------------------------------------------
{

   
   reco::TransientTrack tkTT = theTTrackBuilder->build(tkTracks->front());
   double ipX   = tkTT.impactPointState().globalPosition().x();
   double ipY   = tkTT.impactPointState().globalPosition().y();
   double ipZ   = tkTT.impactPointState().globalPosition().z();
   double eta   = tkTT.impactPointState().globalDirection().eta();
   double phi   = tkTT.impactPointState().globalDirection().phi();
   double d0    = (*tkTracks)[0].d0();

   TrajectoryStateOnSurface tkInner = tkTT.innermostMeasurementState();
   LocalVector diffLocal;   
   TrajectoryStateOnSurface staState;      
   bool isTrack = false;
   
   
   if(!staTracks->empty()){     
    
       //---------------------------------------------------   
       //look for associated muons
       //---------------------------------------------------   
       
    float DR2min = 1000;
    reco::TrackCollection::const_iterator closestTrk = staTracks->end();
       //----------------------loop on tracker tracks:
    for(reco::TrackCollection::const_iterator staTrack = staTracks->begin(); staTrack != staTracks->end(); ++staTrack){
      if(checkSemiCylinder(*staTrack) == TrackEfficiencyMonitor::Up){
      
        reco::TransientTrack staTT = theTTrackBuilder->build(*staTrack); 
        failedToPropagate = 1;  
    staState = thePropagator->propagate(staTT.outermostMeasurementState(),tkInner.surface());   

       
      if(staState.isValid())
      {
            failedToPropagate = 0;
            diffLocal = tkInner.localPosition() - staState.localPosition();
               
        double DR2 = diffLocal.x()*diffLocal.x()+diffLocal.y()*diffLocal.y();
        if (DR2<DR2min) { DR2min = DR2; closestTrk = staTrack;}
        if ( pow(DR2,0.5) < 100. ) isTrack = true;
           }
          }     
      }
    }
   
   if(failedToPropagate == 0)
      {         
     
         trackX->Fill(ipX);
         trackY->Fill(ipY);
         trackZ->Fill(ipZ);
         trackEta->Fill(eta);
         trackPhi->Fill(phi);
     trackD0->Fill(d0);

         if(isTrack)
         { 
        muonX->Fill(ipX);
            muonY->Fill(ipY);
            muonZ->Fill(ipZ);
            muonEta->Fill(eta);
            muonPhi->Fill(phi);
        muonD0->Fill(d0);
          } 
        }



}

















//-----------------------------------------------------------------------------------
TrackEfficiencyMonitor::SemiCylinder TrackEfficiencyMonitor::checkSemiCylinder(const reco::Track& tk)
//-----------------------------------------------------------------------------------
{
  return tk.innerPosition().phi() > 0 ? TrackEfficiencyMonitor::Up : TrackEfficiencyMonitor::Down;
}

  
  
//-----------------------------------------------------------------------------------
bool TrackEfficiencyMonitor::trackerAcceptance( TrajectoryStateOnSurface theTSOS, double theRadius, double theMaxZ)
//-----------------------------------------------------------------------------------
{

  //---------------------------------------------------  
  // check if the muon is in the tracker acceptance    
  // check the compatibility with a cylinder of radius "theRadius"
  //---------------------------------------------------   
      
  //Propagator*  theTmpPropagator = new SteppingHelixPropagator(&*bField,anyDirection);
  
  //Propagator*  theTmpPropagator = &*thePropagator;
  Propagator* theTmpPropagator = &*thePropagator->clone();
  
  if (theTSOS.globalPosition().y() <0 ) theTmpPropagator->setPropagationDirection(oppositeToMomentum);
  else                                  theTmpPropagator->setPropagationDirection(alongMomentum);
  
  Cylinder::PositionType pos0;
  Cylinder::RotationType rot0;
  const Cylinder::CylinderPointer cyl = Cylinder::build(pos0, rot0, theRadius);
  TrajectoryStateOnSurface tsosAtCyl  = theTmpPropagator->propagate(*theTSOS.freeState(), *cyl);  
  double accept = false;
  if ( tsosAtCyl.isValid() ) {
    if (fabs(tsosAtCyl.globalPosition().z())<theMaxZ ) {
      Cylinder::PositionType pos02;
      Cylinder::RotationType rot02;
      const Cylinder::CylinderPointer cyl2 = Cylinder::build(pos02, rot02, theRadius -10);
      TrajectoryStateOnSurface tsosAtCyl2 = theTmpPropagator->propagate(*tsosAtCyl.freeState(), *cyl2);
      if ( tsosAtCyl2.isValid() ) {
        Cylinder::PositionType pos03;
        Cylinder::RotationType rot03;
        const Cylinder::CylinderPointer cyl3 = Cylinder::build(pos03, rot03, theRadius );
        TrajectoryStateOnSurface tsosAtCyl3 = theTmpPropagator->propagate(*tsosAtCyl2.freeState(), *cyl3);
    if ( tsosAtCyl3.isValid() ) {
      accept=true;
        }
      }
      
    }
  }
  delete theTmpPropagator;
  //muon propagated to the barrel cylinder
  return accept;

}



//-----------------------------------------------------------------------------------
int TrackEfficiencyMonitor::compatibleLayers(const NavigationSchool& navigationSchool,  TrajectoryStateOnSurface theTSOS)
//-----------------------------------------------------------------------------------
{  
  //---------------------------------------------------   
  // check the number of compatible layers
  //---------------------------------------------------   
  
  std::vector< const BarrelDetLayer*> barrelTOBLayers = measurementTrackerHandle->geometricSearchTracker()->tobLayers() ;
  
  unsigned int layers = 0;
  for ( unsigned int k=0 ; k < barrelTOBLayers.size() ; k++ ) 
  {  
    const DetLayer* firstLay = barrelTOBLayers[barrelTOBLayers.size() -1 - k];  
    
    //Propagator*  theTmpPropagator = new SteppingHelixPropagator(&*bField,anyDirection);
    
    
    
    Propagator* theTmpPropagator = &*thePropagator->clone();
    theTmpPropagator->setPropagationDirection(alongMomentum); 
    
    TrajectoryStateOnSurface startTSOS = theTmpPropagator->propagate(*theTSOS.freeState(),firstLay->surface());
    
    std::vector< const DetLayer*> trackCompatibleLayers;
    
    findDetLayer      = true;
    bool isUpMuon     = false;
    bool firstdtep    = true; 
    
    if(startTSOS.isValid())
    {
      
      if(firstdtep) layers++; 
      
      int nwhile = 0;
      
      //for other compatible layers 
      while ( startTSOS.isValid() &&  firstLay && findDetLayer)
      {
    
    if(firstdtep && startTSOS.globalPosition().y()> 0) isUpMuon = true;
    
    if(firstdtep){
      std::vector< const DetLayer*> firstCompatibleLayers;
      firstCompatibleLayers.push_back(firstLay);
      std::pair<TrajectoryStateOnSurface, const  DetLayer* > nextLayer = findNextLayer(theTSOS, firstCompatibleLayers, isUpMuon );
      firstdtep = false;
    }
    else{
      trackCompatibleLayers = navigationSchool.nextLayers(*firstLay, *(startTSOS.freeState()),alongMomentum);
          if (!trackCompatibleLayers.empty() ){ 
        std::pair<TrajectoryStateOnSurface, const  DetLayer* > nextLayer = findNextLayer(startTSOS, trackCompatibleLayers, isUpMuon );
        if (firstLay != nextLayer.second ){
          firstLay  = nextLayer.second;
          startTSOS = nextLayer.first;
          layers++;         
        }
        else firstLay=nullptr;
          }
    }
    nwhile++;
    if(nwhile > 100) break;
        
      }
      delete theTmpPropagator;
      break;
    }
    delete theTmpPropagator;
  }
  return layers;
  
}




//-----------------------------------------------------------------------------------
std::pair<TrajectoryStateOnSurface,  const DetLayer*>  TrackEfficiencyMonitor::findNextLayer( TrajectoryStateOnSurface sTSOS, const std::vector< const DetLayer*>& trackCompatibleLayers , bool isUpMuon )
//-----------------------------------------------------------------------------------
{
  
  //Propagator*  theTmpPropagator = new SteppingHelixPropagator(&*bField,anyDirection);
  
  
  Propagator* theTmpPropagator = &*thePropagator->clone();
  theTmpPropagator->setPropagationDirection(alongMomentum); 
  
  
  std::vector<const DetLayer*>::const_iterator itl;
  findDetLayer = false;
  for (itl=trackCompatibleLayers.begin();itl!=trackCompatibleLayers.end();++itl) {
    TrajectoryStateOnSurface tsos = theTmpPropagator->propagate(*(sTSOS.freeState()),(**itl).surface());
    if (tsos.isValid()) 
    {
      sTSOS = tsos;
      findDetLayer = true; 
      
      break;
    }
  }
  std::pair<TrajectoryStateOnSurface,  const DetLayer* >  blabla;
  blabla.first = sTSOS;
  blabla.second = &**itl;
  delete theTmpPropagator;
  return blabla;
}

DEFINE_FWK_MODULE(TrackEfficiencyMonitor);