MkFit/plugins/MkFitSeedConverter.cc

0001 #include "FWCore/Framework/interface/global/EDProducer.h"
0002
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006
0007 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0008 #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h"
0010 #include "DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h"
0011 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0012
0013 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0014 #include "DataFormats/TrackerCommon/interface/TrackerDetSide.h"
0015 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0016
0017 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0018 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0019 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0020
0021 #include "MagneticField/Engine/interface/MagneticField.h"
0022 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0023
0024 #include "RecoTracker/MkFit/interface/MkFitGeometry.h"
0025 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0026
0027 #include "RecoTracker/MkFit/interface/MkFitSeedWrapper.h"
0028
0029 // ROOT
0030 #include "Math/SVector.h"
0031 #include "Math/SMatrix.h"
0032
0033 // mkFit includes
0034 #include "RecoTracker/MkFitCMS/interface/LayerNumberConverter.h"
0035 #include "RecoTracker/MkFitCore/interface/Track.h"
0036
0037 class MkFitSeedConverter : public edm::global::EDProducer<> {
0038 public:
0039   explicit MkFitSeedConverter(edm::ParameterSet const& iConfig);
0040   ~MkFitSeedConverter() override = default;
0041
0042   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0043
0044 private:
0045   void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
0046
0047   mkfit::TrackVec convertSeeds(const edm::View<TrajectorySeed>& seeds,
0048                                const TrackerTopology& ttopo,
0049                                const TransientTrackingRecHitBuilder& ttrhBuilder,
0050                                const MagneticField& mf,
0051                                const MkFitGeometry& mkFitGeom) const;
0052
0053   using SVector3 = ROOT::Math::SVector<float, 3>;
0054   using SMatrixSym33 = ROOT::Math::SMatrix<float, 3, 3, ROOT::Math::MatRepSym<float, 3>>;
0055   using SMatrixSym66 = ROOT::Math::SMatrix<float, 6, 6, ROOT::Math::MatRepSym<float, 6>>;
0056
0057   const edm::EDGetTokenT<edm::View<TrajectorySeed>> seedToken_;
0058   const edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> ttrhBuilderToken_;
0059   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> ttopoToken_;
0060   const edm::ESGetToken<MkFitGeometry, TrackerRecoGeometryRecord> mkFitGeomToken_;
0061   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> mfToken_;
0062   const edm::EDPutTokenT<MkFitSeedWrapper> putToken_;
0063   const unsigned int maxNSeeds_;
0064 };
0065
0066 MkFitSeedConverter::MkFitSeedConverter(edm::ParameterSet const& iConfig)
0067     : seedToken_{consumes<edm::View<TrajectorySeed>>(iConfig.getParameter<edm::InputTag>("seeds"))},
0068       ttrhBuilderToken_{esConsumes<TransientTrackingRecHitBuilder, TransientRecHitRecord>(
0069           iConfig.getParameter<edm::ESInputTag>("ttrhBuilder"))},
0070       ttopoToken_{esConsumes<TrackerTopology, TrackerTopologyRcd>()},
0071       mkFitGeomToken_{esConsumes<MkFitGeometry, TrackerRecoGeometryRecord>()},
0072       mfToken_{esConsumes<MagneticField, IdealMagneticFieldRecord>()},
0073       putToken_{produces<MkFitSeedWrapper>()},
0074       maxNSeeds_{iConfig.getParameter<unsigned int>("maxNSeeds")} {}
0075
0076 void MkFitSeedConverter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0077   edm::ParameterSetDescription desc;
0078
0079   desc.add("seeds", edm::InputTag{"initialStepSeeds"});
0080   desc.add("ttrhBuilder", edm::ESInputTag{"", "WithTrackAngle"});
0081   desc.add("maxNSeeds", 500000U);
0082
0083   descriptions.addWithDefaultLabel(desc);
0084 }
0085
0086 void MkFitSeedConverter::produce(edm::StreamID iID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0087   iEvent.emplace(putToken_,
0088                  convertSeeds(iEvent.get(seedToken_),
0089                               iSetup.getData(ttopoToken_),
0090                               iSetup.getData(ttrhBuilderToken_),
0091                               iSetup.getData(mfToken_),
0092                               iSetup.getData(mkFitGeomToken_)));
0093 }
0094
0095 mkfit::TrackVec MkFitSeedConverter::convertSeeds(const edm::View<TrajectorySeed>& seeds,
0096                                                  const TrackerTopology& ttopo,
0097                                                  const TransientTrackingRecHitBuilder& ttrhBuilder,
0098                                                  const MagneticField& mf,
0099                                                  const MkFitGeometry& mkFitGeom) const {
0100   mkfit::TrackVec ret;
0101   if (seeds.size() > maxNSeeds_) {
0102     edm::LogError("TooManySeeds") << "Exceeded maximum number of seeds! maxNSeeds=" << maxNSeeds_
0103                                   << " nSeed=" << seeds.size();
0104     return ret;
0105   }
0106   ret.reserve(seeds.size());
0107
0108   auto isPlusSide = [&ttopo](const DetId& detid) {
0109     return ttopo.side(detid) == static_cast<unsigned>(TrackerDetSide::PosEndcap);
0110   };
0111
0112   int seed_index = 0;
0113   for (const auto& seed : seeds) {
0114     auto const& hitRange = seed.recHits();
0115     const auto lastRecHit = ttrhBuilder.build(&*(hitRange.end() - 1));
0116     const auto tsos = trajectoryStateTransform::transientState(seed.startingState(), lastRecHit->surface(), &mf);
0117     const auto& stateGlobal = tsos.globalParameters();
0118     const auto& gpos = stateGlobal.position();
0119     const auto& gmom = stateGlobal.momentum();
0120     SVector3 pos(gpos.x(), gpos.y(), gpos.z());
0121     SVector3 mom(gmom.x(), gmom.y(), gmom.z());
0122
0123     const auto& cov = tsos.curvilinearError().matrix();
0124     SMatrixSym66 err;  //fill a sub-matrix, mkfit::TrackState will convert internally
0125     for (int i = 0; i < 5; ++i) {
0126       for (int j = i; j < 5; ++j) {
0127         err.At(i, j) = cov[i][j];
0128       }
0129     }
0130
0131     mkfit::TrackState state(tsos.charge(), pos, mom, err);
0132     state.convertFromGlbCurvilinearToCCS();
0133     ret.emplace_back(state, 0, seed_index, 0, nullptr);
0134     LogTrace("MkFitSeedConverter") << "Inserted seed with index " << seed_index;
0135
0136     // Add hits
0137     for (auto const& recHit : hitRange) {
0138       if (not trackerHitRTTI::isFromDet(recHit)) {
0139         throw cms::Exception("Assert") << "Encountered a seed with a hit which is not trackerHitRTTI::isFromDet()";
0140       }
0141       auto& baseTrkRecHit = static_cast<const BaseTrackerRecHit&>(recHit);
0142       if (!baseTrkRecHit.isMatched()) {
0143         const auto& clusterRef = baseTrkRecHit.firstClusterRef();
0144         const auto detId = recHit.geographicalId();
0145         const auto ilay = mkFitGeom.layerNumberConverter().convertLayerNumber(
0146             detId.subdetId(), ttopo.layer(detId), false, ttopo.isStereo(detId), isPlusSide(detId));
0147         LogTrace("MkFitSeedConverter") << " adding hit detid " << detId.rawId() << " index " << clusterRef.index()
0148                                        << " ilay " << ilay;
0149         ret.back().addHitIdx(clusterRef.index(), ilay, 0);  // per-hit chi2 is not known
0150       } else {
0151         auto& matched2D = dynamic_cast<const SiStripMatchedRecHit2D&>(recHit);
0152         const OmniClusterRef* const clRefs[2] = {&matched2D.monoClusterRef(), &matched2D.stereoClusterRef()};
0153         const DetId detIds[2] = {matched2D.monoId(), matched2D.stereoId()};
0154         for (int ii = 0; ii < 2; ++ii) {
0155           const auto& detId = detIds[ii];
0156           const auto ilay = mkFitGeom.layerNumberConverter().convertLayerNumber(
0157               detId.subdetId(), ttopo.layer(detId), false, ttopo.isStereo(detId), isPlusSide(detId));
0158           LogTrace("MkFitSeedConverter") << " adding matched hit detid " << detId.rawId() << " index "
0159                                          << clRefs[ii]->index() << " ilay " << ilay;
0160           ret.back().addHitIdx(clRefs[ii]->index(), ilay, 0);  // per-hit chi2 is not known
0161         }
0162       }
0163     }
0164     ++seed_index;
0165   }
0166   return ret;
0167 }
0168
0169 DEFINE_FWK_MODULE(MkFitSeedConverter);