Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-27 01:28:01

0001 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0002 #include "DataFormats/Common/interface/OrphanHandle.h"
0003 #include "DataFormats/TrackReco/interface/Track.h"
0004 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0005 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0006 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0007 #include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
0008 #include "DataFormats/GeometrySurface/interface/Plane.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/Framework/interface/global/EDProducer.h"
0014 #include "FWCore/Framework/interface/ConsumesCollector.h"
0015 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "FWCore/Utilities/interface/InputTag.h"
0019 #include "FWCore/PluginManager/interface/ModuleDef.h"
0020 #include "FWCore/Utilities/interface/EDGetToken.h"
0021 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0022 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0023 
0024 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
0025 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0026 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
0027 #include "RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h"
0028 
0029 #include "CUDADataFormats/Common/interface/HostProduct.h"
0030 #include "CUDADataFormats/Track/interface/PixelTrackHeterogeneous.h"
0031 #include "CUDADataFormats/SiPixelCluster/interface/gpuClusteringConstants.h"
0032 
0033 #include "storeTracks.h"
0034 #include "CUDADataFormats/Common/interface/HostProduct.h"
0035 
0036 /**
0037  * This class creates "leagcy"  reco::Track
0038  * objects from the output of SoA CA. 
0039  */
0040 class PixelTrackProducerFromSoA : public edm::global::EDProducer<> {
0041 public:
0042   using IndToEdm = std::vector<uint16_t>;
0043 
0044   explicit PixelTrackProducerFromSoA(const edm::ParameterSet &iConfig);
0045   ~PixelTrackProducerFromSoA() override = default;
0046 
0047   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0048 
0049   //  using HitModuleStart = std::array<uint32_t, gpuClustering::maxNumModules + 1>;
0050   using HMSstorage = HostProduct<uint32_t[]>;
0051 
0052 private:
0053   void produce(edm::StreamID streamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const override;
0054 
0055   // Event Data tokens
0056   const edm::EDGetTokenT<reco::BeamSpot> tBeamSpot_;
0057   const edm::EDGetTokenT<PixelTrackHeterogeneous> tokenTrack_;
0058   const edm::EDGetTokenT<SiPixelRecHitCollectionNew> cpuHits_;
0059   const edm::EDGetTokenT<HMSstorage> hmsToken_;
0060   // Event Setup tokens
0061   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> idealMagneticFieldToken_;
0062   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> ttTopoToken_;
0063 
0064   int32_t const minNumberOfHits_;
0065   pixelTrack::Quality const minQuality_;
0066 };
0067 
0068 PixelTrackProducerFromSoA::PixelTrackProducerFromSoA(const edm::ParameterSet &iConfig)
0069     : tBeamSpot_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))),
0070       tokenTrack_(consumes<PixelTrackHeterogeneous>(iConfig.getParameter<edm::InputTag>("trackSrc"))),
0071       cpuHits_(consumes<SiPixelRecHitCollectionNew>(iConfig.getParameter<edm::InputTag>("pixelRecHitLegacySrc"))),
0072       hmsToken_(consumes<HMSstorage>(iConfig.getParameter<edm::InputTag>("pixelRecHitLegacySrc"))),
0073       idealMagneticFieldToken_(esConsumes()),
0074       ttTopoToken_(esConsumes()),
0075       minNumberOfHits_(iConfig.getParameter<int>("minNumberOfHits")),
0076       minQuality_(pixelTrack::qualityByName(iConfig.getParameter<std::string>("minQuality"))) {
0077   if (minQuality_ == pixelTrack::Quality::notQuality) {
0078     throw cms::Exception("PixelTrackConfiguration")
0079         << iConfig.getParameter<std::string>("minQuality") + " is not a pixelTrack::Quality";
0080   }
0081   if (minQuality_ < pixelTrack::Quality::dup) {
0082     throw cms::Exception("PixelTrackConfiguration")
0083         << iConfig.getParameter<std::string>("minQuality") + " not supported";
0084   }
0085   produces<TrackingRecHitCollection>();
0086   produces<reco::TrackExtraCollection>();
0087   // TrackCollection refers to TrackingRechit and TrackExtra
0088   // collections, need to declare its production after them to work
0089   // around a rare race condition in framework scheduling
0090   produces<reco::TrackCollection>();
0091   produces<IndToEdm>();
0092 }
0093 
0094 void PixelTrackProducerFromSoA::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0095   edm::ParameterSetDescription desc;
0096   desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0097   desc.add<edm::InputTag>("trackSrc", edm::InputTag("pixelTracksSoA"));
0098   desc.add<edm::InputTag>("pixelRecHitLegacySrc", edm::InputTag("siPixelRecHitsPreSplittingLegacy"));
0099   desc.add<int>("minNumberOfHits", 0);
0100   desc.add<std::string>("minQuality", "loose");
0101   descriptions.addWithDefaultLabel(desc);
0102 }
0103 
0104 void PixelTrackProducerFromSoA::produce(edm::StreamID streamID,
0105                                         edm::Event &iEvent,
0106                                         const edm::EventSetup &iSetup) const {
0107   // enum class Quality : uint8_t { bad = 0, edup, dup, loose, strict, tight, highPurity };
0108   reco::TrackBase::TrackQuality recoQuality[] = {reco::TrackBase::undefQuality,
0109                                                  reco::TrackBase::undefQuality,
0110                                                  reco::TrackBase::discarded,
0111                                                  reco::TrackBase::loose,
0112                                                  reco::TrackBase::tight,
0113                                                  reco::TrackBase::tight,
0114                                                  reco::TrackBase::highPurity};
0115   assert(reco::TrackBase::highPurity == recoQuality[int(pixelTrack::Quality::highPurity)]);
0116 
0117   // std::cout << "Converting gpu helix in reco tracks" << std::endl;
0118 
0119   auto indToEdmP = std::make_unique<IndToEdm>();
0120   auto &indToEdm = *indToEdmP;
0121 
0122   auto const &idealField = iSetup.getData(idealMagneticFieldToken_);
0123 
0124   pixeltrackfitting::TracksWithRecHits tracks;
0125 
0126   auto const &httopo = iSetup.getData(ttTopoToken_);
0127 
0128   const auto &bsh = iEvent.get(tBeamSpot_);
0129   GlobalPoint bs(bsh.x0(), bsh.y0(), bsh.z0());
0130 
0131   auto const &rechits = iEvent.get(cpuHits_);
0132   std::vector<TrackingRecHit const *> hitmap;
0133   auto const &rcs = rechits.data();
0134   auto nhits = rcs.size();
0135   hitmap.resize(nhits, nullptr);
0136 
0137   auto const *hitsModuleStart = iEvent.get(hmsToken_).get();
0138   auto fc = hitsModuleStart;
0139 
0140   for (auto const &h : rcs) {
0141     auto const &thit = static_cast<BaseTrackerRecHit const &>(h);
0142     auto detI = thit.det()->index();
0143     auto const &clus = thit.firstClusterRef();
0144     assert(clus.isPixel());
0145     auto i = fc[detI] + clus.pixelCluster().originalId();
0146     if (i >= hitmap.size())
0147       hitmap.resize(i + 256, nullptr);  // only in case of hit overflow in one module
0148     assert(nullptr == hitmap[i]);
0149     hitmap[i] = &h;
0150   }
0151 
0152   std::vector<const TrackingRecHit *> hits;
0153   hits.reserve(5);
0154 
0155   const auto &tsoa = *iEvent.get(tokenTrack_);
0156 
0157   auto const *quality = tsoa.qualityData();
0158   auto const &fit = tsoa.stateAtBS;
0159   auto const &hitIndices = tsoa.hitIndices;
0160   auto nTracks = tsoa.nTracks();
0161 
0162   tracks.reserve(nTracks);
0163 
0164   int32_t nt = 0;
0165 
0166   //sort index by pt
0167   std::vector<int32_t> sortIdxs(nTracks);
0168   std::iota(sortIdxs.begin(), sortIdxs.end(), 0);
0169   std::sort(
0170       sortIdxs.begin(), sortIdxs.end(), [&](int32_t const i1, int32_t const i2) { return tsoa.pt(i1) > tsoa.pt(i2); });
0171 
0172   //store the index of the SoA: indToEdm[index_SoAtrack] -> index_edmTrack (if it exists)
0173   indToEdm.resize(sortIdxs.size(), -1);
0174   for (const auto &it : sortIdxs) {
0175     auto nHits = tsoa.nHits(it);
0176     assert(nHits >= 3);
0177     auto q = quality[it];
0178     if (q < minQuality_)
0179       continue;
0180     if (nHits < minNumberOfHits_)
0181       continue;
0182     indToEdm[it] = nt;
0183     ++nt;
0184 
0185     hits.resize(nHits);
0186     auto b = hitIndices.begin(it);
0187     for (int iHit = 0; iHit < nHits; ++iHit)
0188       hits[iHit] = hitmap[*(b + iHit)];
0189 
0190     // mind: this values are respect the beamspot!
0191 
0192     float chi2 = tsoa.chi2(it);
0193     float phi = tsoa.phi(it);
0194 
0195     riemannFit::Vector5d ipar, opar;
0196     riemannFit::Matrix5d icov, ocov;
0197     fit.copyToDense(ipar, icov, it);
0198     riemannFit::transformToPerigeePlane(ipar, icov, opar, ocov);
0199 
0200     LocalTrajectoryParameters lpar(opar(0), opar(1), opar(2), opar(3), opar(4), 1.);
0201     AlgebraicSymMatrix55 m;
0202     for (int i = 0; i < 5; ++i)
0203       for (int j = i; j < 5; ++j)
0204         m(i, j) = ocov(i, j);
0205 
0206     float sp = std::sin(phi);
0207     float cp = std::cos(phi);
0208     Surface::RotationType rot(sp, -cp, 0, 0, 0, -1.f, cp, sp, 0);
0209 
0210     Plane impPointPlane(bs, rot);
0211     GlobalTrajectoryParameters gp(
0212         impPointPlane.toGlobal(lpar.position()), impPointPlane.toGlobal(lpar.momentum()), lpar.charge(), &idealField);
0213     JacobianLocalToCurvilinear jl2c(impPointPlane, lpar, idealField);
0214 
0215     AlgebraicSymMatrix55 mo = ROOT::Math::Similarity(jl2c.jacobian(), m);
0216 
0217     int ndof = 2 * hits.size() - 5;
0218     chi2 = chi2 * ndof;
0219     GlobalPoint vv = gp.position();
0220     math::XYZPoint pos(vv.x(), vv.y(), vv.z());
0221     GlobalVector pp = gp.momentum();
0222     math::XYZVector mom(pp.x(), pp.y(), pp.z());
0223 
0224     auto track = std::make_unique<reco::Track>(chi2, ndof, pos, mom, gp.charge(), CurvilinearTrajectoryError(mo));
0225 
0226     // bad and edup not supported as fit not present or not reliable
0227     auto tkq = recoQuality[int(q)];
0228     track->setQuality(tkq);
0229     // loose,tight and HP are inclusive
0230     if (reco::TrackBase::highPurity == tkq) {
0231       track->setQuality(reco::TrackBase::tight);
0232       track->setQuality(reco::TrackBase::loose);
0233     } else if (reco::TrackBase::tight == tkq) {
0234       track->setQuality(reco::TrackBase::loose);
0235     }
0236     track->setQuality(tkq);
0237     // filter???
0238     tracks.emplace_back(track.release(), hits);
0239   }
0240   // std::cout << "processed " << nt << " good tuples " << tracks.size() << "out of " << indToEdm.size() << std::endl;
0241 
0242   // store tracks
0243   storeTracks(iEvent, tracks, httopo);
0244   iEvent.put(std::move(indToEdmP));
0245 }
0246 
0247 DEFINE_FWK_MODULE(PixelTrackProducerFromSoA);