EgammaPhotonProducers/src/TrackProducerWithSCAssociation.cc

0001 /** \class  TrackProducerWithSCAssociation
0002  **
0003  **
0004  **  \author Nancy Marinelli, U. of Notre Dame, US
0005  **   Modified version of TrackProducer by Giuseppe Cerati
0006  **   to have super cluster - conversion track association
0007  **
0008  ***/
0009
0010 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0011 #include "DataFormats/EgammaTrackReco/interface/TrackCaloClusterAssociation.h"
0012 #include "DataFormats/EgammaTrackReco/interface/TrackCandidateCaloClusterAssociation.h"
0013 #include "DataFormats/TrackReco/interface/Track.h"
0014 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0015 #include "FWCore/Framework/interface/stream/EDProducer.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0018 #include "RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.h"
0019 #include "RecoTracker/TrackProducer/interface/TrackProducerBase.h"
0020 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0021 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0022 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0023 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0025 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0026 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0027
0028 class TrackProducerWithSCAssociation : public TrackProducerBase<reco::Track>, public edm::stream::EDProducer<> {
0029 public:
0030   explicit TrackProducerWithSCAssociation(const edm::ParameterSet& iConfig);
0031
0032   void produce(edm::Event&, const edm::EventSetup&) override;
0033
0034   std::vector<reco::TransientTrack> getTransient(edm::Event&, const edm::EventSetup&);
0035
0036 private:
0037   std::string myname_;
0038   TrackProducerAlgorithm<reco::Track> theAlgo;
0039   std::string conversionTrackCandidateProducer_;
0040   std::string trackCSuperClusterAssociationCollection_;
0041   std::string trackSuperClusterAssociationCollection_;
0042   edm::EDGetTokenT<reco::TrackCandidateCaloClusterPtrAssociation> assoc_token;
0043   edm::OrphanHandle<reco::TrackCollection> rTracks_;
0044   edm::EDGetTokenT<MeasurementTrackerEvent> measurementTrkToken_;
0045   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> ttopoToken_;
0046   bool myTrajectoryInEvent_;
0047   bool validTrackCandidateSCAssociationInput_;
0048
0049   //Same recipe as Ursula's for electrons. Copy this from TrackProducerBase to get the OrphanHandle
0050   //ugly temporary solution!! I agree !
0051   void putInEvt(edm::Event& evt,
0052                 const Propagator* thePropagator,
0053                 const MeasurementTracker* theMeasTk,
0054                 std::unique_ptr<TrackingRecHitCollection> selHits,
0055                 std::unique_ptr<reco::TrackCollection> selTracks,
0056                 std::unique_ptr<reco::TrackExtraCollection> selTrackExtras,
0057                 std::unique_ptr<std::vector<Trajectory>> selTrajectories,
0058                 AlgoProductCollection& algoResults,
0059                 TransientTrackingRecHitBuilder const* hitBuilder,
0060                 const TrackerTopology* ttopo);
0061 };
0062
0063 #include "FWCore/Framework/interface/MakerMacros.h"
0064 DEFINE_FWK_MODULE(TrackProducerWithSCAssociation);
0065
0066 TrackProducerWithSCAssociation::TrackProducerWithSCAssociation(const edm::ParameterSet& iConfig)
0067     : TrackProducerBase<reco::Track>(iConfig.getParameter<bool>("TrajectoryInEvent")),
0068       theAlgo(iConfig),
0069       ttopoToken_(esConsumes()) {
0070   initTrackProducerBase(
0071       iConfig, consumesCollector(), consumes<TrackCandidateCollection>(iConfig.getParameter<edm::InputTag>("src")));
0072   setAlias(iConfig.getParameter<std::string>("@module_label"));
0073
0074   myname_ = iConfig.getParameter<std::string>("ComponentName");
0075   conversionTrackCandidateProducer_ = iConfig.getParameter<std::string>("producer");
0076   trackCSuperClusterAssociationCollection_ = iConfig.getParameter<std::string>("trackCandidateSCAssociationCollection");
0077   trackSuperClusterAssociationCollection_ = iConfig.getParameter<std::string>("recoTrackSCAssociationCollection");
0078   myTrajectoryInEvent_ = iConfig.getParameter<bool>("TrajectoryInEvent");
0079
0080   assoc_token = consumes<reco::TrackCandidateCaloClusterPtrAssociation>(
0081       edm::InputTag(conversionTrackCandidateProducer_, trackCSuperClusterAssociationCollection_));
0082   measurementTrkToken_ = consumes<MeasurementTrackerEvent>(
0083       edm::InputTag("MeasurementTrackerEvent"));  //hardcoded because the original was and no time to fix (sigh)
0084
0085   //register your products
0086   produces<reco::TrackCollection>().setBranchAlias(alias_ + "Tracks");
0087   produces<reco::TrackExtraCollection>().setBranchAlias(alias_ + "TrackExtras");
0088   produces<TrackingRecHitCollection>().setBranchAlias(alias_ + "RecHits");
0089   produces<std::vector<Trajectory>>();
0090   produces<TrajTrackAssociationCollection>();
0091   //  produces< reco::TrackSuperClusterAssociationCollection > (trackSuperClusterAssociationCollection_ );
0092   produces<reco::TrackCaloClusterPtrAssociation>(trackSuperClusterAssociationCollection_);
0093 }
0094
0095 void TrackProducerWithSCAssociation::produce(edm::Event& theEvent, const edm::EventSetup& setup) {
0096   //edm::LogInfo("TrackProducerWithSCAssociation") << "Analyzing event number: " << theEvent.id() << "\n";
0097
0098   //LogDebug("TrackProducerWithSCAssociation") << "Analyzing event number: " << theEvent.id() << "\n";
0099   //  std::cout << " TrackProducerWithSCAssociation Analyzing event number: " << theEvent.id() << "\n";
0100
0101   //
0102   // create empty output collections
0103   //
0104   auto outputRHColl = std::make_unique<TrackingRecHitCollection>();
0105   auto outputTColl = std::make_unique<reco::TrackCollection>();
0106   auto outputTEColl = std::make_unique<reco::TrackExtraCollection>();
0107   auto outputTrajectoryColl = std::make_unique<std::vector<Trajectory>>();
0108   //   Reco Track - Super Cluster Association
0109   auto scTrkAssoc_p = std::make_unique<reco::TrackCaloClusterPtrAssociation>();
0110
0111   //
0112   //declare and get stuff to be retrieved from ES
0113   //
0114   edm::ESHandle<TrackerGeometry> theG;
0115   edm::ESHandle<MagneticField> theMF;
0116   edm::ESHandle<TrajectoryFitter> theFitter;
0117   edm::ESHandle<Propagator> thePropagator;
0118   edm::ESHandle<TransientTrackingRecHitBuilder> theBuilder;
0119   edm::ESHandle<MeasurementTracker> theMeasTk;
0120   getFromES(setup, theG, theMF, theFitter, thePropagator, theMeasTk, theBuilder);
0121
0122   const TrackerTopology* ttopo = &setup.getData(ttopoToken_);
0123
0124   //
0125   //declare and get TrackColection to be retrieved from the event
0126   edm::Handle<TrackCandidateCollection> theTCCollection;
0127   //// Get the association map between candidate out in tracks and the SC where they originated
0128   validTrackCandidateSCAssociationInput_ = true;
0129   edm::Handle<reco::TrackCandidateCaloClusterPtrAssociation> trkCandidateSCAssocHandle;
0130   theEvent.getByToken(assoc_token, trkCandidateSCAssocHandle);
0131   if (!trkCandidateSCAssocHandle.isValid()) {
0132     //    std::cout << "Error! Can't get the product  "<<trackCSuperClusterAssociationCollection_.c_str() << " but keep running. Empty collection will be produced " << "\n";
0133     edm::LogError("TrackProducerWithSCAssociation")
0134         << "Error! Can't get the product  " << trackCSuperClusterAssociationCollection_.c_str()
0135         << " but keep running. Empty collection will be produced "
0136         << "\n";
0137     validTrackCandidateSCAssociationInput_ = false;
0138   }
0139   reco::TrackCandidateCaloClusterPtrAssociation scTrkCandAssCollection = *(trkCandidateSCAssocHandle.product());
0140   if (scTrkCandAssCollection.empty())
0141     validTrackCandidateSCAssociationInput_ = false;
0142
0143   std::vector<int> tccLocations;
0144   AlgoProductCollection algoResults;
0145   reco::BeamSpot bs;
0146
0147   getFromEvt(theEvent, theTCCollection, bs);
0148
0149   if (theTCCollection.failedToGet()) {
0150     edm::LogError("TrackProducerWithSCAssociation")
0151         << "TrackProducerWithSCAssociation could not get the TrackCandidateCollection.";
0152   } else {
0153     //
0154     //run the algorithm
0155     //
0156     //  LogDebug("TrackProducerWithSCAssociation") << "TrackProducerWithSCAssociation run the algorithm" << "\n";
0157     //    theAlgo.runWithCandidate(theG.product(), theMF.product(), *theTCCollection,
0158     //               theFitter.product(), thePropagator.product(), theBuilder.product(), algoResults);
0159     // we have to copy this method from the algo in order to get the association track-seed
0160     // this is ugly temporary code that should be replaced!!!!!
0161     // start of copied code ======================================================
0162
0163     //    std::cout << "TrackProducerWithSCAssociation  Number of TrackCandidates: " << theTCCollection->size() << "\n";
0164     try {
0165       int cont = 0;
0166       int tcc = 0;
0167
0168       for (TrackCandidateCollection::const_iterator i = theTCCollection->begin(); i != theTCCollection->end(); i++) {
0169         const TrackCandidate* theTC = &(*i);
0170         PTrajectoryStateOnDet state = theTC->trajectoryStateOnDet();
0171         const TrajectorySeed& seed = theTC->seed();
0172
0173         //convert PTrajectoryStateOnDet to TrajectoryStateOnSurface
0174
0175         DetId detId(state.detId());
0176         TrajectoryStateOnSurface theTSOS = trajectoryStateTransform::transientState(
0177             state, &(theG.product()->idToDet(detId)->surface()), theMF.product());
0178
0179         //LogDebug("TrackProducerWithSCAssociation")  << "TrackProducerWithSCAssociation  Initial TSOS\n" << theTSOS << "\n";
0180
0181         //convert the TrackingRecHit vector to a TransientTrackingRecHit vector
0182         //meanwhile computes the number of degrees of freedom
0183         TransientTrackingRecHit::RecHitContainer hits;
0184
0185         float ndof = 0;
0186
0187         for (auto const& recHit : theTC->recHits()) {
0188           hits.push_back(theBuilder.product()->build(&recHit));
0189         }
0190
0191         //build Track
0192         // LogDebug("TrackProducerWithSCAssociation") << "TrackProducerWithSCAssociation going to buildTrack"<< "\n";
0193         FitterCloner fc(theFitter.product(), theBuilder.product());
0194         bool ok = theAlgo.buildTrack(
0195             fc.fitter.get(), thePropagator.product(), algoResults, hits, theTSOS, seed, ndof, bs, theTC->seedRef());
0196         // LogDebug("TrackProducerWithSCAssociation")  << "TrackProducerWithSCAssociation buildTrack result: " << ok << "\n";
0197         if (ok) {
0198           cont++;
0199           tccLocations.push_back(tcc);
0200         }
0201         tcc++;
0202       }
0203       edm::LogInfo("TrackProducerWithSCAssociation") << "Number of Tracks found: " << cont << "\n";
0204       //LogDebug("TrackProducerWithSCAssociation") << "TrackProducerWithSCAssociation Number of Tracks found: " << cont << "\n";
0205       // end of copied code ======================================================
0206
0207     } catch (cms::Exception& e) {
0208       edm::LogInfo("TrackProducerWithSCAssociation") << "cms::Exception caught!!!"
0209                                                      << "\n"
0210                                                      << e << "\n";
0211     }
0212     //
0213     //put everything in the event
0214     // we copy putInEvt to get OrphanHandle filled...
0215     putInEvt(theEvent,
0216              thePropagator.product(),
0217              theMeasTk.product(),
0218              std::move(outputRHColl),
0219              std::move(outputTColl),
0220              std::move(outputTEColl),
0221              std::move(outputTrajectoryColl),
0222              algoResults,
0223              theBuilder.product(),
0224              ttopo);
0225
0226     // now construct associationmap and put it in the  event
0227     if (validTrackCandidateSCAssociationInput_) {
0228       int itrack = 0;
0229       std::vector<edm::Ptr<reco::CaloCluster>> caloPtrVec;
0230       for (AlgoProductCollection::iterator i = algoResults.begin(); i != algoResults.end(); i++) {
0231         edm::Ref<TrackCandidateCollection> trackCRef(theTCCollection, tccLocations[itrack]);
0232         const edm::Ptr<reco::CaloCluster>& aClus = (*trkCandidateSCAssocHandle)[trackCRef];
0233         caloPtrVec.push_back(aClus);
0234         itrack++;
0235       }
0236
0237       edm::ValueMap<reco::CaloClusterPtr>::Filler filler(*scTrkAssoc_p);
0238       filler.insert(rTracks_, caloPtrVec.begin(), caloPtrVec.end());
0239       filler.fill();
0240     }
0241
0242     theEvent.put(std::move(scTrkAssoc_p), trackSuperClusterAssociationCollection_);
0243   }
0244 }
0245
0246 std::vector<reco::TransientTrack> TrackProducerWithSCAssociation::getTransient(edm::Event& theEvent,
0247                                                                                const edm::EventSetup& setup) {
0248   edm::LogInfo("TrackProducerWithSCAssociation") << "Analyzing event number: " << theEvent.id() << "\n";
0249   //
0250   // create empty output collections
0251   //
0252   std::vector<reco::TransientTrack> ttks;
0253
0254   //
0255   //declare and get stuff to be retrieved from ES
0256   //
0257   edm::ESHandle<TrackerGeometry> theG;
0258   edm::ESHandle<MagneticField> theMF;
0259   edm::ESHandle<TrajectoryFitter> theFitter;
0260   edm::ESHandle<Propagator> thePropagator;
0261   edm::ESHandle<TransientTrackingRecHitBuilder> theBuilder;
0262   edm::ESHandle<MeasurementTracker> theMeasTk;
0263   getFromES(setup, theG, theMF, theFitter, thePropagator, theMeasTk, theBuilder);
0264
0265   //
0266   //declare and get TrackColection to be retrieved from the event
0267   //
0268   AlgoProductCollection algoResults;
0269   reco::BeamSpot bs;
0270
0271   try {
0272     edm::Handle<TrackCandidateCollection> theTCCollection;
0273     getFromEvt(theEvent, theTCCollection, bs);
0274
0275     //
0276     //run the algorithm
0277     //
0278     //LogDebug("TrackProducerWithSCAssociation") << "TrackProducerWithSCAssociation run the algorithm" << "\n";
0279     theAlgo.runWithCandidate(theG.product(),
0280                              theMF.product(),
0281                              *theTCCollection,
0282                              theFitter.product(),
0283                              thePropagator.product(),
0284                              theBuilder.product(),
0285                              bs,
0286                              algoResults);
0287
0288   } catch (cms::Exception& e) {
0289     edm::LogInfo("TrackProducerWithSCAssociation") << "cms::Exception caught!!!"
0290                                                    << "\n"
0291                                                    << e << "\n";
0292   }
0293
0294   for (auto& prod : algoResults) {
0295     ttks.emplace_back(*prod.track, thePropagator.product()->magneticField());
0296   }
0297
0298   //LogDebug("TrackProducerWithSCAssociation") << "TrackProducerWithSCAssociation end" << "\n";
0299
0300   return ttks;
0301 }
0302
0303 #include "RecoTracker/TransientTrackingRecHit/interface/Traj2TrackHits.h"
0304
0305 void TrackProducerWithSCAssociation::putInEvt(edm::Event& evt,
0306                                               const Propagator* thePropagator,
0307                                               const MeasurementTracker* theMeasTk,
0308                                               std::unique_ptr<TrackingRecHitCollection> selHits,
0309                                               std::unique_ptr<reco::TrackCollection> selTracks,
0310                                               std::unique_ptr<reco::TrackExtraCollection> selTrackExtras,
0311                                               std::unique_ptr<std::vector<Trajectory>> selTrajectories,
0312                                               AlgoProductCollection& algoResults,
0313                                               TransientTrackingRecHitBuilder const* hitBuilder,
0314                                               const TrackerTopology* ttopo) {
0315   TrackingRecHitRefProd rHits = evt.getRefBeforePut<TrackingRecHitCollection>();
0316   reco::TrackExtraRefProd rTrackExtras = evt.getRefBeforePut<reco::TrackExtraCollection>();
0317
0318   edm::Ref<reco::TrackExtraCollection>::key_type idx = 0;
0319   edm::Ref<reco::TrackCollection>::key_type iTkRef = 0;
0320   edm::Ref<std::vector<Trajectory>>::key_type iTjRef = 0;
0321   std::map<unsigned int, unsigned int> tjTkMap;
0322
0323   for (auto& i : algoResults) {
0324     Trajectory* theTraj = i.trajectory;
0325     if (myTrajectoryInEvent_) {
0326       selTrajectories->push_back(*theTraj);
0327       iTjRef++;
0328     }
0329
0330     reco::Track* theTrack = i.track;
0331     PropagationDirection seedDir = i.pDir;
0332
0333     //LogDebug("TrackProducer") << "In KfTrackProducerBase::putInEvt - seedDir=" << seedDir;
0334
0335     reco::Track t = *theTrack;
0336     selTracks->push_back(t);
0337     iTkRef++;
0338
0339     // Store indices in local map (starts at 0)
0340     if (trajectoryInEvent_)
0341       tjTkMap[iTjRef - 1] = iTkRef - 1;
0342
0343     //sets the outermost and innermost TSOSs
0344
0345     TrajectoryStateOnSurface outertsos;
0346     TrajectoryStateOnSurface innertsos;
0347     unsigned int innerId, outerId;
0348
0349     // ---  NOTA BENE: the convention is to sort hits and measurements "along the momentum".
0350     // This is consistent with innermost and outermost labels only for tracks from LHC collision
0351     if (theTraj->direction() == alongMomentum) {
0352       outertsos = theTraj->lastMeasurement().updatedState();
0353       innertsos = theTraj->firstMeasurement().updatedState();
0354       outerId = theTraj->lastMeasurement().recHit()->geographicalId().rawId();
0355       innerId = theTraj->firstMeasurement().recHit()->geographicalId().rawId();
0356     } else {
0357       outertsos = theTraj->firstMeasurement().updatedState();
0358       innertsos = theTraj->lastMeasurement().updatedState();
0359       outerId = theTraj->firstMeasurement().recHit()->geographicalId().rawId();
0360       innerId = theTraj->lastMeasurement().recHit()->geographicalId().rawId();
0361     }
0362     // ---
0363     //build the TrackExtra
0364     GlobalPoint v = outertsos.globalParameters().position();
0365     GlobalVector p = outertsos.globalParameters().momentum();
0366     math::XYZVector outmom(p.x(), p.y(), p.z());
0367     math::XYZPoint outpos(v.x(), v.y(), v.z());
0368     v = innertsos.globalParameters().position();
0369     p = innertsos.globalParameters().momentum();
0370     math::XYZVector inmom(p.x(), p.y(), p.z());
0371     math::XYZPoint inpos(v.x(), v.y(), v.z());
0372
0373     reco::TrackExtraRef teref = reco::TrackExtraRef(rTrackExtras, idx++);
0374     reco::Track& track = selTracks->back();
0375     track.setExtra(teref);
0376
0377     //======= I want to set the second hitPattern here =============
0378     if (theSchool.isValid()) {
0379       edm::Handle<MeasurementTrackerEvent> mte;
0380       evt.getByToken(measurementTrkToken_, mte);
0381       setSecondHitPattern(theTraj, track, thePropagator, &*mte, ttopo);
0382     }
0383     //==============================================================
0384
0385     selTrackExtras->push_back(reco::TrackExtra(outpos,
0386                                                outmom,
0387                                                true,
0388                                                inpos,
0389                                                inmom,
0390                                                true,
0391                                                outertsos.curvilinearError(),
0392                                                outerId,
0393                                                innertsos.curvilinearError(),
0394                                                innerId,
0395                                                seedDir,
0396                                                theTraj->seedRef()));
0397
0398     reco::TrackExtra& tx = selTrackExtras->back();
0399     // ---  NOTA BENE: the convention is to sort hits and measurements "along the momentum".
0400     // This is consistent with innermost and outermost labels only for tracks from LHC collisions
0401     reco::TrackExtra::TrajParams trajParams;
0402     reco::TrackExtra::Chi2sFive chi2s;
0403     Traj2TrackHits t2t;
0404     auto ih = selHits->size();
0405     t2t(*theTraj, *selHits, trajParams, chi2s);
0406     auto ie = selHits->size();
0407     tx.setHits(rHits, ih, ie - ih);
0408     tx.setTrajParams(std::move(trajParams), std::move(chi2s));
0409     for (; ih < ie; ++ih) {
0410       auto const& hit = (*selHits)[ih];
0411       track.appendHitPattern(hit, *ttopo);
0412     }
0413     // ----
0414
0415     delete theTrack;
0416     delete theTraj;
0417   }
0418
0419   //LogTrace("TrackingRegressionTest") << "========== TrackProducer Info ===================";
0420   //LogDebug("TrackProducerWithSCAssociation") << "number of finalTracks: " << selTracks->size() << std::endl;
0421   //for (reco::TrackCollection::const_iterator it = selTracks->begin(); it != selTracks->end(); it++) {
0422   //LogDebug("TrackProducerWithSCAssociation")  << "track's n valid and invalid hit, chi2, pt : "
0423   //                                  << it->found() << " , "
0424   //                                  << it->lost()  <<" , "
0425   //                                  << it->normalizedChi2() << " , "
0426   //           << it->pt() << std::endl;
0427   // }
0428   //LogTrace("TrackingRegressionTest") << "=================================================";
0429
0430   rTracks_ = evt.put(std::move(selTracks));
0431
0432   evt.put(std::move(selTrackExtras));
0433   evt.put(std::move(selHits));
0434
0435   if (myTrajectoryInEvent_) {
0436     edm::OrphanHandle<std::vector<Trajectory>> rTrajs = evt.put(std::move(selTrajectories));
0437
0438     // Now Create traj<->tracks association map
0439     auto trajTrackMap = std::make_unique<TrajTrackAssociationCollection>(rTrajs, rTracks_);
0440     for (std::map<unsigned int, unsigned int>::iterator i = tjTkMap.begin(); i != tjTkMap.end(); i++) {
0441       edm::Ref<std::vector<Trajectory>> trajRef(rTrajs, (*i).first);
0442       edm::Ref<reco::TrackCollection> tkRef(rTracks_, (*i).second);
0443       trajTrackMap->insert(edm::Ref<std::vector<Trajectory>>(rTrajs, (*i).first),
0444                            edm::Ref<reco::TrackCollection>(rTracks_, (*i).second));
0445     }
0446     evt.put(std::move(trajTrackMap));
0447   }
0448 }