SiPixelHeterogeneous/plugins/SiPixelCompareVertices.cc

0001 // TODO: change file name to SiPixelCompareVerticesSoA.cc when CUDA code is removed
0002
0003 // -*- C++ -*-
0004 // Package:    SiPixelCompareVertices
0005 // Class:      SiPixelCompareVertices
0006 //
0007 /**\class SiPixelCompareVertices SiPixelCompareVertices.cc
0008 */
0009 //
0010 // Author: Suvankar Roy Chowdhury
0011 //
0012
0013 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0014 #include "DQMServices/Core/interface/DQMStore.h"
0015 #include "DQMServices/Core/interface/MonitorElement.h"
0016 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0017 #include "DataFormats/Common/interface/Handle.h"
0018 #include "DataFormats/VertexSoA/interface/ZVertexHost.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/Frameworkfwd.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025
0026 // TODO: change class name to SiPixelCompareVerticesSoA when CUDA code is removed
0027 class SiPixelCompareVertices : public DQMEDAnalyzer {
0028 public:
0029   using IndToEdm = std::vector<uint16_t>;
0030   explicit SiPixelCompareVertices(const edm::ParameterSet&);
0031   ~SiPixelCompareVertices() override = default;
0032   void bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) override;
0033   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0034   // analyzeSeparate is templated to accept distinct types of SoAs
0035   // The default use case is to use vertices from Alpaka reconstructed on CPU and GPU;
0036   template <typename U, typename V>
0037   void analyzeSeparate(U tokenRef, V tokenTar, const edm::Event& iEvent);
0038   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0039
0040 private:
0041   // these two are both on Host but originally they have been produced on Host or on Device
0042   const edm::EDGetTokenT<ZVertexHost> tokenSoAVertexReferenceSoA_;
0043   const edm::EDGetTokenT<ZVertexHost> tokenSoAVertexTargetSoA_;
0044   const edm::EDGetTokenT<reco::BeamSpot> tokenBeamSpot_;
0045   const std::string topFolderName_;
0046   const float dzCut_;
0047   MonitorElement* hnVertex_;
0048   MonitorElement* hx_;
0049   MonitorElement* hy_;
0050   MonitorElement* hz_;
0051   MonitorElement* hchi2_;
0052   MonitorElement* hchi2oNdof_;
0053   MonitorElement* hptv2_;
0054   MonitorElement* hntrks_;
0055   MonitorElement* hxdiff_;
0056   MonitorElement* hydiff_;
0057   MonitorElement* hzdiff_;
0058 };
0059
0060 //
0061 // constructors
0062 //
0063
0064 SiPixelCompareVertices::SiPixelCompareVertices(const edm::ParameterSet& iConfig)
0065     : tokenSoAVertexReferenceSoA_(
0066           consumes<ZVertexHost>(iConfig.getParameter<edm::InputTag>("pixelVertexReferenceSoA"))),
0067       tokenSoAVertexTargetSoA_(consumes<ZVertexHost>(iConfig.getParameter<edm::InputTag>("pixelVertexTargetSoA"))),
0068       tokenBeamSpot_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpotSrc"))),
0069       topFolderName_(iConfig.getParameter<std::string>("topFolderName")),
0070       dzCut_(iConfig.getParameter<double>("dzCut")) {}
0071
0072 template <typename U, typename V>
0073 void SiPixelCompareVertices::analyzeSeparate(U tokenRef, V tokenTar, const edm::Event& iEvent) {
0074   const auto& vsoaHandleRef = iEvent.getHandle(tokenRef);
0075   const auto& vsoaHandleTar = iEvent.getHandle(tokenTar);
0076
0077   if (not vsoaHandleRef or not vsoaHandleTar) {
0078     edm::LogWarning out("SiPixelCompareVertices");
0079     if (not vsoaHandleRef) {
0080       out << "reference vertices not found; ";
0081     }
0082     if (not vsoaHandleTar) {
0083       out << "target vertices not found; ";
0084     }
0085     out << "the comparison will not run.";
0086     return;
0087   }
0088
0089   auto const& vsoaRef = *vsoaHandleRef;
0090   int nVerticesRef = vsoaRef.view().nvFinal();
0091   auto const& vsoaTar = *vsoaHandleTar;
0092   int nVerticesTar = vsoaTar.view().nvFinal();
0093
0094   auto bsHandle = iEvent.getHandle(tokenBeamSpot_);
0095   float x0 = 0., y0 = 0., z0 = 0., dxdz = 0., dydz = 0.;
0096   if (not bsHandle.isValid()) {
0097     edm::LogWarning("SiPixelCompareVertices") << "No beamspot found, returning vertexes with (0,0,Z).";
0098   } else {
0099     const reco::BeamSpot& bs = *bsHandle;
0100     x0 = bs.x0();
0101     y0 = bs.y0();
0102     z0 = bs.z0();
0103     dxdz = bs.dxdz();
0104     dydz = bs.dydz();
0105   }
0106
0107   for (int ivc = 0; ivc < nVerticesRef; ivc++) {
0108     auto sic = vsoaRef.view()[ivc].sortInd();
0109     auto zc = vsoaRef.view()[sic].zv();
0110     auto xc = x0 + dxdz * zc;
0111     auto yc = y0 + dydz * zc;
0112     zc += z0;
0113
0114     auto ndofRef = vsoaRef.template view<reco::ZVertexTracksSoA>()[sic].ndof();
0115     auto chi2Ref = vsoaRef.view()[sic].chi2();
0116
0117     const int32_t notFound = -1;
0118     int32_t closestVtxidx = notFound;
0119     float mindz = dzCut_;
0120
0121     for (int ivg = 0; ivg < nVerticesTar; ivg++) {
0122       auto sig = vsoaTar.view()[ivg].sortInd();
0123       auto zgc = vsoaTar.view()[sig].zv() + z0;
0124       auto zDist = std::abs(zc - zgc);
0125       //insert some matching condition
0126       if (zDist > dzCut_)
0127         continue;
0128       if (mindz > zDist) {
0129         mindz = zDist;
0130         closestVtxidx = sig;
0131       }
0132     }
0133     if (closestVtxidx == notFound)
0134       continue;
0135
0136     auto zg = vsoaTar.view()[closestVtxidx].zv();
0137     auto xg = x0 + dxdz * zg;
0138     auto yg = y0 + dydz * zg;
0139     zg += z0;
0140     auto ndofTar = vsoaTar.template view<reco::ZVertexTracksSoA>()[closestVtxidx].ndof();
0141     auto chi2Tar = vsoaTar.view()[closestVtxidx].chi2();
0142
0143     hx_->Fill(xc - x0, xg - x0);
0144     hy_->Fill(yc - y0, yg - y0);
0145     hz_->Fill(zc, zg);
0146     hxdiff_->Fill(xc - xg);
0147     hydiff_->Fill(yc - yg);
0148     hzdiff_->Fill(zc - zg);
0149     hchi2_->Fill(chi2Ref, chi2Tar);
0150     hchi2oNdof_->Fill(chi2Ref / ndofRef, chi2Tar / ndofTar);
0151     hptv2_->Fill(vsoaRef.view()[sic].ptv2(), vsoaTar.view()[closestVtxidx].ptv2());
0152     hntrks_->Fill(ndofRef + 1, ndofTar + 1);
0153   }
0154   hnVertex_->Fill(nVerticesRef, nVerticesTar);
0155 }
0156
0157 //
0158 // -- Analyze
0159 //
0160 void SiPixelCompareVertices::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0161   // The default use case is to use vertices from Alpaka reconstructed on CPU and GPU;
0162   // The function is left templated if any other cases need to be added
0163   analyzeSeparate(tokenSoAVertexReferenceSoA_, tokenSoAVertexTargetSoA_, iEvent);
0164 }
0165
0166 //
0167 // -- Book Histograms
0168 //
0169 void SiPixelCompareVertices::bookHistograms(DQMStore::IBooker& ibooker,
0170                                             edm::Run const& iRun,
0171                                             edm::EventSetup const& iSetup) {
0172   ibooker.cd();
0173   ibooker.setCurrentFolder(topFolderName_);
0174
0175   // FIXME: all the 2D correlation plots are quite heavy in terms of memory consumption, so a as soon as DQM supports either TH2I or THnSparse
0176   // these should be moved to a less resource consuming format
0177   hnVertex_ = ibooker.book2I("nVertex", "# of Vertices;Reference;Target", 101, -0.5, 100.5, 101, -0.5, 100.5);
0178   hx_ = ibooker.book2I("vx", "Vertez x - Beamspot x;Reference;Target", 50, -0.1, 0.1, 50, -0.1, 0.1);
0179   hy_ = ibooker.book2I("vy", "Vertez y - Beamspot y;Reference;Target", 50, -0.1, 0.1, 50, -0.1, 0.1);
0180   hz_ = ibooker.book2I("vz", "Vertez z;Reference;Target", 30, -30., 30., 30, -30., 30.);
0181   hchi2_ = ibooker.book2I("chi2", "Vertex chi-squared;Reference;Target", 40, 0., 20., 40, 0., 20.);
0182   hchi2oNdof_ = ibooker.book2I("chi2oNdof", "Vertex chi-squared/Ndof;Reference;Target", 40, 0., 20., 40, 0., 20.);
0183   hptv2_ = ibooker.book2I("ptsq", "Vertex #sum (p_{T})^{2};Reference;Target", 200, 0., 200., 200, 0., 200.);
0184   hntrks_ = ibooker.book2I("ntrk", "#tracks associated;Reference;Target", 100, -0.5, 99.5, 100, -0.5, 99.5);
0185   hntrks_ = ibooker.book2I("ntrk", "#tracks associated;Reference;Target", 100, -0.5, 99.5, 100, -0.5, 99.5);
0186   hxdiff_ = ibooker.book1D("vxdiff", ";Vertex x difference (Reference - Target);#entries", 100, -0.001, 0.001);
0187   hydiff_ = ibooker.book1D("vydiff", ";Vertex y difference (Reference - Target);#entries", 100, -0.001, 0.001);
0188   hzdiff_ = ibooker.book1D("vzdiff", ";Vertex z difference (Reference - Target);#entries", 100, -2.5, 2.5);
0189 }
0190
0191 void SiPixelCompareVertices::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0192   // monitorpixelVertexSoA
0193   edm::ParameterSetDescription desc;
0194   desc.add<edm::InputTag>("pixelVertexReferenceSoA", edm::InputTag("pixelVerticesAlpakaSerial"));
0195   desc.add<edm::InputTag>("pixelVertexTargetSoA", edm::InputTag("pixelVerticesAlpaka"));
0196   desc.add<edm::InputTag>("beamSpotSrc", edm::InputTag("offlineBeamSpot"));
0197   desc.add<std::string>("topFolderName", "SiPixelHeterogeneous/PixelVertexCompareSoADeviceVSHost");
0198   desc.add<double>("dzCut", 1.);
0199   descriptions.addWithDefaultLabel(desc);
0200 }
0201
0202 // TODO: change module name to SiPixelCompareVerticesSoA when CUDA code is removed
0203 DEFINE_FWK_MODULE(SiPixelCompareVertices);