MCTruth/plugins/SeedToTrackProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    SeedToTrackProducer
0004 // Class:      SeedToTrackProducer
0005 //
0006 /**\class SeedToTrackProducer SeedToTrackProducer.cc
0007  hugues/SeedToTrackProducer/plugins/SeedToTrackProducer.cc
0008
0009  Description: [one line class summary]
0010
0011  Implementation:
0012      [Notes on implementation]
0013 */
0014 //
0015 // Original Author:  Hugues Brun
0016 //         Created:  Tue, 05 Nov 2013 13:42:04 GMT
0017 // $Id$
0018 //
0019 //
0020
0021 #include "SeedToTrackProducer.h"
0022
0023 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0024 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0025
0026 //
0027 // constructors and destructor
0028 //
0029 SeedToTrackProducer::SeedToTrackProducer(const edm::ParameterSet &iConfig)
0030     : theMGFieldToken(esConsumes()), theTrackingGeometryToken(esConsumes()), theTopoToken(esConsumes()) {
0031   L2seedsTagT_ = consumes<TrajectorySeedCollection>(iConfig.getParameter<edm::InputTag>("L2seedsCollection"));
0032   L2seedsTagS_ = consumes<edm::View<TrajectorySeed>>(iConfig.getParameter<edm::InputTag>("L2seedsCollection"));
0033
0034   produces<reco::TrackCollection>();
0035   produces<reco::TrackExtraCollection>();
0036   produces<TrackingRecHitCollection>();
0037 }
0038
0039 //
0040 // member functions
0041 //
0042
0043 // ------------ method called to produce the data  ------------
0044 void SeedToTrackProducer::produce(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0045   using namespace edm;
0046   using namespace std;
0047
0048   std::unique_ptr<reco::TrackCollection> selectedTracks(new reco::TrackCollection);
0049   std::unique_ptr<reco::TrackExtraCollection> selectedTrackExtras(new reco::TrackExtraCollection());
0050   std::unique_ptr<TrackingRecHitCollection> selectedTrackHits(new TrackingRecHitCollection());
0051
0052   reco::TrackRefProd rTracks = iEvent.getRefBeforePut<reco::TrackCollection>();
0053   reco::TrackExtraRefProd rTrackExtras = iEvent.getRefBeforePut<reco::TrackExtraCollection>();
0054   TrackingRecHitRefProd rHits = iEvent.getRefBeforePut<TrackingRecHitCollection>();
0055
0056   edm::Ref<reco::TrackExtraCollection>::key_type hidx = 0;
0057   edm::Ref<reco::TrackExtraCollection>::key_type idx = 0;
0058
0059   // magnetic fied and detector geometry
0060   auto const &mgField = iSetup.getData(theMGFieldToken);
0061   auto const &trackingGeometry = iSetup.getData(theTrackingGeometryToken);
0062
0063   const TrackerTopology &ttopo = iSetup.getData(theTopoToken);
0064
0065   // now read the L2 seeds collection :
0066   edm::Handle<TrajectorySeedCollection> L2seedsCollection;
0067   iEvent.getByToken(L2seedsTagT_, L2seedsCollection);
0068   const std::vector<TrajectorySeed> *L2seeds = nullptr;
0069   if (L2seedsCollection.isValid())
0070     L2seeds = L2seedsCollection.product();
0071   else
0072     edm::LogError("SeedToTrackProducer") << "L2 seeds collection not found !! " << endl;
0073
0074   edm::Handle<edm::View<TrajectorySeed>> seedHandle;
0075   iEvent.getByToken(L2seedsTagS_, seedHandle);
0076
0077   // now  loop on the seeds :
0078   for (unsigned int i = 0; i < L2seeds->size(); i++) {
0079     // get the kinematic extrapolation from the seed
0080     TrajectoryStateOnSurface theTrajectory = seedTransientState(L2seeds->at(i), mgField, trackingGeometry);
0081     float seedEta = theTrajectory.globalMomentum().eta();
0082     float seedPhi = theTrajectory.globalMomentum().phi();
0083     float seedPt = theTrajectory.globalMomentum().perp();
0084     CovarianceMatrix matrixSeedErr = theTrajectory.curvilinearError().matrix();
0085     edm::LogVerbatim("SeedToTrackProducer")
0086         << "seedPt=" << seedPt << " seedEta=" << seedEta << " seedPhi=" << seedPhi << endl;
0087     /*AlgebraicSymMatrix66 errors = theTrajectory.cartesianError().matrix();
0088     double partialPterror =
0089     errors(3,3)*pow(theTrajectory.globalMomentum().x(),2) +
0090     errors(4,4)*pow(theTrajectory.globalMomentum().y(),2);
0091     edm::LogVerbatim("SeedToTrackProducer") <<  "seedPtError=" <<
0092     sqrt(partialPterror)/theTrajectory.globalMomentum().perp() <<
0093     "seedPhiError=" << theTrajectory.curvilinearError().matrix()(2,2) << endl;*/
0094     // fill the track in a way that its pt, phi and eta will be the same as the
0095     // seed
0096     math::XYZPoint initPoint(0, 0, 0);
0097     math::XYZVector initMom(seedPt * cos(seedPhi), seedPt * sin(seedPhi), seedPt * sinh(seedEta));
0098     reco::Track theTrack(1,
0099                          1,  // dummy Chi2 and ndof
0100                          initPoint,
0101                          initMom,
0102                          1,
0103                          matrixSeedErr,
0104                          reco::TrackBase::TrackAlgorithm::globalMuon,
0105                          reco::TrackBase::TrackQuality::tight);
0106
0107     // fill the extra track with dummy information
0108     math::XYZPoint dummyFinalPoint(1, 1, 1);
0109     math::XYZVector dummyFinalMom(0, 0, 10);
0110     edm::RefToBase<TrajectorySeed> seed(seedHandle, i);
0111     CovarianceMatrix matrixExtra = ROOT::Math::SMatrixIdentity();
0112     reco::TrackExtra theTrackExtra(dummyFinalPoint,
0113                                    dummyFinalMom,
0114                                    true,
0115                                    initPoint,
0116                                    initMom,
0117                                    true,
0118                                    matrixSeedErr,
0119                                    1,
0120                                    matrixExtra,
0121                                    2,
0122                                    (L2seeds->at(i)).direction(),
0123                                    seed);
0124     theTrack.setExtra(reco::TrackExtraRef(rTrackExtras, idx++));
0125     edm::LogVerbatim("SeedToTrackProducer")
0126         << "trackPt=" << theTrack.pt() << " trackEta=" << theTrack.eta() << " trackPhi=" << theTrack.phi() << endl;
0127     edm::LogVerbatim("SeedToTrackProducer")
0128         << "trackPtError=" << theTrack.ptError() << "trackPhiError=" << theTrack.phiError() << endl;
0129
0130     // fill the seed segments in the track
0131     unsigned int nHitsAdded = 0;
0132     for (auto const &recHit : L2seeds->at(i).recHits()) {
0133       TrackingRecHit *hit = recHit.clone();
0134       theTrack.appendHitPattern(*hit, ttopo);
0135       selectedTrackHits->push_back(hit);
0136       nHitsAdded++;
0137     }
0138     theTrackExtra.setHits(rHits, hidx, nHitsAdded);
0139     hidx += nHitsAdded;
0140     selectedTracks->push_back(theTrack);
0141     selectedTrackExtras->push_back(theTrackExtra);
0142   }
0143   iEvent.put(std::move(selectedTracks));
0144   iEvent.put(std::move(selectedTrackExtras));
0145   iEvent.put(std::move(selectedTrackHits));
0146 }
0147
0148 TrajectoryStateOnSurface SeedToTrackProducer::seedTransientState(const TrajectorySeed &tmpSeed,
0149                                                                  const MagneticField &mgField,
0150                                                                  const GlobalTrackingGeometry &trackingGeometry) const {
0151   PTrajectoryStateOnDet tmpTSOD = tmpSeed.startingState();
0152   DetId tmpDetId(tmpTSOD.detId());
0153   const GeomDet *tmpGeomDet = trackingGeometry.idToDet(tmpDetId);
0154   TrajectoryStateOnSurface tmpTSOS =
0155       trajectoryStateTransform::transientState(tmpTSOD, &(tmpGeomDet->surface()), &mgField);
0156   return tmpTSOS;
0157 }
0158
0159 // define this as a plug-in
0160 DEFINE_FWK_MODULE(SeedToTrackProducer);