PixelVertexFinding/plugins/PixelVertexProducerFromSoA.cc

0001 #include "CUDADataFormats/Vertex/interface/ZVertexSoAHeterogeneousHost.h"
0002 #include "CUDADataFormats/Vertex/interface/ZVertexSoAHeterogeneousDevice.h"
0003 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0004 #include "DataFormats/Common/interface/OrphanHandle.h"
0005 #include "DataFormats/TrackReco/interface/Track.h"
0006 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0007 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0008 #include "DataFormats/VertexReco/interface/Vertex.h"
0009 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0010 #include "FWCore/Framework/interface/ConsumesCollector.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/Framework/interface/global/EDProducer.h"
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0019 #include "FWCore/PluginManager/interface/ModuleDef.h"
0020 #include "FWCore/Utilities/interface/EDGetToken.h"
0021 #include "FWCore/Utilities/interface/InputTag.h"
0022 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0023 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0024
0025 #undef PIXVERTEX_DEBUG_PRODUCE
0026
0027 class PixelVertexProducerFromSoA : public edm::global::EDProducer<> {
0028 public:
0029   using IndToEdm = std::vector<uint32_t>;
0030
0031   explicit PixelVertexProducerFromSoA(const edm::ParameterSet &iConfig);
0032   ~PixelVertexProducerFromSoA() override = default;
0033
0034   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0035
0036 private:
0037   void produce(edm::StreamID streamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const override;
0038
0039   edm::EDGetTokenT<ZVertexSoAHost> tokenVertex_;
0040   edm::EDGetTokenT<reco::BeamSpot> tokenBeamSpot_;
0041   edm::EDGetTokenT<reco::TrackCollection> tokenTracks_;
0042   edm::EDGetTokenT<IndToEdm> tokenIndToEdm_;
0043 };
0044
0045 PixelVertexProducerFromSoA::PixelVertexProducerFromSoA(const edm::ParameterSet &conf)
0046     : tokenVertex_(consumes(conf.getParameter<edm::InputTag>("src"))),
0047       tokenBeamSpot_(consumes(conf.getParameter<edm::InputTag>("beamSpot"))),
0048       tokenTracks_(consumes(conf.getParameter<edm::InputTag>("TrackCollection"))),
0049       tokenIndToEdm_(consumes(conf.getParameter<edm::InputTag>("TrackCollection"))) {
0050   produces<reco::VertexCollection>();
0051 }
0052
0053 void PixelVertexProducerFromSoA::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0054   edm::ParameterSetDescription desc;
0055
0056   desc.add<edm::InputTag>("TrackCollection", edm::InputTag("pixelTracks"));
0057   desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0058   desc.add<edm::InputTag>("src", edm::InputTag("pixelVerticesSoA"));
0059
0060   descriptions.add("pixelVertexFromSoA", desc);
0061 }
0062
0063 void PixelVertexProducerFromSoA::produce(edm::StreamID streamID, edm::Event &iEvent, const edm::EventSetup &) const {
0064   auto vertexes = std::make_unique<reco::VertexCollection>();
0065
0066   auto tracksHandle = iEvent.getHandle(tokenTracks_);
0067   auto tracksSize = tracksHandle->size();
0068   auto const &indToEdm = iEvent.get(tokenIndToEdm_);
0069   auto bsHandle = iEvent.getHandle(tokenBeamSpot_);
0070
0071   float x0 = 0, y0 = 0, z0 = 0, dxdz = 0, dydz = 0;
0072   std::vector<int32_t> itrk;
0073   itrk.reserve(64);  // avoid first relocations
0074   if (!bsHandle.isValid()) {
0075     edm::LogWarning("PixelVertexProducer") << "No beamspot found. returning vertexes with (0,0,Z) ";
0076   } else {
0077     const reco::BeamSpot &bs = *bsHandle;
0078     x0 = bs.x0();
0079     y0 = bs.y0();
0080     z0 = bs.z0();
0081     dxdz = bs.dxdz();
0082     dydz = bs.dydz();
0083   }
0084
0085   auto const &soa = iEvent.get(tokenVertex_);
0086
0087   int nv = soa.view().nvFinal();
0088
0089 #ifdef PIXVERTEX_DEBUG_PRODUCE
0090   std::cout << "converting " << nv << " vertices "
0091             << " from " << indToEdm.size() << " tracks" << std::endl;
0092 #endif  // PIXVERTEX_DEBUG_PRODUCE
0093
0094   std::set<uint32_t> uind;  // for verifing index consistency
0095   for (int j = nv - 1; j >= 0; --j) {
0096     auto i = soa.view()[j].sortInd();  // on gpu sorted in ascending order....
0097     assert(i < nv);
0098     uind.insert(i);
0099     assert(itrk.empty());
0100     auto z = soa.view()[i].zv();
0101     auto x = x0 + dxdz * z;
0102     auto y = y0 + dydz * z;
0103     z += z0;
0104     reco::Vertex::Error err;
0105     err(2, 2) = 1.f / soa.view()[i].wv();
0106     err(2, 2) *= 2.;  // artifically inflate error
0107     //Copy also the tracks (no intention to be efficient....)
0108     for (auto k = 0U; k < indToEdm.size(); ++k) {
0109       if (soa.view()[k].idv() == int16_t(i))
0110         itrk.push_back(k);
0111     }
0112     auto nt = itrk.size();
0113     if (nt == 0) {
0114 #ifdef PIXVERTEX_DEBUG_PRODUCE
0115       std::cout << "vertex " << i << " with no tracks..." << std::endl;
0116 #endif  // PIXVERTEX_DEBUG_PRODUCE
0117       continue;
0118     }
0119     if (nt < 2) {
0120       itrk.clear();
0121       continue;
0122     }  // remove outliers
0123     (*vertexes).emplace_back(reco::Vertex::Point(x, y, z), err, soa.view()[i].chi2(), soa.view()[i].ndof(), nt);
0124     auto &v = (*vertexes).back();
0125     v.reserve(itrk.size());
0126     for (auto it : itrk) {
0127       assert(it < int(indToEdm.size()));
0128       auto k = indToEdm[it];
0129       if (k > tracksSize) {
0130         edm::LogWarning("PixelVertexProducer") << "oops track " << it << " does not exists on CPU " << k;
0131         continue;
0132       }
0133       auto tk = reco::TrackRef(tracksHandle, k);
0134       v.add(tk);
0135     }
0136     itrk.clear();
0137   }
0138
0139   LogDebug("PixelVertexProducer") << ": Found " << vertexes->size() << " vertexes\n";
0140   for (unsigned int i = 0; i < vertexes->size(); ++i) {
0141     LogDebug("PixelVertexProducer") << "Vertex number " << i << " has " << (*vertexes)[i].tracksSize()
0142                                     << " tracks with a position of " << (*vertexes)[i].z() << " +- "
0143                                     << std::sqrt((*vertexes)[i].covariance(2, 2));
0144   }
0145
0146   // legacy logic....
0147   if (vertexes->empty() && bsHandle.isValid()) {
0148     const reco::BeamSpot &bs = *bsHandle;
0149
0150     GlobalError bse(bs.rotatedCovariance3D());
0151     if ((bse.cxx() <= 0.) || (bse.cyy() <= 0.) || (bse.czz() <= 0.)) {
0152       AlgebraicSymMatrix33 we;
0153       we(0, 0) = 10000;
0154       we(1, 1) = 10000;
0155       we(2, 2) = 10000;
0156       vertexes->push_back(reco::Vertex(bs.position(), we, 0., 0., 0));
0157
0158       edm::LogInfo("PixelVertexProducer") << "No vertices found. Beamspot with invalid errors " << bse.matrix()
0159                                           << "\nWill put Vertex derived from dummy-fake BeamSpot into Event.\n"
0160                                           << (*vertexes)[0].x() << "\n"
0161                                           << (*vertexes)[0].y() << "\n"
0162                                           << (*vertexes)[0].z() << "\n";
0163     } else {
0164       vertexes->push_back(reco::Vertex(bs.position(), bs.rotatedCovariance3D(), 0., 0., 0));
0165
0166       edm::LogInfo("PixelVertexProducer") << "No vertices found. Will put Vertex derived from BeamSpot into Event:\n"
0167                                           << (*vertexes)[0].x() << "\n"
0168                                           << (*vertexes)[0].y() << "\n"
0169                                           << (*vertexes)[0].z() << "\n";
0170     }
0171   } else if (vertexes->empty() && !bsHandle.isValid()) {
0172     edm::LogWarning("PixelVertexProducer") << "No beamspot and no vertex found. No vertex returned.";
0173   }
0174
0175   iEvent.put(std::move(vertexes));
0176 }
0177
0178 DEFINE_FWK_MODULE(PixelVertexProducerFromSoA);