TrackExtenderWithMTD.cc

Go to the documentation of this file.

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h"
#include "RecoMTD/Records/interface/MTDRecoGeometryRecord.h"

#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"

#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"

#include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h"

#include "RecoMTD/DetLayers/interface/MTDTrayBarrelLayer.h"
#include "RecoMTD/DetLayers/interface/MTDDetTray.h"
#include "RecoMTD/DetLayers/interface/MTDRingForwardDoubleLayer.h"
#include "RecoMTD/DetLayers/interface/MTDDetRing.h"

#include "DataFormats/ForwardDetId/interface/BTLDetId.h"
#include "DataFormats/ForwardDetId/interface/ETLDetId.h"
#include "DataFormats/ForwardDetId/interface/MTDChannelIdentifier.h"
#include "Geometry/CommonTopologies/interface/PixelTopology.h"

#include "TrackingTools/PatternTools/interface/Trajectory.h"

#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"

#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"

#include "RecoMTD/TransientTrackingRecHit/interface/MTDTransientTrackingRecHitBuilder.h"
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"

#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"

#include "TrackingTools/PatternTools/interface/TSCBLBuilderWithPropagator.h"

#include "RecoTracker/TransientTrackingRecHit/interface/Traj2TrackHits.h"
#include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"

#include <sstream>

#include "Geometry/CommonTopologies/interface/Topology.h"

#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "CLHEP/Units/GlobalPhysicalConstants.h"
#include "CLHEP/Units/GlobalSystemOfUnits.h"

using namespace std;
using namespace edm;




template<class TrackCollection>
class TrackExtenderWithMTDT : public edm::stream::EDProducer<> {  
 public:
  typedef typename TrackCollection:: value_type TrackType;
  typedef edm::View<TrackType> InputCollection;
  
  TrackExtenderWithMTDT(const ParameterSet& pset); 

  template <class H, class T>
  void fillValueMap(edm::Event& iEvent, const H& handle, const std::vector<T>& vec, const std::string& name) const;
  
  void produce(edm::Event& ev, const edm::EventSetup& es) final;

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

  TransientTrackingRecHit::ConstRecHitContainer tryBTLLayers(const TrackType&, 
                                 const MTDTrackingDetSetVector&,
                                 const MTDDetLayerGeometry*, 
                                 const MagneticField* field,
                                 const Propagator* prop) const;

  TransientTrackingRecHit::ConstRecHitContainer tryETLLayers(const TrackType&,
                                 const MTDTrackingDetSetVector&,
                                 const MTDDetLayerGeometry*, 
                                 const MagneticField* field,
                                 const Propagator* prop) const;
  
  RefitDirection::GeometricalDirection
  checkRecHitsOrdering(TransientTrackingRecHit::ConstRecHitContainer const & recHits) const {
    
    if (!recHits.empty()){
      GlobalPoint first = gtg->idToDet(recHits.front()->geographicalId())->position();
      GlobalPoint last = gtg->idToDet(recHits.back()->geographicalId())->position();
      
      // maybe perp2?
      auto rFirst = first.mag2();
      auto rLast  = last.mag2();
      if(rFirst < rLast) return RefitDirection::insideOut;
      if(rFirst > rLast) return RefitDirection::outsideIn;
    }
    LogDebug("TrackExtenderWithMTD") << "Impossible to determine the rechits order" <<endl;
    return RefitDirection::undetermined;
  }

  reco::Track buildTrack(const reco::Track&, const Trajectory&, const Trajectory &, const reco::BeamSpot&, const MagneticField* field,const Propagator* prop, bool hasMTD, float& pathLength, float& tmtdOut, float& sigmatmtdOut) const;
  reco::TrackExtra buildTrackExtra(const Trajectory& trajectory) const;

  string dumpLayer(const DetLayer* layer) const;

 private:
  static constexpr char pathLengthName[] = "pathLength";
  static constexpr char tmtdName[] = "tmtd";
  static constexpr char sigmatmtdName[] = "sigmatmtd";  
  static constexpr char pOrigTrkName[] = "generalTrackp";
  static constexpr char betaOrigTrkName[] = "generalTrackBeta";
  static constexpr char t0OrigTrkName[] = "generalTrackt0";
  static constexpr char sigmat0OrigTrkName[] = "generalTracksigmat0";
  static constexpr char pathLengthOrigTrkName[] = "generalTrackPathLength";
  static constexpr char tmtdOrigTrkName[] = "generalTracktmtd";
  static constexpr char sigmatmtdOrigTrkName[] = "generalTracksigmatmtd";
   
  edm::EDGetTokenT<InputCollection> tracksToken_;
  edm::EDGetTokenT<MTDTrackingDetSetVector> hitsToken_;
  edm::EDGetTokenT<reco::BeamSpot> bsToken_;
  const bool updateTraj_, updateExtra_, updatePattern_;
  const std::string mtdRecHitBuilder_,propagator_, transientTrackBuilder_;
  std::unique_ptr<MeasurementEstimator> theEstimator;
  std::unique_ptr<TrackTransformer> theTransformer;
  edm::ESHandle<TransientTrackBuilder> builder;
  edm::ESHandle<TransientTrackingRecHitBuilder> hitbuilder;
  edm::ESHandle<GlobalTrackingGeometry> gtg;
  edm::ESHandle<Propagator> prop;
};


template<class TrackCollection>  
TrackExtenderWithMTDT<TrackCollection>::TrackExtenderWithMTDT(const ParameterSet& iConfig) :
  tracksToken_(consumes<InputCollection>(iConfig.getParameter<edm::InputTag>("tracksSrc"))),
  hitsToken_(consumes<MTDTrackingDetSetVector>(iConfig.getParameter<edm::InputTag>("hitsSrc"))),
  bsToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpotSrc"))),
  updateTraj_(iConfig.getParameter<bool>("updateTrackTrajectory")),
  updateExtra_(iConfig.getParameter<bool>("updateTrackExtra")),
  updatePattern_(iConfig.getParameter<bool>("updateTrackHitPattern")),
  mtdRecHitBuilder_(iConfig.getParameter<std::string>("MTDRecHitBuilder")),
  propagator_(iConfig.getParameter<std::string>("Propagator")),
  transientTrackBuilder_(iConfig.getParameter<std::string>("TransientTrackBuilder")) {
  constexpr float maxChi2=500.;
  constexpr float nSigma=10.;
  theEstimator = std::make_unique<Chi2MeasurementEstimator>(maxChi2,nSigma);
  
  theTransformer = std::make_unique<TrackTransformer>(iConfig.getParameterSet("TrackTransformer"));
  
  produces<edm::ValueMap<float> >(pathLengthName);
  produces<edm::ValueMap<float> >(tmtdName);
  produces<edm::ValueMap<float> >(sigmatmtdName); 
  produces<edm::ValueMap<float> >(pOrigTrkName);
  produces<edm::ValueMap<float> >(betaOrigTrkName);
  produces<edm::ValueMap<float> >(t0OrigTrkName);
  produces<edm::ValueMap<float> >(sigmat0OrigTrkName);
  produces<edm::ValueMap<float> >(pathLengthOrigTrkName);
  produces<edm::ValueMap<float> >(tmtdOrigTrkName);
  produces<edm::ValueMap<float> >(sigmatmtdOrigTrkName);
 
  produces<edm::OwnVector<TrackingRecHit>>();
  produces<reco::TrackExtraCollection>();
  produces<TrackCollection>();
}

template<class TrackCollection>
void TrackExtenderWithMTDT<TrackCollection>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc,transDesc;
  desc.add<edm::InputTag>("tracksSrc",edm::InputTag("generalTracks"));
  desc.add<edm::InputTag>("hitsSrc",edm::InputTag("mtdTrackingRecHits"));
  desc.add<edm::InputTag>("beamSpotSrc",edm::InputTag("offlineBeamSpot"));
  desc.add<bool>("updateTrackTrajectory",true);
  desc.add<bool>("updateTrackExtra",true);
  desc.add<bool>("updateTrackHitPattern",true);
  desc.add<std::string>("TransientTrackBuilder","TransientTrackBuilder");
  desc.add<std::string>("MTDRecHitBuilder","MTDRecHitBuilder");
  desc.add<std::string>("Propagator","PropagatorWithMaterialForMTD");  
  TrackTransformer::fillPSetDescription(transDesc,
                    false,
                    "KFFitterForRefitInsideOut",
                    "KFSmootherForRefitInsideOut",
                    "PropagatorWithMaterialForMTD",
                    "alongMomentum",
                    true,
                    "WithTrackAngle",
                    "MuonRecHitBuilder",
                    "MTDRecHitBuilder"
                    );
  desc.add<edm::ParameterSetDescription>("TrackTransformer",transDesc);
  descriptions.add("trackExtenderWithMTDBase", desc);
}

template<class TrackCollection>
template <class H, class T>
void TrackExtenderWithMTDT<TrackCollection>::fillValueMap(edm::Event& iEvent, const H& handle, const std::vector<T>& vec, const std::string& name) const {
  auto out = std::make_unique<edm::ValueMap<T>>();
  typename edm::ValueMap<T>::Filler filler(*out);
  filler.insert(handle, vec.begin(), vec.end());
  filler.fill();
  iEvent.put(std::move(out),name);
}

template<class TrackCollection>
void TrackExtenderWithMTDT<TrackCollection>::produce( edm::Event& ev,
                              const edm::EventSetup& es ) {  
  //this produces pieces of the track extra
  Traj2TrackHits t2t;

  theTransformer->setServices(es);

  TrackingRecHitRefProd hitsRefProd = ev.getRefBeforePut<TrackingRecHitCollection>();
  reco::TrackExtraRefProd extrasRefProd = ev.getRefBeforePut<reco::TrackExtraCollection>();
  
  es.get<GlobalTrackingGeometryRecord>().get(gtg);

  edm::ESHandle<MTDDetLayerGeometry> geo;
  es.get<MTDRecoGeometryRecord>().get(geo);

  edm::ESHandle<MagneticField> magfield;
  es.get<IdealMagneticFieldRecord>().get(magfield);  
    
  es.get<TransientTrackRecord>().get(transientTrackBuilder_, builder);
  es.get<TransientRecHitRecord>().get(mtdRecHitBuilder_,hitbuilder);

  edm::ESHandle<Propagator> prop;
  es.get<TrackingComponentsRecord>().get(propagator_,prop);

  edm::ESHandle<TrackerTopology> httopo;
  es.get<TrackerTopologyRcd>().get(httopo);
  const TrackerTopology& ttopo = *httopo;
    
  auto output  = std::make_unique<TrackCollection>();
  auto extras  = std::make_unique<reco::TrackExtraCollection>();
  auto outhits = std::make_unique<edm::OwnVector<TrackingRecHit>>();

  std::vector<float> pathLengthsRaw;
  std::vector<float> tmtdRaw;
  std::vector<float> sigmatmtdRaw;
  std::vector<float> pOrigTrkRaw;
  std::vector<float> betaOrigTrkRaw;
  std::vector<float> t0OrigTrkRaw;
  std::vector<float> sigmat0OrigTrkRaw;
  std::vector<float> pathLengthsOrigTrkRaw;
  std::vector<float> tmtdOrigTrkRaw;
  std::vector<float> sigmatmtdOrigTrkRaw;  
  
  edm::Handle<InputCollection> tracksH;  
  ev.getByToken(tracksToken_,tracksH);
  const auto& tracks = *tracksH;

  edm::Handle<MTDTrackingDetSetVector> hitsH;  
  ev.getByToken(hitsToken_,hitsH);
  const auto& hits = *hitsH;

  edm::Handle<reco::BeamSpot> bsH;  
  ev.getByToken(bsToken_,bsH);
  const auto& bs = *bsH;

  std::vector<unsigned> track_indices;
  unsigned itrack = 0;
  for( const auto& track : tracks ) { 
    
    reco::TransientTrack ttrack(track,magfield.product(),gtg);
    const auto& trajs = theTransformer->transform(track);
    auto thits = theTransformer->getTransientRecHits(ttrack);

    TransientTrackingRecHit::ConstRecHitContainer mtdthits;
    const auto& btlhits = tryBTLLayers(track,hits,geo.product(),magfield.product(),prop.product());
    mtdthits.insert(mtdthits.end(),btlhits.begin(),btlhits.end());
    
    // in the future this should include an intermediate refit before propagating to the ETL
    // for now it is ok
    const auto& etlhits = tryETLLayers(track,hits,geo.product(),magfield.product(),prop.product());
    mtdthits.insert(mtdthits.end(),etlhits.begin(),etlhits.end());
    
    auto ordering = checkRecHitsOrdering(thits);
    if( ordering == RefitDirection::insideOut) {
      thits.insert(thits.end(),mtdthits.begin(),mtdthits.end());
    } else {
      std::reverse(mtdthits.begin(),mtdthits.end());
      mtdthits.insert(mtdthits.end(),thits.begin(),thits.end());
      thits.swap(mtdthits);
    }
    const auto& trajwithmtd = theTransformer->transform(ttrack,thits);
    float pMap = 0.f, betaMap = 0.f, t0Map = 0.f, sigmat0Map = -1.f, pathLengthMap = -1.f, tmtdMap = 0.f, sigmatmtdMap = -1.f;
    
    for( const auto& trj : trajwithmtd ) {      
      const auto& thetrj = (updateTraj_ ? trj : trajs.front());
      float pathLength = 0.f, tmtd = 0.f, sigmatmtd = -1.f;
      reco::Track result = buildTrack(track, thetrj, trj, bs, magfield.product(), 
                      prop.product(), !mtdthits.empty(),pathLength,tmtd,sigmatmtd);
      if( result.ndof() >= 0 ) {
        reco::TrackExtra::TrajParams trajParams;
        reco::TrackExtra::Chi2sFive chi2s; 
        size_t hitsstart = outhits->size();
        if( updatePattern_ ) {
          t2t(trj,*outhits,trajParams,chi2s); // this fills the output hit collection
        } else {
          t2t(thetrj,*outhits,trajParams,chi2s);
        }
        size_t hitsend = outhits->size();
        extras->push_back(buildTrackExtra(trj)); // always push back the fully built extra, update by setting in track
        extras->back().setHits(hitsRefProd,hitsstart,hitsend-hitsstart);    
        extras->back().setTrajParams(trajParams,chi2s);        
    //create the track
    output->push_back(result);
    pathLengthsRaw.push_back(pathLength);
    tmtdRaw.push_back(tmtd);
    sigmatmtdRaw.push_back(sigmatmtd);
    pathLengthMap = pathLength;
    tmtdMap = tmtd;
    sigmatmtdMap = sigmatmtd;
        auto& backtrack = output->back();
    pMap = backtrack.p();
    betaMap = backtrack.beta();
    t0Map = backtrack.t0();
    sigmat0Map = std::copysign(std::sqrt(std::abs(backtrack.covt0t0())),backtrack.covt0t0());
        reco::TrackExtraRef extraRef(extrasRefProd,extras->size()-1);
        backtrack.setExtra( (updateExtra_ ? extraRef : track.extra()) );
        for(unsigned ihit = hitsstart; ihit < hitsend; ++ihit) {
          backtrack.appendHitPattern((*outhits)[ihit],ttopo);
        }
      }      
    }
    pOrigTrkRaw.push_back(pMap);
    betaOrigTrkRaw.push_back(betaMap);
    t0OrigTrkRaw.push_back(t0Map);
    sigmat0OrigTrkRaw.push_back(sigmat0Map);
    pathLengthsOrigTrkRaw.push_back(pathLengthMap);
    tmtdOrigTrkRaw.push_back(tmtdMap);
    sigmatmtdOrigTrkRaw.push_back(sigmatmtdMap);
    ++itrack;
  }

  auto outTrksHandle = ev.put(std::move(output));
  ev.put(std::move(extras));
  ev.put(std::move(outhits));

  fillValueMap(ev, outTrksHandle, pathLengthsRaw, pathLengthName);
  fillValueMap(ev, outTrksHandle, tmtdRaw, tmtdName);
  fillValueMap(ev, outTrksHandle, sigmatmtdRaw, sigmatmtdName);
  fillValueMap(ev, tracksH, pOrigTrkRaw, pOrigTrkName);
  fillValueMap(ev, tracksH, betaOrigTrkRaw, betaOrigTrkName);
  fillValueMap(ev, tracksH, t0OrigTrkRaw, t0OrigTrkName);
  fillValueMap(ev, tracksH, sigmat0OrigTrkRaw, sigmat0OrigTrkName);
  fillValueMap(ev, tracksH, pathLengthsOrigTrkRaw, pathLengthOrigTrkName);
  fillValueMap(ev, tracksH, tmtdOrigTrkRaw, tmtdOrigTrkName);
  fillValueMap(ev, tracksH, sigmatmtdOrigTrkRaw, sigmatmtdOrigTrkName);
}

namespace {
  bool cmp_for_detset(const unsigned one, const unsigned two) { return one < two; };
  
  void find_hits_in_dets(const MTDTrackingDetSetVector& hits, const DetLayer* layer,
             const TrajectoryStateOnSurface& tsos, const Propagator* prop,
             const MeasurementEstimator& theEstimator,             
             const TransientTrackingRecHitBuilder& hitbuilder,
             TransientTrackingRecHit::ConstRecHitContainer& output) {
    pair<bool, TrajectoryStateOnSurface> comp = layer->compatible(tsos,*prop,theEstimator);
    if( comp.first ) {    
      vector<DetLayer::DetWithState> compDets = layer->compatibleDets(tsos,*prop,theEstimator);
      if (!compDets.empty()) {
    MTDTrackingRecHit* best = nullptr;
    double best_chi2 = std::numeric_limits<double>::max();
    for( const auto& detWithState : compDets ) {    
      auto range = hits.equal_range(detWithState.first->geographicalId(),cmp_for_detset);     
      for( auto detitr = range.first; detitr != range.second; ++detitr ) {
        for( auto itr = detitr->begin(); itr != detitr->end(); ++itr ) {
          auto est =  theEstimator.estimate(detWithState.second,*itr);
          if( est.first && est.second < best_chi2 ) { // just take the best chi2
        best = &(*itr);
        best_chi2 = est.second;
          }
          
        }
      }
    }
    if( best ) 
      output.push_back(hitbuilder.build(best));
      }
    }
  }
}

template<class TrackCollection>
TransientTrackingRecHit::ConstRecHitContainer
TrackExtenderWithMTDT<TrackCollection>::tryBTLLayers(const TrackType& track,
                             const MTDTrackingDetSetVector& hits,
                             const MTDDetLayerGeometry* geo,
                             const MagneticField* field,
                             const Propagator* prop) const {
  
  TransientTrackingRecHit::ConstRecHitContainer output;
  const vector<const DetLayer*>& layers = geo->allBTLLayers();
  auto tTrack = builder->build(track);

  for (const DetLayer* ilay : layers) {
    // get the outermost trajectory point on the track    
    TrajectoryStateOnSurface tsos = tTrack.outermostMeasurementState();
    find_hits_in_dets(hits,ilay,tsos,prop,*theEstimator,*hitbuilder,output);
  }
  return output;
}

template<class TrackCollection>
TransientTrackingRecHit::ConstRecHitContainer
TrackExtenderWithMTDT<TrackCollection>::tryETLLayers(const TrackType& track, 
                             const MTDTrackingDetSetVector& hits,
                             const MTDDetLayerGeometry* geo,
                             const MagneticField* field,
                             const Propagator* prop) const {

  TransientTrackingRecHit::ConstRecHitContainer output;
  const vector<const DetLayer*>& layers = geo->allETLLayers();

  auto tTrack = builder->build(track);
  
  for (const DetLayer* ilay : layers) {
    const BoundDisk& disk = static_cast<const MTDRingForwardDoubleLayer*>(ilay)->specificSurface();
    const double diskZ = disk.position().z();

    // get the outermost trajectory point on the track    
    TrajectoryStateOnSurface tsos = tTrack.outermostMeasurementState();
    if( tsos.globalPosition().z() * diskZ < 0 ) continue; // only propagate to the disk that's on the same side
    find_hits_in_dets(hits,ilay,tsos,prop,*theEstimator,*hitbuilder,output);    
  }
  return output;
}


//below is unfortunately ripped from other places but 
//since track producer doesn't know about MTD we have to do this
template<class TrackCollection>
reco::Track TrackExtenderWithMTDT<TrackCollection>::buildTrack(const reco::Track& orig,
                                   const Trajectory& traj,
                                   const Trajectory& trajWithMtd,
                                   const reco::BeamSpot& bs,
                                   const MagneticField* field,
                                   const Propagator* thePropagator,
                                   bool hasMTD,
                                   float& pathLengthOut,
                                   float& tmtdOut,
                                   float& sigmatmtdOut) const {  

  // get the state closest to the beamline
  TrajectoryStateOnSurface stateForProjectionToBeamLineOnSurface = 
    traj.closestMeasurement(GlobalPoint(bs.x0(),bs.y0(),bs.z0())).updatedState();
  
  if UNLIKELY(!stateForProjectionToBeamLineOnSurface.isValid()) {
    edm::LogError("CannotPropagateToBeamLine")<<"the state on the closest measurement isnot valid. skipping track.";
    return reco::Track();
  }

  constexpr double m_pi = 0.13957018;
  constexpr double m_pi_inv2 = 1.0/m_pi/m_pi;
  constexpr double m_p = 0.9382720813;
  constexpr double m_p_inv2 = 1.0/m_p/m_p;
  constexpr double c_cm_ns = CLHEP::c_light*CLHEP::ns/CLHEP::cm; //[cm/ns]
  constexpr double c_inv = 1.0/c_cm_ns;
  
  const FreeTrajectoryState & stateForProjectionToBeamLine=*stateForProjectionToBeamLineOnSurface.freeState();
  
  TSCBLBuilderWithPropagator tscblBuilder(*thePropagator);
  TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine,bs);
  
  if UNLIKELY(!tscbl.isValid()) {
    return reco::Track();
  }

  GlobalPoint v = tscbl.trackStateAtPCA().position();
  math::XYZPoint  pos( v.x(), v.y(), v.z() );
  GlobalVector p = tscbl.trackStateAtPCA().momentum();
  math::XYZVector mom( p.x(), p.y(), p.z() );

  int ndof = traj.ndof();
  
  double t0 = 0.;
  double covt0t0 = -1.;
  pathLengthOut = -1.f; // if there is no MTD flag the pathlength with -1
  tmtdOut = 0.f;
  sigmatmtdOut = -1.f;
  double betaOut = 0.;
  double covbetabeta = -1.;

  //compute path length for time backpropagation, using first MTD hit for the momentum
  if (hasMTD) {
    
    bool validpropagation = true;
    double pathlength = 0.;
    double pathlength1 = 0.;
    double pathlength2 = 0.;
    for (auto it=trajWithMtd.measurements().begin(); it!=trajWithMtd.measurements().end()-1; ++it) {
       const auto &propresult = thePropagator->propagateWithPath(it->updatedState(), (it+1)->updatedState().surface());
       double layerpathlength = std::abs(propresult.second);
       if (layerpathlength==0.) {
         validpropagation = false;
       }
       pathlength1 += layerpathlength;
    }
     
    double thit = 0.;
    double thiterror = -1.;
    bool validmtd = false;
    if (trajWithMtd.direction() == alongMomentum) {
      for (auto it=trajWithMtd.measurements().begin(); it!=trajWithMtd.measurements().end(); ++it) {
        bool ismtd = it->recHit()->geographicalId().det() == DetId::Forward && ForwardSubdetector(it->recHit()->geographicalId().subdetId()) == FastTime;
        if (ismtd) {
          const auto &propresult2 = thePropagator->propagateWithPath(tscbl.trackStateAtPCA(), trajWithMtd.firstMeasurement().updatedState().surface());
      pathlength2 = propresult2.second;
          if (pathlength2 == 0.) {
            validpropagation = false;
          }
          pathlength = pathlength1 + pathlength2; 
          const MTDTrackingRecHit *mtdhit = static_cast<const MTDTrackingRecHit*>(it->recHit()->hit());
          thit = mtdhit->time();
          thiterror = mtdhit->timeError();
          validmtd = true;
          break;          
        }
      }
    }
    else {
      for (auto it=trajWithMtd.measurements().rbegin(); it!=trajWithMtd.measurements().rend(); ++it) {
        bool ismtd = it->recHit()->geographicalId().det() == DetId::Forward && ForwardSubdetector(it->recHit()->geographicalId().subdetId()) == FastTime;
        if (ismtd) {
          const auto &propresult2 = thePropagator->propagateWithPath(tscbl.trackStateAtPCA(), trajWithMtd.lastMeasurement().updatedState().surface());
          pathlength2 = propresult2.second;
          if (pathlength2 == 0.) {
            validpropagation = false;
          }
          pathlength = pathlength1 + pathlength2;
          const MTDTrackingRecHit *mtdhit = static_cast<const MTDTrackingRecHit*>(it->recHit()->hit());
          thit = mtdhit->time();
          thiterror = mtdhit->timeError();
          validmtd = true;
          break;          
        }
      }
    }
        
    if (validmtd && validpropagation) {
      double magp2 = p.mag2();

      double gammasq_pi = 1. + magp2*m_pi_inv2;
      double beta_pi = std::sqrt(1.-1./gammasq_pi);
      double dt_pi = pathlength/beta_pi*c_inv;
      
      double gammasq_p = 1. + magp2*m_p_inv2;
      double beta_p = std::sqrt(1.-1./gammasq_p);
      double dt_p = pathlength/beta_p*c_inv;
      
      double dterror = dt_p - dt_pi;
      double betaerror = beta_p - beta_pi;
      
      pathLengthOut = pathlength; // set path length if we've got a timing hit
      tmtdOut = thit;
      sigmatmtdOut = thiterror;
      t0 = thit - dt_pi;
      covt0t0 = thiterror*thiterror + dterror*dterror;
      betaOut = beta_pi;
      covbetabeta  = betaerror*betaerror;
    }
    
  }
  
  return reco::Track(traj.chiSquared(),
             int(ndof),
             pos, mom, tscbl.trackStateAtPCA().charge(), 
             tscbl.trackStateAtPCA().curvilinearError(),
             orig.algo(),reco::TrackBase::undefQuality,t0,betaOut,covt0t0,covbetabeta);
}

template<class TrackCollection>
reco::TrackExtra TrackExtenderWithMTDT<TrackCollection>::buildTrackExtra(const Trajectory& trajectory) const {

  static const string metname = "TrackExtenderWithMTD";

  const Trajectory::RecHitContainer transRecHits = trajectory.recHits();
  
  // put the collection of TrackingRecHit in the event
  
  // sets the outermost and innermost TSOSs
  // ToDo: validation for track states with MTD
  TrajectoryStateOnSurface outerTSOS;
  TrajectoryStateOnSurface innerTSOS;
  unsigned int innerId=0, outerId=0;
  TrajectoryMeasurement::ConstRecHitPointer outerRecHit;
  DetId outerDetId;

  if (trajectory.direction() == alongMomentum) {
    LogTrace(metname)<<"alongMomentum";
    outerTSOS = trajectory.lastMeasurement().updatedState();
    innerTSOS = trajectory.firstMeasurement().updatedState();
    outerId = trajectory.lastMeasurement().recHit()->geographicalId().rawId();
    innerId = trajectory.firstMeasurement().recHit()->geographicalId().rawId();
    outerRecHit =  trajectory.lastMeasurement().recHit();
    outerDetId =   trajectory.lastMeasurement().recHit()->geographicalId();
  } 
  else if (trajectory.direction() == oppositeToMomentum) {
    LogTrace(metname)<<"oppositeToMomentum";
    outerTSOS = trajectory.firstMeasurement().updatedState();
    innerTSOS = trajectory.lastMeasurement().updatedState();
    outerId = trajectory.firstMeasurement().recHit()->geographicalId().rawId();
    innerId = trajectory.lastMeasurement().recHit()->geographicalId().rawId();
    outerRecHit =  trajectory.firstMeasurement().recHit();
    outerDetId =   trajectory.firstMeasurement().recHit()->geographicalId();
  }
  else LogError(metname)<<"Wrong propagation direction!";
  
  const GeomDet *outerDet = gtg->idToDet(outerDetId);
  GlobalPoint outerTSOSPos = outerTSOS.globalParameters().position();
  bool inside = outerDet->surface().bounds().inside(outerDet->toLocal(outerTSOSPos));

  
  GlobalPoint hitPos = (outerRecHit->isValid()) ? outerRecHit->globalPosition() :  outerTSOS.globalParameters().position() ;
  
  if(!inside) {
    LogTrace(metname)
      << "The Global Muon outerMostMeasurementState is not compatible with the recHit detector!"
      << " Setting outerMost postition to recHit position if recHit isValid: " << outerRecHit->isValid();
    LogTrace(metname)<<"From " << outerTSOSPos << " to " <<  hitPos;
  }
  
  
  //build the TrackExtra
  GlobalPoint v = (inside) ? outerTSOSPos : hitPos ;
  GlobalVector p = outerTSOS.globalParameters().momentum();
  math::XYZPoint  outpos( v.x(), v.y(), v.z() );   
  math::XYZVector outmom( p.x(), p.y(), p.z() );
  
  v = innerTSOS.globalParameters().position();
  p = innerTSOS.globalParameters().momentum();
  math::XYZPoint  inpos( v.x(), v.y(), v.z() );   
  math::XYZVector inmom( p.x(), p.y(), p.z() );

  reco::TrackExtra trackExtra(outpos, outmom, true, inpos, inmom, true,
                              outerTSOS.curvilinearError(), outerId,
                              innerTSOS.curvilinearError(), innerId,
                  trajectory.direction(),trajectory.seedRef());
  
  return trackExtra;
 
}

template<class TrackCollection>
string TrackExtenderWithMTDT<TrackCollection>::dumpLayer(const DetLayer* layer) const {
  stringstream output;
  
  const BoundSurface* sur=nullptr;
  const BoundCylinder* bc=nullptr;
  const BoundDisk* bd=nullptr;

  sur = &(layer->surface());
  if ( (bc = dynamic_cast<const BoundCylinder*>(sur)) ) {
    output << "  Cylinder of radius: " << bc->radius() << endl;
  }
  else if ( (bd = dynamic_cast<const BoundDisk*>(sur)) ) {
    output << "  Disk at: " <<  bd->position().z() << endl;
  }
  return output.str();
}

//define this as a plug-in
#include <FWCore/Framework/interface/MakerMacros.h>
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
typedef TrackExtenderWithMTDT<reco::TrackCollection> TrackExtenderWithMTD;

DEFINE_FWK_MODULE(TrackExtenderWithMTD);