TrackingTools/src/MuonTrackLoader.cc

0001 /** \class MuonTrackLoader
0002  *  Class to load the product in the event
0003  *
0004
0005
0006  *  \author R. Bellan - INFN Torino <riccardo.bellan@cern.ch>
0007  */
0008
0009 #include <memory>
0010
0011 #include "RecoMuon/TrackingTools/interface/MuonTrackLoader.h"
0012
0013 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0014 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0015 #include "RecoMuon/TrackingTools/interface/MuonUpdatorAtVertex.h"
0016
0017 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0018 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0019 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0020 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0021 #include "TrackingTools/GeomPropagators/interface/TrackerBounds.h"
0022 #include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
0023 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0024 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0025 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0026
0027 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0032
0033 #include "DataFormats/Math/interface/LorentzVector.h"
0034 #include "DataFormats/TrackReco/interface/Track.h"
0035 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0036 #include "DataFormats/TrackReco/interface/TrackToTrackMap.h"
0037
0038 #include "DataFormats/MuonReco/interface/MuonTrackLinks.h"
0039 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0040
0041 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0042 #include "RecoTracker/TransientTrackingRecHit/interface/Traj2TrackHits.h"
0043
0044 using namespace edm;
0045 using namespace std;
0046 using namespace reco;
0047
0048 std::vector<const TrackingRecHit*> MuonTrackLoader::unpackHit(const TrackingRecHit& hit) {
0049   // get rec hit det id and rec hit type
0050   DetId id = hit.geographicalId();
0051   uint16_t detid = id.det();
0052
0053   std::vector<const TrackingRecHit*> hits;
0054
0055   if (detid == DetId::Tracker) {
0056     hits.push_back(&hit);
0057   } else if (detid == DetId::Muon) {
0058     uint16_t subdet = id.subdetId();
0059     if (subdet == (uint16_t)MuonSubdetId::DT) {
0060       if (hit.dimension() == 1) {  // DT rechit (granularity 2)
0061         hits.push_back(&hit);
0062       } else if (hit.dimension() > 1) {  // 4D segment (granularity 0).
0063         // Both 2 and 4 dim cases. MB4s have 2D, but formatted in 4D segment
0064         // 4D --> 2D
0065         std::vector<const TrackingRecHit*> seg2D = hit.recHits();
0066         // load 1D hits (2D --> 1D)
0067         for (std::vector<const TrackingRecHit*>::const_iterator it = seg2D.begin(); it != seg2D.end(); ++it) {
0068           std::vector<const TrackingRecHit*> hits1D = (*it)->recHits();
0069           copy(hits1D.begin(), hits1D.end(), back_inserter(hits));
0070         }
0071       }
0072     } else if (subdet == (uint16_t)MuonSubdetId::CSC) {
0073       if (hit.dimension() == 2) {  // CSC rechit (granularity 2)
0074         hits.push_back(&hit);
0075       } else if (hit.dimension() == 4) {  // 4D segment (granularity 0)
0076         // load 2D hits (4D --> 1D)
0077         hits = hit.recHits();
0078       }
0079     } else if (subdet == (uint16_t)MuonSubdetId::RPC) {
0080       hits.push_back(&hit);
0081     } else if (subdet == (uint16_t)MuonSubdetId::GEM) {
0082       if (hit.dimension() == 2) {  // GEM rechit
0083         hits.push_back(&hit);
0084       } else if (hit.dimension() == 4) {  // GEM segment
0085         hits = hit.recHits();
0086       }
0087     } else if (subdet == (uint16_t)MuonSubdetId::ME0) {  //segment
0088       hits = hit.recHits();
0089     }
0090   }
0091   return hits;
0092 }
0093
0094 // constructor (obsolete, should use eventSetUp not service..)
0095 MuonTrackLoader::MuonTrackLoader(ParameterSet& parameterSet,
0096                                  edm::ConsumesCollector& iC,
0097                                  const MuonServiceProxy* service)
0098     : theService(service) {
0099   // option to do or not the smoothing step.
0100   // the trajectories which are passed to the track loader are supposed to be non-smoothed
0101   theSmoothingStep = parameterSet.getParameter<bool>("DoSmoothing");
0102   if (theSmoothingStep) {
0103     auto smootherName = parameterSet.getParameter<string>("Smoother");
0104     theSmootherToken = iC.esConsumes(edm::ESInputTag("", smootherName));
0105
0106     auto trackerRecHitBuilderName = parameterSet.getParameter<std::string>("TTRHBuilder");
0107     theTrackerRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", trackerRecHitBuilderName));
0108   }
0109
0110   // update at vertex
0111   theUpdatingAtVtx = parameterSet.getParameter<bool>("VertexConstraint");
0112
0113   // beam spot input tag
0114   theBeamSpotInputTag = parameterSet.getParameter<edm::InputTag>("beamSpot");
0115   theBeamSpotToken = iC.consumes<reco::BeamSpot>(theBeamSpotInputTag);
0116
0117   // Flag to put the trajectory into the event
0118   theTrajectoryFlag = parameterSet.getUntrackedParameter<bool>("PutTrajectoryIntoEvent", true);
0119
0120   theL2SeededTkLabel = parameterSet.getUntrackedParameter<string>("MuonSeededTracksInstance", string());
0121
0122   ParameterSet updatorPar = parameterSet.getParameter<ParameterSet>("MuonUpdatorAtVertexParameters");
0123   theUpdatorAtVtx = std::make_unique<MuonUpdatorAtVertex>(updatorPar, service);
0124
0125   thePutTkTrackFlag = parameterSet.getUntrackedParameter<bool>("PutTkTrackIntoEvent", false);
0126   theSmoothTkTrackFlag = parameterSet.getUntrackedParameter<bool>("SmoothTkTrack", false);
0127   theAllowNoVtxFlag = parameterSet.getUntrackedParameter<bool>("AllowNoVertex", false);
0128 }
0129
0130 MuonTrackLoader::~MuonTrackLoader() {}
0131
0132 OrphanHandle<reco::TrackCollection> MuonTrackLoader::loadTracks(TrajectoryContainer& trajectories,
0133                                                                 Event& event,
0134                                                                 const TrackerTopology& ttopo,
0135                                                                 const string& instance,
0136                                                                 bool reallyDoSmoothing) {
0137   std::vector<bool> dummyVecBool;
0138   return loadTracks(trajectories, event, dummyVecBool, ttopo, instance, reallyDoSmoothing);
0139 }
0140
0141 OrphanHandle<reco::TrackCollection> MuonTrackLoader::loadTracks(TrajectoryContainer& trajectories,
0142                                                                 Event& event,
0143                                                                 std::vector<bool>& tkBoolVec,
0144                                                                 const TrackerTopology& ttopo,
0145                                                                 const string& instance,
0146                                                                 bool reallyDoSmoothing) {
0147   const bool doSmoothing = theSmoothingStep && reallyDoSmoothing;
0148
0149   const string metname = "Muon|RecoMuon|MuonTrackLoader";
0150
0151   // the track collectios; they will be loaded in the event
0152   auto trackCollection = std::make_unique<reco::TrackCollection>();
0153   // ... and its reference into the event
0154   reco::TrackRefProd trackCollectionRefProd = event.getRefBeforePut<reco::TrackCollection>(instance);
0155
0156   // track collection for the tracks updated at vertex
0157   auto updatedAtVtxTrackCollection = std::make_unique<reco::TrackCollection>();
0158   // ... and its (eventually) reference into the event
0159   reco::TrackRefProd trackUpdatedCollectionRefProd;
0160   if (theUpdatingAtVtx)
0161     trackUpdatedCollectionRefProd = event.getRefBeforePut<reco::TrackCollection>(instance + "UpdatedAtVtx");
0162
0163   // Association map between updated and non updated at vtx tracks
0164   auto trackToTrackmap = std::make_unique<reco::TrackToTrackMap>(trackCollectionRefProd, trackUpdatedCollectionRefProd);
0165
0166   // the track extra collection, it will be loaded in the event
0167   auto trackExtraCollection = std::make_unique<reco::TrackExtraCollection>();
0168   // ... and its reference into the event
0169   reco::TrackExtraRefProd trackExtraCollectionRefProd = event.getRefBeforePut<reco::TrackExtraCollection>(instance);
0170
0171   // the rechit collection, it will be loaded in the event
0172   auto recHitCollection = std::make_unique<TrackingRecHitCollection>();
0173   // ... and its reference into the event
0174   TrackingRecHitRefProd recHitCollectionRefProd = event.getRefBeforePut<TrackingRecHitCollection>(instance);
0175
0176   // Collection of Trajectory
0177   auto trajectoryCollection = std::make_unique<vector<Trajectory>>();
0178
0179   // don't waste any time...
0180   if (trajectories.empty()) {
0181     event.put(std::move(recHitCollection), instance);
0182     event.put(std::move(trackExtraCollection), instance);
0183     if (theTrajectoryFlag) {
0184       event.put(std::move(trajectoryCollection), instance);
0185
0186       // Association map between track and trajectory
0187       auto trajTrackMap = std::make_unique<TrajTrackAssociationCollection>();
0188       event.put(std::move(trajTrackMap), instance);
0189     }
0190     if (theUpdatingAtVtx) {
0191       event.put(std::move(trackToTrackmap));
0192       event.put(std::move(updatedAtVtxTrackCollection), instance + "UpdatedAtVtx");
0193     }
0194     return event.put(std::move(trackCollection), instance);
0195   }
0196
0197   edm::Handle<reco::BeamSpot> beamSpot;
0198   event.getByToken(theBeamSpotToken, beamSpot);
0199
0200   LogTrace(metname) << "Create the collection of Tracks";
0201
0202   reco::TrackRef::key_type trackIndex = 0;
0203   reco::TrackRef::key_type trackUpdatedIndex = 0;
0204
0205   reco::TrackExtraRef::key_type trackExtraIndex = 0;
0206
0207   edm::Ref<reco::TrackCollection>::key_type iTkRef = 0;
0208   edm::Ref<std::vector<Trajectory>>::key_type iTjRef = 0;
0209   std::map<unsigned int, unsigned int> tjTkMap;
0210
0211   if (doSmoothing) {
0212     TrajectorySmoother const& aSmoother = theService->eventSetup().getData(theSmootherToken);
0213     theSmoother.reset(aSmoother.clone());
0214     TransientTrackingRecHitBuilder const& theTrackerRecHitBuilder =
0215         theService->eventSetup().getData(theTrackerRecHitBuilderToken);
0216     hitCloner = static_cast<TkTransientTrackingRecHitBuilder const&>(theTrackerRecHitBuilder).cloner();
0217     theSmoother->setHitCloner(&hitCloner);
0218   }
0219
0220   unsigned int tjCnt = 0;
0221   for (TrajectoryContainer::iterator itRawTrajectory = trajectories.begin(); itRawTrajectory != trajectories.end();
0222        ++itRawTrajectory, ++tjCnt) {
0223     auto rawTrajectory = std::move(*itRawTrajectory);
0224     Trajectory& trajectory = *rawTrajectory;
0225
0226     if (doSmoothing) {
0227       vector<Trajectory> trajectoriesSM = theSmoother->trajectories(*rawTrajectory);
0228
0229       if (!trajectoriesSM.empty()) {
0230         const edm::RefToBase<TrajectorySeed> tmpSeedRef = (*rawTrajectory).seedRef();
0231         trajectory = trajectoriesSM.front();
0232         trajectory.setSeedRef(tmpSeedRef);
0233         LogDebug(metname) << "theSeedRef.isNonnull " << trajectory.seedRef().isNonnull();
0234       } else
0235         LogInfo(metname) << "The trajectory has not been smoothed!" << endl;
0236     }
0237
0238     if (theTrajectoryFlag) {
0239       trajectoryCollection->push_back(trajectory);
0240       iTjRef++;
0241     }
0242
0243     // build the "bare" track from the trajectory.
0244     // This track has the parameters defined at PCA (no update)
0245     pair<bool, reco::Track> resultOfTrackExtrapAtPCA = buildTrackAtPCA(trajectory, *beamSpot);
0246
0247     // Check if the extrapolation went well
0248     if (!resultOfTrackExtrapAtPCA.first) {
0249       continue;
0250     }
0251
0252     // take the "bare" track at PCA
0253     reco::Track& track = resultOfTrackExtrapAtPCA.second;
0254
0255     // build the "bare" track extra from the trajectory
0256     reco::TrackExtra trackExtra = buildTrackExtra(trajectory);
0257
0258     // get the TrackExtraRef (persitent reference of the track extra)
0259     reco::TrackExtraRef trackExtraRef(trackExtraCollectionRefProd, trackExtraIndex++);
0260
0261     // set the persistent track-extra reference to the Track
0262     track.setExtra(trackExtraRef);
0263
0264     // build the updated-at-vertex track, starting from the previous track
0265     pair<bool, reco::Track> updateResult(false, reco::Track());
0266
0267     if (theUpdatingAtVtx) {
0268       // build the "bare" track UPDATED at vtx
0269       updateResult = buildTrackUpdatedAtPCA(track, *beamSpot);
0270
0271       if (!updateResult.first)
0272         ++trackIndex;
0273       else {
0274         // set the persistent track-extra reference to the Track
0275         updateResult.second.setExtra(trackExtraRef);
0276
0277         // Fill the map
0278         trackToTrackmap->insert(reco::TrackRef(trackCollectionRefProd, trackIndex++),
0279                                 reco::TrackRef(trackUpdatedCollectionRefProd, trackUpdatedIndex++));
0280       }
0281     }
0282
0283     // get the transient rechit and co from the trajectory
0284     reco::TrackExtra::TrajParams trajParams;
0285     reco::TrackExtra::Chi2sFive chi2s;
0286     Traj2TrackHits t2t;
0287     auto ih = recHitCollection->size();
0288     t2t(trajectory, *recHitCollection, trajParams, chi2s);
0289     auto ie = recHitCollection->size();
0290     // set the TrackingRecHitRef (persitent reference of the tracking rec hits)
0291     trackExtra.setHits(recHitCollectionRefProd, ih, ie - ih);
0292     trackExtra.setTrajParams(std::move(trajParams), std::move(chi2s));
0293     assert(trackExtra.trajParams().size() == trackExtra.recHitsSize());
0294
0295     // Fill the hit pattern
0296     for (; ih < ie; ++ih) {
0297       auto const& hit = (*recHitCollection)[ih];
0298       auto hits = MuonTrackLoader::unpackHit(hit);
0299       for (auto hh : hits) {
0300         if UNLIKELY (!track.appendHitPattern(*hh, ttopo))
0301           break;
0302       }
0303
0304       if (theUpdatingAtVtx && updateResult.first) {
0305         for (auto hh : hits) {
0306           if UNLIKELY (!updateResult.second.appendHitPattern(*hh, ttopo))
0307             break;
0308         }
0309       }
0310     }
0311
0312     // fill the TrackExtraCollection
0313     trackExtraCollection->push_back(trackExtra);
0314
0315     // fill the TrackCollection
0316     trackCollection->push_back(track);
0317     iTkRef++;
0318     LogTrace(metname) << "Debug Track being loaded pt " << track.pt();
0319     // fill the TrackCollection updated at vtx
0320     if (theUpdatingAtVtx && updateResult.first)
0321       updatedAtVtxTrackCollection->push_back(updateResult.second);
0322
0323     if (tkBoolVec.size() > tjCnt)
0324       tkBoolVec[tjCnt] = true;
0325     if (theTrajectoryFlag)
0326       tjTkMap[iTjRef - 1] = iTkRef - 1;
0327   }
0328
0329   // Put the Collections in the event
0330   LogTrace(metname) << "put the Collections in the event";
0331   event.put(std::move(recHitCollection), instance);
0332   event.put(std::move(trackExtraCollection), instance);
0333
0334   OrphanHandle<reco::TrackCollection> returnTrackHandle;
0335   OrphanHandle<reco::TrackCollection> nonUpdatedHandle;
0336   if (theUpdatingAtVtx) {
0337     nonUpdatedHandle = event.put(std::move(trackCollection), instance);
0338     event.put(std::move(trackToTrackmap));
0339     returnTrackHandle = event.put(std::move(updatedAtVtxTrackCollection), instance + "UpdatedAtVtx");
0340   } else {
0341     returnTrackHandle = event.put(std::move(trackCollection), instance);
0342     nonUpdatedHandle = returnTrackHandle;
0343   }
0344
0345   if (theTrajectoryFlag) {
0346     OrphanHandle<std::vector<Trajectory>> rTrajs = event.put(std::move(trajectoryCollection), instance);
0347
0348     // Association map between track and trajectory
0349     auto trajTrackMap = std::make_unique<TrajTrackAssociationCollection>(rTrajs, nonUpdatedHandle);
0350
0351     // Now Create traj<->tracks association map
0352     for (std::map<unsigned int, unsigned int>::iterator i = tjTkMap.begin(); i != tjTkMap.end(); i++) {
0353       trajTrackMap->insert(edm::Ref<std::vector<Trajectory>>(rTrajs, (*i).first),
0354                            edm::Ref<reco::TrackCollection>(nonUpdatedHandle, (*i).second));
0355     }
0356     event.put(std::move(trajTrackMap), instance);
0357   }
0358
0359   return returnTrackHandle;
0360 }
0361
0362 OrphanHandle<reco::MuonTrackLinksCollection> MuonTrackLoader::loadTracks(CandidateContainer& muonCands,
0363                                                                          Event& event,
0364                                                                          const TrackerTopology& ttopo) {
0365   const string metname = "Muon|RecoMuon|MuonTrackLoader";
0366
0367   // the muon collection, it will be loaded in the event
0368   auto trackLinksCollection = std::make_unique<reco::MuonTrackLinksCollection>();
0369
0370   // don't waste any time...
0371   if (muonCands.empty()) {
0372     auto trackExtraCollection = std::make_unique<reco::TrackExtraCollection>();
0373     auto recHitCollection = std::make_unique<TrackingRecHitCollection>();
0374     auto trackCollection = std::make_unique<reco::TrackCollection>();
0375
0376     event.put(std::move(recHitCollection));
0377     event.put(std::move(trackExtraCollection));
0378     event.put(std::move(trackCollection));
0379
0380     //need to also put the tracker tracks collection if requested
0381     if (thePutTkTrackFlag) {
0382       //will take care of putting nothing in the event but the empty collection
0383       TrajectoryContainer trackerTrajs;
0384       loadTracks(trackerTrajs, event, ttopo, theL2SeededTkLabel, theSmoothTkTrackFlag);
0385     }
0386
0387     return event.put(std::move(trackLinksCollection));
0388   }
0389
0390   // get combined Trajectories
0391   TrajectoryContainer combinedTrajs;
0392   TrajectoryContainer trackerTrajs;
0393   for (CandidateContainer::iterator it = muonCands.begin(); it != muonCands.end(); ++it) {
0394     LogDebug(metname) << "Loader glbSeedRef " << (*it)->trajectory()->seedRef().isNonnull();
0395     if ((*it)->trackerTrajectory())
0396       LogDebug(metname) << " "
0397                         << "tkSeedRef " << (*it)->trackerTrajectory()->seedRef().isNonnull();
0398
0399     combinedTrajs.push_back((*it)->releaseTrajectory());
0400     {
0401       auto tt = (*it)->releaseTrackerTrajectory();
0402       if (thePutTkTrackFlag)
0403         trackerTrajs.push_back(std::move(tt));
0404     }
0405
0406     // // Create the links between sta and tracker tracks
0407     // reco::MuonTrackLinks links;
0408     // links.setStandAloneTrack((*it)->muonTrack());
0409     // links.setTrackerTrack((*it)->trackerTrack());
0410     // trackLinksCollection->push_back(links);
0411     // delete *it;
0412   }
0413
0414   // create the TrackCollection of combined Trajectories
0415   // FIXME: could this be done one track at a time in the previous loop?
0416   LogTrace(metname) << "Build combinedTracks";
0417   std::vector<bool> combTksVec(combinedTrajs.size(), false);
0418   OrphanHandle<reco::TrackCollection> combinedTracks = loadTracks(combinedTrajs, event, combTksVec, ttopo);
0419
0420   OrphanHandle<reco::TrackCollection> trackerTracks;
0421   std::vector<bool> trackerTksVec(trackerTrajs.size(), false);
0422   if (thePutTkTrackFlag) {
0423     LogTrace(metname) << "Build trackerTracks: " << trackerTrajs.size();
0424     trackerTracks = loadTracks(trackerTrajs, event, trackerTksVec, ttopo, theL2SeededTkLabel, theSmoothTkTrackFlag);
0425   }
0426
0427   trackerTrajs.clear();
0428
0429   LogTrace(metname) << "Set the final links in the MuonTrackLinks collection";
0430
0431   unsigned int candposition(0), position(0), tkposition(0);
0432   //reco::TrackCollection::const_iterator glIt = combinedTracks->begin(),
0433   //  glEnd = combinedTracks->end();
0434
0435   for (CandidateContainer::const_iterator it = muonCands.begin(); it != muonCands.end(); ++it, ++candposition) {
0436     // The presence of the global track determines whether to fill the MuonTrackLinks or not
0437     // N.B. We are assuming here that the global tracks in "combinedTracks"
0438     //      have the same order as the muon candidates in "muonCands"
0439     //      (except for possible missing tracks), which should always be the case...
0440     //if( glIt == glEnd ) break;
0441     if (combTksVec[candposition]) {
0442       reco::TrackRef combinedTR(combinedTracks, position++);
0443       //++glIt;
0444
0445       // Create the links between sta and tracker tracks
0446       reco::MuonTrackLinks links;
0447       links.setStandAloneTrack((*it)->muonTrack());
0448       links.setTrackerTrack((*it)->trackerTrack());
0449       links.setGlobalTrack(combinedTR);
0450
0451       if (thePutTkTrackFlag && trackerTksVec[candposition]) {
0452         reco::TrackRef trackerTR(trackerTracks, tkposition++);
0453         links.setTrackerTrack(trackerTR);
0454       }
0455
0456       trackLinksCollection->push_back(links);
0457     }
0458
0459     else {  // if no global track, still increment the tracker-track counter when appropriate
0460       if (thePutTkTrackFlag && trackerTksVec[candposition])
0461         tkposition++;
0462     }
0463   }
0464
0465   if (thePutTkTrackFlag && trackerTracks.isValid() && !(!combinedTracks->empty() && !trackerTracks->empty()))
0466     LogWarning(metname) << "The MuonTrackLinkCollection is incomplete";
0467
0468   // put the MuonCollection in the event
0469   LogTrace(metname) << "put the MuonCollection in the event"
0470                     << "\n";
0471
0472   return event.put(std::move(trackLinksCollection));
0473 }
0474
0475 OrphanHandle<reco::TrackCollection> MuonTrackLoader::loadTracks(
0476     TrajectoryContainer& trajectories,
0477     Event& event,
0478     const std::vector<std::pair<Trajectory*, reco::TrackRef>>& miniMap,
0479     Handle<reco::TrackCollection> const& trackHandle,
0480     const TrackerTopology& ttopo,
0481     const string& instance,
0482     bool reallyDoSmoothing) {
0483   const bool doSmoothing = theSmoothingStep && reallyDoSmoothing;
0484
0485   const string metname = "Muon|RecoMuon|MuonTrackLoader|TevMuonTrackLoader";
0486
0487   LogDebug(metname) << "TeV LoadTracks instance: " << instance;
0488
0489   // the track collectios; they will be loaded in the event
0490   auto trackCollection = std::make_unique<reco::TrackCollection>();
0491   // ... and its reference into the event
0492   reco::TrackRefProd trackCollectionRefProd = event.getRefBeforePut<reco::TrackCollection>(instance);
0493
0494   // Association map between GlobalMuons and TeVMuons
0495   auto trackToTrackmap = std::make_unique<reco::TrackToTrackMap>(trackHandle, trackCollectionRefProd);
0496
0497   // the track extra collection, it will be loaded in the event
0498   auto trackExtraCollection = std::make_unique<reco::TrackExtraCollection>();
0499   // ... and its reference into the event
0500   reco::TrackExtraRefProd trackExtraCollectionRefProd = event.getRefBeforePut<reco::TrackExtraCollection>(instance);
0501
0502   // the rechit collection, it will be loaded in the event
0503   auto recHitCollection = std::make_unique<TrackingRecHitCollection>();
0504   // ... and its reference into the event
0505   TrackingRecHitRefProd recHitCollectionRefProd = event.getRefBeforePut<TrackingRecHitCollection>(instance);
0506
0507   // Collection of Trajectory
0508   auto trajectoryCollection = std::make_unique<std::vector<Trajectory>>();
0509
0510   // don't waste any time...
0511   if (trajectories.empty()) {
0512     event.put(std::move(recHitCollection), instance);
0513     event.put(std::move(trackExtraCollection), instance);
0514     if (theTrajectoryFlag) {
0515       event.put(std::move(trajectoryCollection), instance);
0516
0517       // Association map between track and trajectory
0518       auto trajTrackMap = std::make_unique<TrajTrackAssociationCollection>();
0519       event.put(std::move(trajTrackMap), instance);
0520     }
0521     event.put(std::move(trackToTrackmap), instance);
0522     return event.put(std::move(trackCollection), instance);
0523   }
0524
0525   LogTrace(metname) << "Create the collection of Tracks";
0526
0527   edm::Handle<reco::BeamSpot> beamSpot;
0528   event.getByToken(theBeamSpotToken, beamSpot);
0529
0530   reco::TrackRef::key_type trackIndex = 0;
0531
0532   reco::TrackExtraRef::key_type trackExtraIndex = 0;
0533
0534   edm::Ref<reco::TrackCollection>::key_type iTkRef = 0;
0535   edm::Ref<std::vector<Trajectory>>::key_type iTjRef = 0;
0536   std::map<unsigned int, unsigned int> tjTkMap;
0537
0538   if (doSmoothing) {
0539     TrajectorySmoother const& aSmoother = theService->eventSetup().getData(theSmootherToken);
0540     theSmoother.reset(aSmoother.clone());
0541     TransientTrackingRecHitBuilder const& theTrackerRecHitBuilder =
0542         theService->eventSetup().getData(theTrackerRecHitBuilderToken);
0543     hitCloner = static_cast<TkTransientTrackingRecHitBuilder const&>(theTrackerRecHitBuilder).cloner();
0544     theSmoother->setHitCloner(&hitCloner);
0545   }
0546
0547   for (TrajectoryContainer::iterator itRawTrajectory = trajectories.begin(); itRawTrajectory != trajectories.end();
0548        ++itRawTrajectory) {
0549     auto rawTrajectory = std::move(*itRawTrajectory);
0550     reco::TrackRef glbRef;
0551     std::vector<std::pair<Trajectory*, reco::TrackRef>>::const_iterator mmit;
0552     for (mmit = miniMap.begin(); mmit != miniMap.end(); ++mmit) {
0553       if (mmit->first == rawTrajectory.get())
0554         glbRef = mmit->second;
0555     }
0556
0557     Trajectory& trajectory = *rawTrajectory;
0558
0559     if (doSmoothing) {
0560       vector<Trajectory> trajectoriesSM = theSmoother->trajectories(*rawTrajectory);
0561
0562       if (!trajectoriesSM.empty()) {
0563         const edm::RefToBase<TrajectorySeed> tmpSeedRef = (*rawTrajectory).seedRef();
0564         trajectory = trajectoriesSM.front();
0565         trajectory.setSeedRef(tmpSeedRef);
0566       } else
0567         LogInfo(metname) << "The trajectory has not been smoothed!" << endl;
0568     }
0569
0570     if (theTrajectoryFlag) {
0571       trajectoryCollection->push_back(trajectory);
0572       iTjRef++;
0573     }
0574
0575     // build the "bare" track from the trajectory.
0576     // This track has the parameters defined at PCA (no update)
0577     pair<bool, reco::Track> resultOfTrackExtrapAtPCA = buildTrackAtPCA(trajectory, *beamSpot);
0578
0579     // Check if the extrapolation went well
0580     if (!resultOfTrackExtrapAtPCA.first) {
0581       continue;
0582     }
0583
0584     // take the "bare" track at PCA
0585     reco::Track& track = resultOfTrackExtrapAtPCA.second;
0586
0587     // build the "bare" track extra from the trajectory
0588     reco::TrackExtra trackExtra = buildTrackExtra(trajectory);
0589
0590     // get the TrackExtraRef (persitent reference of the track extra)
0591     reco::TrackExtraRef trackExtraRef(trackExtraCollectionRefProd, trackExtraIndex++);
0592
0593     // set the persistent track-extra reference to the Track
0594     track.setExtra(trackExtraRef);
0595
0596     // Fill the map
0597     trackToTrackmap->insert(glbRef, reco::TrackRef(trackCollectionRefProd, trackIndex++));
0598
0599     // build the updated-at-vertex track, starting from the previous track
0600     pair<bool, reco::Track> updateResult(false, reco::Track());
0601
0602     // get the transient rechit and co from the trajectory
0603     reco::TrackExtra::TrajParams trajParams;
0604     reco::TrackExtra::Chi2sFive chi2s;
0605     Traj2TrackHits t2t;
0606     auto ih = recHitCollection->size();
0607     t2t(trajectory, *recHitCollection, trajParams, chi2s);
0608     auto ie = recHitCollection->size();
0609     // set the TrackingRecHitRef (persitent reference of the tracking rec hits)
0610     trackExtra.setHits(recHitCollectionRefProd, ih, ie - ih);
0611     trackExtra.setTrajParams(std::move(trajParams), std::move(chi2s));
0612     assert(trackExtra.trajParams().size() == trackExtra.recHitsSize());
0613
0614     // Fill the hit pattern
0615     for (; ih < ie; ++ih) {
0616       auto const& hit = (*recHitCollection)[ih];
0617       auto hits = MuonTrackLoader::unpackHit(hit);
0618       for (auto hh : hits) {
0619         if UNLIKELY (!track.appendHitPattern(*hh, ttopo))
0620           break;
0621       }
0622
0623       if (theUpdatingAtVtx && updateResult.first) {
0624         for (auto hh : hits) {
0625           if UNLIKELY (!updateResult.second.appendHitPattern(*hh, ttopo))
0626             break;
0627         }
0628       }
0629     }
0630
0631     // fill the TrackExtraCollection
0632     trackExtraCollection->push_back(trackExtra);
0633
0634     // fill the TrackCollection
0635     trackCollection->push_back(track);
0636     iTkRef++;
0637     LogTrace(metname) << "Debug Track being loaded pt " << track.pt();
0638
0639     if (theTrajectoryFlag)
0640       tjTkMap[iTjRef - 1] = iTkRef - 1;
0641   }
0642
0643   // Put the Collections in the event
0644   LogTrace(metname) << "put the Collections in the event";
0645   event.put(std::move(recHitCollection), instance);
0646   event.put(std::move(trackExtraCollection), instance);
0647
0648   OrphanHandle<reco::TrackCollection> returnTrackHandle;
0649   OrphanHandle<reco::TrackCollection> nonUpdatedHandle;
0650   if (theUpdatingAtVtx) {
0651   } else {
0652     event.put(std::move(trackToTrackmap), instance);
0653     returnTrackHandle = event.put(std::move(trackCollection), instance);
0654     nonUpdatedHandle = returnTrackHandle;
0655   }
0656
0657   if (theTrajectoryFlag) {
0658     OrphanHandle<std::vector<Trajectory>> rTrajs = event.put(std::move(trajectoryCollection), instance);
0659
0660     // Association map between track and trajectory
0661     auto trajTrackMap = std::make_unique<TrajTrackAssociationCollection>(rTrajs, nonUpdatedHandle);
0662
0663     // Now Create traj<->tracks association map
0664     for (std::map<unsigned int, unsigned int>::iterator i = tjTkMap.begin(); i != tjTkMap.end(); i++) {
0665       trajTrackMap->insert(edm::Ref<std::vector<Trajectory>>(rTrajs, (*i).first),
0666                            edm::Ref<reco::TrackCollection>(nonUpdatedHandle, (*i).second));
0667     }
0668     event.put(std::move(trajTrackMap), instance);
0669   }
0670
0671   return returnTrackHandle;
0672 }
0673
0674 pair<bool, reco::Track> MuonTrackLoader::buildTrackAtPCA(const Trajectory& trajectory,
0675                                                          const reco::BeamSpot& beamSpot) const {
0676   const string metname = "Muon|RecoMuon|MuonTrackLoader";
0677
0678   MuonPatternRecoDumper debug;
0679
0680   // FIXME: check the prop direction
0681   TrajectoryStateOnSurface innerTSOS = trajectory.geometricalInnermostState();
0682
0683   // This is needed to extrapolate the tsos at vertex
0684   LogTrace(metname) << "Propagate to PCA...";
0685   pair<bool, FreeTrajectoryState> extrapolationResult = theUpdatorAtVtx->propagate(innerTSOS, beamSpot);
0686   FreeTrajectoryState ftsAtVtx;
0687
0688   if (extrapolationResult.first)
0689     ftsAtVtx = extrapolationResult.second;
0690   else {
0691     if (TrackerBounds::isInside(innerTSOS.globalPosition())) {
0692       LogInfo(metname) << "Track in the Tracker: taking the innermost state instead of the state at PCA";
0693       ftsAtVtx = *innerTSOS.freeState();
0694     } else {
0695       if (theAllowNoVtxFlag) {
0696         LogInfo(metname) << "Propagation to PCA failed, taking the innermost state instead of the state at PCA";
0697         ftsAtVtx = *innerTSOS.freeState();
0698       } else {
0699         LogInfo(metname) << "Stand Alone track: this track will be rejected";
0700         return pair<bool, reco::Track>(false, reco::Track());
0701       }
0702     }
0703   }
0704
0705   LogTrace(metname) << "TSOS after the extrapolation at vtx";
0706   LogTrace(metname) << debug.dumpFTS(ftsAtVtx);
0707
0708   GlobalPoint pca = ftsAtVtx.position();
0709   math::XYZPoint persistentPCA(pca.x(), pca.y(), pca.z());
0710   GlobalVector p = ftsAtVtx.momentum();
0711   math::XYZVector persistentMomentum(p.x(), p.y(), p.z());
0712
0713   bool bon = true;
0714   if (fabs(theService->magneticField()->inTesla(GlobalPoint(0, 0, 0)).z()) < 0.01)
0715     bon = false;
0716   double ndof = trajectory.ndof(bon);
0717
0718   reco::Track track(
0719       trajectory.chiSquared(), ndof, persistentPCA, persistentMomentum, ftsAtVtx.charge(), ftsAtVtx.curvilinearError());
0720
0721   return pair<bool, reco::Track>(true, track);
0722 }
0723
0724 pair<bool, reco::Track> MuonTrackLoader::buildTrackUpdatedAtPCA(const reco::Track& track,
0725                                                                 const reco::BeamSpot& beamSpot) const {
0726   const string metname = "Muon|RecoMuon|MuonTrackLoader";
0727   MuonPatternRecoDumper debug;
0728
0729   // build the transient track
0730   reco::TransientTrack transientTrack(track, &*theService->magneticField(), theService->trackingGeometry());
0731
0732   LogTrace(metname) << "Apply the vertex constraint";
0733   pair<bool, FreeTrajectoryState> updateResult = theUpdatorAtVtx->update(transientTrack, beamSpot);
0734
0735   if (!updateResult.first) {
0736     return pair<bool, reco::Track>(false, reco::Track());
0737   }
0738
0739   LogTrace(metname) << "FTS after the vertex constraint";
0740   FreeTrajectoryState& ftsAtVtx = updateResult.second;
0741
0742   LogTrace(metname) << debug.dumpFTS(ftsAtVtx);
0743
0744   GlobalPoint pca = ftsAtVtx.position();
0745   math::XYZPoint persistentPCA(pca.x(), pca.y(), pca.z());
0746   GlobalVector p = ftsAtVtx.momentum();
0747   math::XYZVector persistentMomentum(p.x(), p.y(), p.z());
0748
0749   reco::Track updatedTrack(
0750       track.chi2(), track.ndof(), persistentPCA, persistentMomentum, ftsAtVtx.charge(), ftsAtVtx.curvilinearError());
0751
0752   return pair<bool, reco::Track>(true, updatedTrack);
0753 }
0754
0755 reco::TrackExtra MuonTrackLoader::buildTrackExtra(const Trajectory& trajectory) const {
0756   const string metname = "Muon|RecoMuon|MuonTrackLoader";
0757
0758   const Trajectory::RecHitContainer transRecHits = trajectory.recHits();
0759
0760   // put the collection of TrackingRecHit in the event
0761
0762   // sets the outermost and innermost TSOSs
0763   // FIXME: check it!
0764   TrajectoryStateOnSurface outerTSOS;
0765   TrajectoryStateOnSurface innerTSOS;
0766   unsigned int innerId = 0, outerId = 0;
0767   TrajectoryMeasurement::ConstRecHitPointer outerRecHit;
0768   DetId outerDetId;
0769
0770   if (trajectory.direction() == alongMomentum) {
0771     LogTrace(metname) << "alongMomentum";
0772     outerTSOS = trajectory.lastMeasurement().updatedState();
0773     innerTSOS = trajectory.firstMeasurement().updatedState();
0774     outerId = trajectory.lastMeasurement().recHit()->geographicalId().rawId();
0775     innerId = trajectory.firstMeasurement().recHit()->geographicalId().rawId();
0776     outerRecHit = trajectory.lastMeasurement().recHit();
0777     outerDetId = trajectory.lastMeasurement().recHit()->geographicalId();
0778   } else if (trajectory.direction() == oppositeToMomentum) {
0779     LogTrace(metname) << "oppositeToMomentum";
0780     outerTSOS = trajectory.firstMeasurement().updatedState();
0781     innerTSOS = trajectory.lastMeasurement().updatedState();
0782     outerId = trajectory.firstMeasurement().recHit()->geographicalId().rawId();
0783     innerId = trajectory.lastMeasurement().recHit()->geographicalId().rawId();
0784     outerRecHit = trajectory.firstMeasurement().recHit();
0785     outerDetId = trajectory.firstMeasurement().recHit()->geographicalId();
0786   } else
0787     LogError(metname) << "Wrong propagation direction!";
0788
0789   const GeomDet* outerDet = theService->trackingGeometry()->idToDet(outerDetId);
0790   GlobalPoint outerTSOSPos = outerTSOS.globalParameters().position();
0791   bool inside = outerDet->surface().bounds().inside(outerDet->toLocal(outerTSOSPos));
0792
0793   GlobalPoint hitPos =
0794       (outerRecHit->isValid()) ? outerRecHit->globalPosition() : outerTSOS.globalParameters().position();
0795
0796   if (!inside) {
0797     LogTrace(metname) << "The Global Muon outerMostMeasurementState is not compatible with the recHit detector! "
0798                          "Setting outerMost postition to recHit position if recHit isValid: "
0799                       << outerRecHit->isValid();
0800     LogTrace(metname) << "From " << outerTSOSPos << " to " << hitPos;
0801   }
0802
0803   //build the TrackExtra
0804   GlobalPoint v = (inside) ? outerTSOSPos : hitPos;
0805   GlobalVector p = outerTSOS.globalParameters().momentum();
0806   math::XYZPoint outpos(v.x(), v.y(), v.z());
0807   math::XYZVector outmom(p.x(), p.y(), p.z());
0808
0809   v = innerTSOS.globalParameters().position();
0810   p = innerTSOS.globalParameters().momentum();
0811   math::XYZPoint inpos(v.x(), v.y(), v.z());
0812   math::XYZVector inmom(p.x(), p.y(), p.z());
0813
0814   reco::TrackExtra trackExtra(outpos,
0815                               outmom,
0816                               true,
0817                               inpos,
0818                               inmom,
0819                               true,
0820                               outerTSOS.curvilinearError(),
0821                               outerId,
0822                               innerTSOS.curvilinearError(),
0823                               innerId,
0824                               trajectory.direction(),
0825                               trajectory.seedRef());
0826
0827   return trackExtra;
0828 }