MCTruth/plugins/SeedToTrackProducerBase.cc

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