Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-04-02 03:11:04

0001 #include "CUDADataFormats/Track/interface/PixelTrackHeterogeneous.h"
0002 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0003 #include "DataFormats/GeometrySurface/interface/Plane.h"
0004 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0005 #include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
0006 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0007 #include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
0008 #include "FWCore/Framework/interface/ConsumesCollector.h"
0009 #include "FWCore/Framework/interface/ESHandle.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/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0017 #include "FWCore/PluginManager/interface/ModuleDef.h"
0018 #include "FWCore/Utilities/interface/EDGetToken.h"
0019 #include "FWCore/Utilities/interface/InputTag.h"
0020 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0021 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0022 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0023 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0024 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0025 #include "RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h"
0026 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
0027 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0028 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0029 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
0030 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0031 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0032 
0033 /*
0034   produces seeds directly from cuda produced tuples
0035 */
0036 class SeedProducerFromSoA : public edm::global::EDProducer<> {
0037 public:
0038   explicit SeedProducerFromSoA(const edm::ParameterSet& iConfig);
0039   ~SeedProducerFromSoA() override = default;
0040 
0041   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0042 
0043 private:
0044   void produce(edm::StreamID streamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
0045 
0046   // Event data tokens
0047   const edm::EDGetTokenT<reco::BeamSpot> tBeamSpot_;
0048   const edm::EDGetTokenT<PixelTrackHeterogeneous> tokenTrack_;
0049   // Event setup tokens
0050   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> idealMagneticFieldToken_;
0051   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerDigiGeometryToken_;
0052   const edm::ESGetToken<Propagator, TrackingComponentsRecord> trackerPropagatorToken_;
0053   int32_t minNumberOfHits_;
0054 };
0055 
0056 SeedProducerFromSoA::SeedProducerFromSoA(const edm::ParameterSet& iConfig)
0057     : tBeamSpot_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))),
0058       tokenTrack_(consumes<PixelTrackHeterogeneous>(iConfig.getParameter<edm::InputTag>("src"))),
0059       idealMagneticFieldToken_(esConsumes()),
0060       trackerDigiGeometryToken_(esConsumes()),
0061       trackerPropagatorToken_(esConsumes(edm::ESInputTag("PropagatorWithMaterial"))),
0062       minNumberOfHits_(iConfig.getParameter<int>("minNumberOfHits"))
0063 
0064 {
0065   produces<TrajectorySeedCollection>();
0066 }
0067 
0068 void SeedProducerFromSoA::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0069   edm::ParameterSetDescription desc;
0070   desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0071   desc.add<edm::InputTag>("src", edm::InputTag("pixelTrackSoA"));
0072   desc.add<int>("minNumberOfHits", 0);
0073 
0074   descriptions.addWithDefaultLabel(desc);
0075 }
0076 
0077 void SeedProducerFromSoA::produce(edm::StreamID streamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0078   // std::cout << "Converting gpu helix to trajectory seed" << std::endl;
0079   auto result = std::make_unique<TrajectorySeedCollection>();
0080 
0081   auto const& fieldESH = iSetup.getHandle(idealMagneticFieldToken_);
0082   auto const& tracker = iSetup.getHandle(trackerDigiGeometryToken_);
0083   auto const& dus = tracker->detUnits();
0084 
0085   auto const& propagatorHandle = iSetup.getHandle(trackerPropagatorToken_);
0086   const Propagator* propagator = &(*propagatorHandle);
0087 
0088   const auto& bsh = iEvent.get(tBeamSpot_);
0089   // std::cout << "beamspot " << bsh.x0() << ' ' << bsh.y0() << ' ' << bsh.z0() << std::endl;
0090   GlobalPoint bs(bsh.x0(), bsh.y0(), bsh.z0());
0091 
0092   const auto& tsoa = *(iEvent.get(tokenTrack_));
0093 
0094   auto const* quality = tsoa.qualityData();
0095   auto const& fit = tsoa.stateAtBS;
0096   auto const& detIndices = tsoa.detIndices;
0097   auto maxTracks = tsoa.stride();
0098 
0099   int32_t nt = 0;
0100   for (int32_t it = 0; it < maxTracks; ++it) {
0101     auto nHits = tsoa.nHits(it);
0102     if (nHits == 0)
0103       break;  // this is a guard: maybe we need to move to nTracks...
0104 
0105     auto q = quality[it];
0106     if (q != pixelTrack::Quality::loose)
0107       continue;  // FIXME
0108     if (nHits < minNumberOfHits_)
0109       continue;
0110     ++nt;
0111 
0112     // fill hits with invalid just to hold the detId
0113     auto b = detIndices.begin(it);
0114     edm::OwnVector<TrackingRecHit> hits;
0115     for (int iHit = 0; iHit < nHits; ++iHit) {
0116       auto const* det = dus[*(b + iHit)];
0117       // FIXME at some point get a proper type ...
0118       hits.push_back(new InvalidTrackingRecHit(*det, TrackingRecHit::bad));
0119     }
0120 
0121     // mind: this values are respect the beamspot!
0122 
0123     float phi = tsoa.phi(it);
0124 
0125     riemannFit::Vector5d ipar, opar;
0126     riemannFit::Matrix5d icov, ocov;
0127     fit.copyToDense(ipar, icov, it);
0128     riemannFit::transformToPerigeePlane(ipar, icov, opar, ocov);
0129 
0130     LocalTrajectoryParameters lpar(opar(0), opar(1), opar(2), opar(3), opar(4), 1.);
0131     AlgebraicSymMatrix55 m;
0132     for (int i = 0; i < 5; ++i)
0133       for (int j = i; j < 5; ++j)
0134         m(i, j) = ocov(i, j);
0135 
0136     float sp = std::sin(phi);
0137     float cp = std::cos(phi);
0138     Surface::RotationType rot(sp, -cp, 0, 0, 0, -1.f, cp, sp, 0);
0139 
0140     Plane impPointPlane(bs, rot);
0141     GlobalTrajectoryParameters gp(impPointPlane.toGlobal(lpar.position()),
0142                                   impPointPlane.toGlobal(lpar.momentum()),
0143                                   lpar.charge(),
0144                                   fieldESH.product());
0145 
0146     JacobianLocalToCurvilinear jl2c(impPointPlane, lpar, *fieldESH.product());
0147 
0148     AlgebraicSymMatrix55 mo = ROOT::Math::Similarity(jl2c.jacobian(), m);
0149 
0150     FreeTrajectoryState fts(gp, CurvilinearTrajectoryError(mo));
0151 
0152     auto const& lastHit = hits.back();
0153 
0154     TrajectoryStateOnSurface outerState = propagator->propagate(fts, *lastHit.surface());
0155 
0156     if (!outerState.isValid()) {
0157       edm::LogError("SeedFromGPU") << " was trying to create a seed from:\n"
0158                                    << fts << "\n propagating to: " << lastHit.geographicalId().rawId();
0159       continue;
0160     }
0161 
0162     auto const& pTraj = trajectoryStateTransform::persistentState(outerState, lastHit.geographicalId().rawId());
0163 
0164     result->emplace_back(pTraj, hits, alongMomentum);
0165   }
0166 
0167   iEvent.put(std::move(result));
0168 }
0169 
0170 DEFINE_FWK_MODULE(SeedProducerFromSoA);