MuonTrackLoader.cc

Go to the documentation of this file.

#include "RecoMuon/TrackingTools/interface/MuonTrackLoader.h"

#include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
#include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
#include "RecoMuon/TrackingTools/interface/MuonUpdatorAtVertex.h"

#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/GeomPropagators/interface/TrackerBounds.h"
#include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
#include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"

#include "Geometry/CommonDetUnit/interface/GeomDet.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackExtra.h"
#include "DataFormats/TrackReco/interface/TrackToTrackMap.h"

#include "DataFormats/MuonReco/interface/MuonTrackLinks.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"

#include "TrackingTools/Records/interface/TransientRecHitRecord.h"
#include "RecoTracker/TransientTrackingRecHit/interface/Traj2TrackHits.h"


using namespace edm;
using namespace std;
using namespace reco;

std::vector<const TrackingRecHit*> MuonTrackLoader::unpackHit(const TrackingRecHit &hit)
{
    // get rec hit det id and rec hit type
    DetId id = hit.geographicalId();
    uint16_t detid = id.det();

    std::vector<const TrackingRecHit*> hits;

    if (detid == DetId::Tracker) {
        hits.push_back(&hit);
    } else if (detid == DetId::Muon) {
        uint16_t subdet = id.subdetId();
        if (subdet == (uint16_t) MuonSubdetId::DT) {
            if (hit.dimension() == 1) { // DT rechit (granularity 2)
                hits.push_back(&hit);
            } else if (hit.dimension() > 1) { // 4D segment (granularity 0).
                // Both 2 and 4 dim cases. MB4s have 2D, but formatted in 4D segment
                // 4D --> 2D
                std::vector<const TrackingRecHit*> seg2D = hit.recHits();
                // load 1D hits (2D --> 1D)
                for (std::vector<const TrackingRecHit*>::const_iterator it = seg2D.begin(); it != seg2D.end(); ++it) {
                    std::vector<const TrackingRecHit*> hits1D = (*it)->recHits();
                    copy(hits1D.begin(), hits1D.end(), back_inserter(hits));
                }
            }
        } else if (subdet == (uint16_t) MuonSubdetId::CSC) {
            if (hit.dimension() == 2) { // CSC rechit (granularity 2)
                hits.push_back(&hit);
            } else if (hit.dimension() == 4) { // 4D segment (granularity 0)
                // load 2D hits (4D --> 1D)
                hits = hit.recHits();
            }
        } else if (subdet == (uint16_t) MuonSubdetId::RPC) {
            hits.push_back(&hit);
        }
        else if (subdet == (uint16_t) MuonSubdetId::GEM) {
        if (hit.dimension() == 2) { // GEM rechit
                hits.push_back(&hit);
            } else if (hit.dimension() == 4) { // GEM segment
                hits = hit.recHits();
            }
        }
    else if (subdet == (uint16_t) MuonSubdetId::ME0) { //segment
            hits = hit.recHits();
        }
    }
    return hits;
}

// constructor (obsolete, should use eventSetUp not service..)
MuonTrackLoader::MuonTrackLoader(ParameterSet &parameterSet, edm::ConsumesCollector& iC, const MuonServiceProxy *service): 
  theService(service){


  // option to do or not the smoothing step.
  // the trajectories which are passed to the track loader are supposed to be non-smoothed
  theSmoothingStep = parameterSet.getParameter<bool>("DoSmoothing");
  if(theSmoothingStep)
    theSmootherName = parameterSet.getParameter<string>("Smoother");  

  theTrackerRecHitBuilderName = parameterSet.getParameter<std::string>("TTRHBuilder");
  
  // update at vertex
  theUpdatingAtVtx = parameterSet.getParameter<bool>("VertexConstraint");

  // beam spot input tag
  theBeamSpotInputTag = parameterSet.getParameter<edm::InputTag>("beamSpot");
  theBeamSpotToken = iC.consumes<reco::BeamSpot>(theBeamSpotInputTag);


  // Flag to put the trajectory into the event
  theTrajectoryFlag = parameterSet.getUntrackedParameter<bool>("PutTrajectoryIntoEvent",true);

  theL2SeededTkLabel = parameterSet.getUntrackedParameter<string>("MuonSeededTracksInstance",string());
  
  ParameterSet updatorPar = parameterSet.getParameter<ParameterSet>("MuonUpdatorAtVertexParameters");
  theUpdatorAtVtx = new MuonUpdatorAtVertex(updatorPar,service);

  thePutTkTrackFlag = parameterSet.getUntrackedParameter<bool>("PutTkTrackIntoEvent",false);
  theSmoothTkTrackFlag = parameterSet.getUntrackedParameter<bool>("SmoothTkTrack",false);
  theAllowNoVtxFlag = parameterSet.getUntrackedParameter<bool>("AllowNoVertex",false);
}

MuonTrackLoader::~MuonTrackLoader(){
  if(theUpdatorAtVtx) delete theUpdatorAtVtx;
}

OrphanHandle<reco::TrackCollection> 
MuonTrackLoader::loadTracks(const TrajectoryContainer& trajectories,
                Event& event, const TrackerTopology& ttopo, const string& instance, bool reallyDoSmoothing) {
  std::vector<bool> dummyVecBool;
  return loadTracks(trajectories, event, dummyVecBool, ttopo, instance, reallyDoSmoothing);
}
  
OrphanHandle<reco::TrackCollection> 
MuonTrackLoader::loadTracks(const TrajectoryContainer& trajectories,
                Event& event,  std::vector<bool>& tkBoolVec, 
                            const TrackerTopology& ttopo,
                const string& instance, bool reallyDoSmoothing) {
  
  const bool doSmoothing = theSmoothingStep && reallyDoSmoothing;
  
  const string metname = "Muon|RecoMuon|MuonTrackLoader";
  
  // the track collectios; they will be loaded in the event  
  auto trackCollection = std::make_unique<reco::TrackCollection>();
  // ... and its reference into the event
  reco::TrackRefProd trackCollectionRefProd = event.getRefBeforePut<reco::TrackCollection>(instance);
  
  // track collection for the tracks updated at vertex
  auto updatedAtVtxTrackCollection = std::make_unique<reco::TrackCollection>();
  // ... and its (eventually) reference into the event
  reco::TrackRefProd trackUpdatedCollectionRefProd;
  if(theUpdatingAtVtx)  trackUpdatedCollectionRefProd = event.getRefBeforePut<reco::TrackCollection>(instance+"UpdatedAtVtx");
  
  // Association map between updated and non updated at vtx tracks
  auto trackToTrackmap = std::make_unique<reco::TrackToTrackMap>(trackCollectionRefProd, trackUpdatedCollectionRefProd);
  
  // the track extra collection, it will be loaded in the event  
  auto trackExtraCollection = std::make_unique<reco::TrackExtraCollection>();
  // ... and its reference into the event
  reco::TrackExtraRefProd trackExtraCollectionRefProd = event.getRefBeforePut<reco::TrackExtraCollection>(instance);
  
  // the rechit collection, it will be loaded in the event  
  auto recHitCollection = std::make_unique<TrackingRecHitCollection>();
  // ... and its reference into the event
  TrackingRecHitRefProd recHitCollectionRefProd = event.getRefBeforePut<TrackingRecHitCollection>(instance);
  
  // Collection of Trajectory
  auto trajectoryCollection = std::make_unique<vector<Trajectory>>();
  
  // don't waste any time...
  if ( trajectories.empty() ) { 
    event.put(std::move(recHitCollection),instance);
    event.put(std::move(trackExtraCollection),instance);
    if(theTrajectoryFlag) {
      event.put(std::move(trajectoryCollection),instance);

      // Association map between track and trajectory
      auto trajTrackMap = std::make_unique<TrajTrackAssociationCollection>();
      event.put(std::move(trajTrackMap), instance );
    }
    if(theUpdatingAtVtx){
      event.put(std::move(trackToTrackmap));
      event.put(std::move(updatedAtVtxTrackCollection),instance+"UpdatedAtVtx");
    }
    return event.put(std::move(trackCollection),instance);
  }
  
  edm::Handle<reco::BeamSpot> beamSpot;
  event.getByToken(theBeamSpotToken, beamSpot);

  LogTrace(metname) << "Create the collection of Tracks";
  
  reco::TrackRef::key_type trackIndex = 0;
  reco::TrackRef::key_type trackUpdatedIndex = 0;
  
  reco::TrackExtraRef::key_type trackExtraIndex = 0;
  
  edm::Ref<reco::TrackCollection>::key_type iTkRef = 0;
  edm::Ref< std::vector<Trajectory> >::key_type iTjRef = 0;
  std::map<unsigned int, unsigned int> tjTkMap;
  
  if(doSmoothing) {
    edm::ESHandle<TrajectorySmoother> aSmoother;
    theService->eventSetup().get<TrajectoryFitter::Record>().get(theSmootherName,aSmoother);
    theSmoother.reset(aSmoother->clone());
    edm::ESHandle<TransientTrackingRecHitBuilder> theTrackerRecHitBuilder;
    theService->eventSetup().get<TransientRecHitRecord>().get(theTrackerRecHitBuilderName,theTrackerRecHitBuilder);
    theTrackerRecHitBuilder.product();
    hitCloner = static_cast<TkTransientTrackingRecHitBuilder const *>(theTrackerRecHitBuilder.product())->cloner();
    theSmoother->setHitCloner(&hitCloner);
  }  


  unsigned int tjCnt = 0;
  for(TrajectoryContainer::const_iterator rawTrajectory = trajectories.begin();
      rawTrajectory != trajectories.end(); ++rawTrajectory, ++tjCnt){
    
    Trajectory &trajectory = **rawTrajectory;
    
    if(doSmoothing){
      vector<Trajectory> trajectoriesSM = theSmoother->trajectories(**rawTrajectory);
      
      if(!trajectoriesSM.empty()) {
    const edm::RefToBase<TrajectorySeed> tmpSeedRef = (**rawTrajectory).seedRef();
    trajectory = trajectoriesSM.front();
        trajectory.setSeedRef(tmpSeedRef);
    LogDebug(metname) << "theSeedRef.isNonnull " << trajectory.seedRef().isNonnull();
      } else
    LogInfo(metname)<<"The trajectory has not been smoothed!"<<endl;
    }
    
    if(theTrajectoryFlag) {
      trajectoryCollection->push_back(trajectory);
      iTjRef++;
    }    
    

    // build the "bare" track from the trajectory.
    // This track has the parameters defined at PCA (no update)
    pair<bool,reco::Track> resultOfTrackExtrapAtPCA = buildTrackAtPCA(trajectory, *beamSpot);

    // Check if the extrapolation went well
    if(!resultOfTrackExtrapAtPCA.first) {
      delete *rawTrajectory;
      continue;
    }


    // take the "bare" track at PCA
    reco::Track &track = resultOfTrackExtrapAtPCA.second;
    
    // build the "bare" track extra from the trajectory
    reco::TrackExtra trackExtra = buildTrackExtra( trajectory );
    
    // get the TrackExtraRef (persitent reference of the track extra)
    reco::TrackExtraRef trackExtraRef(trackExtraCollectionRefProd, trackExtraIndex++ );
    
    // set the persistent track-extra reference to the Track
    track.setExtra(trackExtraRef);
    
    // build the updated-at-vertex track, starting from the previous track
    pair<bool,reco::Track> updateResult(false,reco::Track());
    
    if(theUpdatingAtVtx){
      // build the "bare" track UPDATED at vtx
      updateResult = buildTrackUpdatedAtPCA(track, *beamSpot);

      if(!updateResult.first) ++trackIndex;
      else{
    
    // set the persistent track-extra reference to the Track
    updateResult.second.setExtra(trackExtraRef);
    
    // Fill the map
    trackToTrackmap->insert(reco::TrackRef(trackCollectionRefProd,trackIndex++),
                reco::TrackRef(trackUpdatedCollectionRefProd,trackUpdatedIndex++));
      }
    }
    

    // get the transient rechit and co from the trajectory
    reco::TrackExtra::TrajParams trajParams;
    reco::TrackExtra::Chi2sFive chi2s;
    Traj2TrackHits t2t;
    auto ih = recHitCollection->size();
    t2t(trajectory,*recHitCollection,trajParams,chi2s);
    auto ie = recHitCollection->size();
    // set the TrackingRecHitRef (persitent reference of the tracking rec hits)
    trackExtra.setHits(recHitCollectionRefProd,ih,ie-ih);
    trackExtra.setTrajParams(std::move(trajParams),std::move(chi2s));
    assert(trackExtra.trajParams().size()==trackExtra.recHitsSize());


    // Fill the hit pattern
    for (;ih<ie; ++ih) {
      auto const & hit = (*recHitCollection)[ih];
      auto hits = MuonTrackLoader::unpackHit(hit);
      for (auto hh : hits) {
          if UNLIKELY(!track.appendHitPattern(*hh, ttopo)) break;
      }
 
      if(theUpdatingAtVtx && updateResult.first){
        for (auto hh : hits) {
              if UNLIKELY(!updateResult.second.appendHitPattern(*hh, ttopo)) break;
        }
      }
    }

    // fill the TrackExtraCollection
    trackExtraCollection->push_back(trackExtra);
    
    // fill the TrackCollection
    trackCollection->push_back(track);
    iTkRef++;
    LogTrace(metname) << "Debug Track being loaded pt "<< track.pt();
    // fill the TrackCollection updated at vtx
    if(theUpdatingAtVtx && updateResult.first) 
      updatedAtVtxTrackCollection->push_back(updateResult.second);
    
    // We don't need the original trajectory anymore.
    // It has been copied by value in the trajectoryCollection, if 
    // it is required to put it into the event.
    delete *rawTrajectory;

    if(tkBoolVec.size()>tjCnt) tkBoolVec[tjCnt] = true;
    if(theTrajectoryFlag) tjTkMap[iTjRef-1] = iTkRef-1;
  }
  

  
  // Put the Collections in the event
  LogTrace(metname) << "put the Collections in the event";
  event.put(std::move(recHitCollection),instance);
  event.put(std::move(trackExtraCollection),instance);

  OrphanHandle<reco::TrackCollection> returnTrackHandle;
  OrphanHandle<reco::TrackCollection> nonUpdatedHandle;
  if(theUpdatingAtVtx){
    nonUpdatedHandle = event.put(std::move(trackCollection),instance);
    event.put(std::move(trackToTrackmap));
    returnTrackHandle = event.put(std::move(updatedAtVtxTrackCollection),instance+"UpdatedAtVtx");
  }
  else {
    returnTrackHandle = event.put(std::move(trackCollection),instance);
    nonUpdatedHandle = returnTrackHandle;
  }

  if ( theTrajectoryFlag ) {
    OrphanHandle<std::vector<Trajectory> > rTrajs = event.put(std::move(trajectoryCollection),instance);

    // Association map between track and trajectory
    auto trajTrackMap = std::make_unique<TrajTrackAssociationCollection>(rTrajs, nonUpdatedHandle);
  
    // Now Create traj<->tracks association map
    for ( std::map<unsigned int, unsigned int>::iterator i = tjTkMap.begin(); 
          i != tjTkMap.end(); i++ ) {
      trajTrackMap->insert( edm::Ref<std::vector<Trajectory> >( rTrajs, (*i).first ),
                            edm::Ref<reco::TrackCollection>( nonUpdatedHandle, (*i).second ) );
    }
    event.put(std::move(trajTrackMap), instance );
  }
  
  return returnTrackHandle;
}

OrphanHandle<reco::MuonTrackLinksCollection> 
MuonTrackLoader::loadTracks(const CandidateContainer& muonCands,
                Event& event,
                            const TrackerTopology& ttopo) {

  const string metname = "Muon|RecoMuon|MuonTrackLoader";
  
  // the muon collection, it will be loaded in the event
  auto trackLinksCollection = std::make_unique<reco::MuonTrackLinksCollection>();
  
  // don't waste any time...
  if ( muonCands.empty() ) {
    auto trackExtraCollection = std::make_unique<reco::TrackExtraCollection>();
    auto recHitCollection = std::make_unique<TrackingRecHitCollection>();
    auto trackCollection = std::make_unique<reco::TrackCollection>();

    event.put(std::move(recHitCollection));
    event.put(std::move(trackExtraCollection));
    event.put(std::move(trackCollection));

    //need to also put the tracker tracks collection if requested
    if(thePutTkTrackFlag){
      //will take care of putting nothing in the event but the empty collection
      TrajectoryContainer trackerTrajs;
      loadTracks(trackerTrajs, event, ttopo, theL2SeededTkLabel, theSmoothTkTrackFlag);
    } 

    return event.put(std::move(trackLinksCollection));
  }
  
  // get combined Trajectories
  TrajectoryContainer combinedTrajs;
  TrajectoryContainer trackerTrajs;
  for (CandidateContainer::const_iterator it = muonCands.begin(); it != muonCands.end(); ++it) {
    LogDebug(metname) << "Loader glbSeedRef " << (*it)->trajectory()->seedRef().isNonnull();
    if ((*it)->trackerTrajectory() )  LogDebug(metname) << " " << "tkSeedRef " << (*it)->trackerTrajectory()->seedRef().isNonnull();

    combinedTrajs.push_back((*it)->trajectory());
    if ( thePutTkTrackFlag ) trackerTrajs.push_back((*it)->trackerTrajectory());

    else {
      if ((*it)->trackerTrajectory()) delete ((*it)->trackerTrajectory());
    }
  
    // // Create the links between sta and tracker tracks
    // reco::MuonTrackLinks links;
    // links.setStandAloneTrack((*it)->muonTrack());
    // links.setTrackerTrack((*it)->trackerTrack());
    // trackLinksCollection->push_back(links);
    // delete *it;
  }
  
  // create the TrackCollection of combined Trajectories
  // FIXME: could this be done one track at a time in the previous loop?
  LogTrace(metname) << "Build combinedTracks";
  std::vector<bool> combTksVec(combinedTrajs.size(), false); 
  OrphanHandle<reco::TrackCollection> combinedTracks = loadTracks(combinedTrajs, event, combTksVec, ttopo);

  OrphanHandle<reco::TrackCollection> trackerTracks;
  std::vector<bool> trackerTksVec(trackerTrajs.size(), false); 
  if(thePutTkTrackFlag) {
    LogTrace(metname) << "Build trackerTracks: "
              << trackerTrajs.size();
    trackerTracks = loadTracks(trackerTrajs, event, trackerTksVec, ttopo, theL2SeededTkLabel, theSmoothTkTrackFlag);
  } else {
    for (TrajectoryContainer::iterator it = trackerTrajs.begin(); it != trackerTrajs.end(); ++it) {
        if(*it) delete *it;
    }
  }

  LogTrace(metname) << "Set the final links in the MuonTrackLinks collection";

  unsigned int candposition(0), position(0), tkposition(0);
  //reco::TrackCollection::const_iterator glIt = combinedTracks->begin(), 
  //  glEnd = combinedTracks->end();

  for (CandidateContainer::const_iterator it = muonCands.begin(); it != muonCands.end(); ++it, ++candposition) {

    // The presence of the global track determines whether to fill the MuonTrackLinks or not
    // N.B. We are assuming here that the global tracks in "combinedTracks" 
    //      have the same order as the muon candidates in "muonCands"
    //      (except for possible missing tracks), which should always be the case...
    //if( glIt == glEnd ) break;
    if(combTksVec[candposition]) {
      reco::TrackRef combinedTR(combinedTracks, position++);
      //++glIt;

      // Create the links between sta and tracker tracks
      reco::MuonTrackLinks links;
      links.setStandAloneTrack((*it)->muonTrack());
      links.setTrackerTrack((*it)->trackerTrack());
      links.setGlobalTrack(combinedTR);

      if(thePutTkTrackFlag && trackerTksVec[candposition]) {
    reco::TrackRef trackerTR(trackerTracks, tkposition++);
    links.setTrackerTrack(trackerTR);
      }

      trackLinksCollection->push_back(links);
    }

    else { // if no global track, still increment the tracker-track counter when appropriate
      if(thePutTkTrackFlag && trackerTksVec[candposition]) tkposition++; 
    }

    delete *it;
  }

  if( thePutTkTrackFlag && trackerTracks.isValid() && !(!combinedTracks->empty() && !trackerTracks->empty() ) )
    LogWarning(metname)<<"The MuonTrackLinkCollection is incomplete"; 
  
  // put the MuonCollection in the event
  LogTrace(metname) << "put the MuonCollection in the event" << "\n";
  
  return event.put(std::move(trackLinksCollection));
} 

OrphanHandle<reco::TrackCollection> 
MuonTrackLoader::loadTracks(const TrajectoryContainer& trajectories,
                Event& event, const std::vector<std::pair<Trajectory*,reco::TrackRef> >& miniMap, Handle<reco::TrackCollection> const& trackHandle, const TrackerTopology& ttopo, const string& instance, bool reallyDoSmoothing) {
  
  const bool doSmoothing = theSmoothingStep && reallyDoSmoothing;
  
  const string metname = "Muon|RecoMuon|MuonTrackLoader|TevMuonTrackLoader";

  LogDebug(metname)<<"TeV LoadTracks instance: " << instance;
  
  // the track collectios; they will be loaded in the event  
  auto trackCollection = std::make_unique<reco::TrackCollection>();
  // ... and its reference into the event
  reco::TrackRefProd trackCollectionRefProd = event.getRefBeforePut<reco::TrackCollection>(instance);
    
  // Association map between GlobalMuons and TeVMuons
  auto trackToTrackmap = std::make_unique<reco::TrackToTrackMap>(trackHandle,trackCollectionRefProd);
   
  // the track extra collection, it will be loaded in the event  
  auto trackExtraCollection = std::make_unique<reco::TrackExtraCollection>();
  // ... and its reference into the event
  reco::TrackExtraRefProd trackExtraCollectionRefProd = event.getRefBeforePut<reco::TrackExtraCollection>(instance);
  
  // the rechit collection, it will be loaded in the event  
  auto recHitCollection = std::make_unique<TrackingRecHitCollection>();
  // ... and its reference into the event
  TrackingRecHitRefProd recHitCollectionRefProd = event.getRefBeforePut<TrackingRecHitCollection>(instance);
  
  // Collection of Trajectory
  auto trajectoryCollection = std::make_unique<std::vector<Trajectory>>();
  
  // don't waste any time...
  if ( trajectories.empty() ) { 
    event.put(std::move(recHitCollection),instance);
    event.put(std::move(trackExtraCollection),instance);
    if(theTrajectoryFlag) {
      event.put(std::move(trajectoryCollection),instance);

      // Association map between track and trajectory
      auto trajTrackMap = std::make_unique<TrajTrackAssociationCollection>();
      event.put(std::move(trajTrackMap), instance );
    }
    event.put(std::move(trackToTrackmap), instance);
    return event.put(std::move(trackCollection),instance);
  }
  
  LogTrace(metname) << "Create the collection of Tracks";
  
  edm::Handle<reco::BeamSpot> beamSpot;
  event.getByToken(theBeamSpotToken,beamSpot);

  reco::TrackRef::key_type trackIndex = 0;
  
  reco::TrackExtraRef::key_type trackExtraIndex = 0;
  
  edm::Ref<reco::TrackCollection>::key_type iTkRef = 0;
  edm::Ref< std::vector<Trajectory> >::key_type iTjRef = 0;
  std::map<unsigned int, unsigned int> tjTkMap;
  
  if(doSmoothing) {
    edm::ESHandle<TrajectorySmoother> aSmoother;
    theService->eventSetup().get<TrajectoryFitter::Record>().get(theSmootherName,aSmoother);
    theSmoother.reset(aSmoother->clone());
    edm::ESHandle<TransientTrackingRecHitBuilder> theTrackerRecHitBuilder;
    theService->eventSetup().get<TransientRecHitRecord>().get(theTrackerRecHitBuilderName,theTrackerRecHitBuilder);
    hitCloner = static_cast<TkTransientTrackingRecHitBuilder const *>(theTrackerRecHitBuilder.product())->cloner();
    theSmoother->setHitCloner(&hitCloner);
  }

  for(TrajectoryContainer::const_iterator rawTrajectory = trajectories.begin();
      rawTrajectory != trajectories.end(); ++rawTrajectory){

    reco::TrackRef glbRef;
    std::vector<std::pair<Trajectory*,reco::TrackRef> >::const_iterator mmit;
    for(mmit = miniMap.begin();mmit!=miniMap.end();++mmit){
      if(mmit->first == *rawTrajectory) glbRef = mmit->second;
    }
    
    Trajectory &trajectory = **rawTrajectory;
    
    if(doSmoothing){
      vector<Trajectory> trajectoriesSM = theSmoother->trajectories(**rawTrajectory);
      
      if(!trajectoriesSM.empty()) {
    const edm::RefToBase<TrajectorySeed> tmpSeedRef = (**rawTrajectory).seedRef();
    trajectory = trajectoriesSM.front();
        trajectory.setSeedRef(tmpSeedRef);
      } else
    LogInfo(metname)<<"The trajectory has not been smoothed!"<<endl;
    }
    
    if(theTrajectoryFlag) {
      trajectoryCollection->push_back(trajectory);
      iTjRef++;
    }    
    
    // build the "bare" track from the trajectory.
    // This track has the parameters defined at PCA (no update)
    pair<bool,reco::Track> resultOfTrackExtrapAtPCA = buildTrackAtPCA(trajectory, *beamSpot);
    
    // Check if the extrapolation went well    
    if(!resultOfTrackExtrapAtPCA.first) {
      delete *rawTrajectory;
      continue;
    }
    
    // take the "bare" track at PCA
    reco::Track &track = resultOfTrackExtrapAtPCA.second;
    
    // build the "bare" track extra from the trajectory
    reco::TrackExtra trackExtra = buildTrackExtra( trajectory );
    
    // get the TrackExtraRef (persitent reference of the track extra)
    reco::TrackExtraRef trackExtraRef(trackExtraCollectionRefProd, trackExtraIndex++ );
    
    // set the persistent track-extra reference to the Track
    track.setExtra(trackExtraRef);

    // Fill the map
    trackToTrackmap->insert(glbRef,
                reco::TrackRef(trackCollectionRefProd,trackIndex++));

    // build the updated-at-vertex track, starting from the previous track
    pair<bool,reco::Track> updateResult(false,reco::Track());
            
    // get the transient rechit and co from the trajectory
    reco::TrackExtra::TrajParams trajParams;
    reco::TrackExtra::Chi2sFive chi2s;
    Traj2TrackHits t2t;
    auto ih = recHitCollection->size();
    t2t(trajectory,*recHitCollection,trajParams,chi2s);
    auto ie = recHitCollection->size();
    // set the TrackingRecHitRef (persitent reference of the tracking rec hits)
    trackExtra.setHits(recHitCollectionRefProd,ih,ie-ih);
    trackExtra.setTrajParams(std::move(trajParams),std::move(chi2s));
    assert(trackExtra.trajParams().size()==trackExtra.recHitsSize());

    // Fill the hit pattern
    for (;ih<ie; ++ih) {
      auto const & hit = (*recHitCollection)[ih];
      auto hits = MuonTrackLoader::unpackHit(hit);
      for (auto hh : hits) {
          if UNLIKELY(!track.appendHitPattern(*hh, ttopo)) break;
      }
 
      if(theUpdatingAtVtx && updateResult.first){
        for (auto hh : hits) {
              if UNLIKELY(!updateResult.second.appendHitPattern(*hh, ttopo)) break;
        }
      }
    }

    // fill the TrackExtraCollection
    trackExtraCollection->push_back(trackExtra);
    
    // fill the TrackCollection
    trackCollection->push_back(track);
    iTkRef++;
    LogTrace(metname) << "Debug Track being loaded pt "<< track.pt();
    
    // We don't need the original trajectory anymore.
    // It has been copied by value in the trajectoryCollection, if 
    // it is required to put it into the event.
    delete *rawTrajectory;

    if(theTrajectoryFlag) tjTkMap[iTjRef-1] = iTkRef-1;
  }
  

  
  // Put the Collections in the event
  LogTrace(metname) << "put the Collections in the event";
  event.put(std::move(recHitCollection),instance);
  event.put(std::move(trackExtraCollection),instance);

  OrphanHandle<reco::TrackCollection> returnTrackHandle;
  OrphanHandle<reco::TrackCollection> nonUpdatedHandle;
  if(theUpdatingAtVtx){
  }
  else {
    event.put(std::move(trackToTrackmap),instance);
    returnTrackHandle = event.put(std::move(trackCollection),instance);
    nonUpdatedHandle = returnTrackHandle;
  }

  if ( theTrajectoryFlag ) {
    OrphanHandle<std::vector<Trajectory> > rTrajs = event.put(std::move(trajectoryCollection),instance);

    // Association map between track and trajectory
    auto trajTrackMap = std::make_unique<TrajTrackAssociationCollection>(rTrajs, nonUpdatedHandle);

    // Now Create traj<->tracks association map
    for ( std::map<unsigned int, unsigned int>::iterator i = tjTkMap.begin(); 
          i != tjTkMap.end(); i++ ) {
      trajTrackMap->insert( edm::Ref<std::vector<Trajectory> >( rTrajs, (*i).first ),
                            edm::Ref<reco::TrackCollection>( nonUpdatedHandle, (*i).second ) );
    }
    event.put(std::move(trajTrackMap), instance );
  }

  return returnTrackHandle;
}


pair<bool,reco::Track> MuonTrackLoader::buildTrackAtPCA(const Trajectory& trajectory, const reco::BeamSpot &beamSpot) const {

  const string metname = "Muon|RecoMuon|MuonTrackLoader";

  MuonPatternRecoDumper debug;
  
  // FIXME: check the prop direction
  TrajectoryStateOnSurface innerTSOS = trajectory.geometricalInnermostState();
  
  // This is needed to extrapolate the tsos at vertex
  LogTrace(metname) << "Propagate to PCA...";
  pair<bool,FreeTrajectoryState> 
    extrapolationResult = theUpdatorAtVtx->propagate(innerTSOS, beamSpot);  
  FreeTrajectoryState ftsAtVtx;
  
  if(extrapolationResult.first)
    ftsAtVtx = extrapolationResult.second;
  else{    
    if(TrackerBounds::isInside(innerTSOS.globalPosition())){
      LogInfo(metname) << "Track in the Tracker: taking the innermost state instead of the state at PCA";
      ftsAtVtx = *innerTSOS.freeState();
    }
    else{
      if ( theAllowNoVtxFlag ) {
        LogInfo(metname) << "Propagation to PCA failed, taking the innermost state instead of the state at PCA";
        ftsAtVtx = *innerTSOS.freeState();
      } else {
        LogInfo(metname) << "Stand Alone track: this track will be rejected";
        return pair<bool,reco::Track>(false,reco::Track());
      }
    }
  }
    
  LogTrace(metname) << "TSOS after the extrapolation at vtx";
  LogTrace(metname) << debug.dumpFTS(ftsAtVtx);
  
  GlobalPoint pca = ftsAtVtx.position();
  math::XYZPoint persistentPCA(pca.x(),pca.y(),pca.z());
  GlobalVector p = ftsAtVtx.momentum();
  math::XYZVector persistentMomentum(p.x(),p.y(),p.z());

  bool bon = true;
  if(fabs(theService->magneticField()->inTesla(GlobalPoint(0,0,0)).z()) < 0.01) bon=false;   
  double ndof = trajectory.ndof(bon);
  
  reco::Track track(trajectory.chiSquared(), 
            ndof,
            persistentPCA,
            persistentMomentum,
            ftsAtVtx.charge(),
            ftsAtVtx.curvilinearError());
  
  return pair<bool,reco::Track>(true,track);
}


pair<bool,reco::Track> MuonTrackLoader::buildTrackUpdatedAtPCA(const reco::Track &track, const reco::BeamSpot &beamSpot) const {

  const string metname = "Muon|RecoMuon|MuonTrackLoader";
  MuonPatternRecoDumper debug;
 
  // build the transient track
  reco::TransientTrack transientTrack(track,
                      &*theService->magneticField(),
                      theService->trackingGeometry());

  LogTrace(metname) << "Apply the vertex constraint";
  pair<bool,FreeTrajectoryState> updateResult = theUpdatorAtVtx->update(transientTrack,beamSpot);

  if(!updateResult.first){
    return pair<bool,reco::Track>(false,reco::Track());
  }

  LogTrace(metname) << "FTS after the vertex constraint";
  FreeTrajectoryState &ftsAtVtx = updateResult.second;

  LogTrace(metname) << debug.dumpFTS(ftsAtVtx);
  
  GlobalPoint pca = ftsAtVtx.position();
  math::XYZPoint persistentPCA(pca.x(),pca.y(),pca.z());
  GlobalVector p = ftsAtVtx.momentum();
  math::XYZVector persistentMomentum(p.x(),p.y(),p.z());
  
  reco::Track updatedTrack(track.chi2(), 
               track.ndof(),
               persistentPCA,
               persistentMomentum,
               ftsAtVtx.charge(),
               ftsAtVtx.curvilinearError());
  
  return pair<bool,reco::Track>(true,updatedTrack);
}


reco::TrackExtra MuonTrackLoader::buildTrackExtra(const Trajectory& trajectory) const {

  const string metname = "Muon|RecoMuon|MuonTrackLoader";

  const Trajectory::RecHitContainer transRecHits = trajectory.recHits();
  
  // put the collection of TrackingRecHit in the event
  
  // sets the outermost and innermost TSOSs
  // FIXME: check it!
  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 = theService->trackingGeometry()->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;
 
}