Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-06-28 02:36:13

0001 // TODO: change file name to SiPixelCompareRecHitsSoA.cc when CUDA code is removed
0002 
0003 #include "DQMServices/Core/interface/MonitorElement.h"
0004 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0005 #include "DQMServices/Core/interface/DQMStore.h"
0006 #include "DataFormats/Math/interface/approx_atan2.h"
0007 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0008 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0009 #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsHost.h"
0010 #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0018 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0019 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0020 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0021 
0022 // TODO: change class name to SiPixelCompareRecHitsSoA when CUDA code is removed
0023 template <typename T>
0024 class SiPixelCompareRecHits : public DQMEDAnalyzer {
0025 public:
0026   using HitsSoA = TrackingRecHitHost<T>;
0027 
0028   explicit SiPixelCompareRecHits(const edm::ParameterSet&);
0029   ~SiPixelCompareRecHits() override = default;
0030   void dqmBeginRun(const edm::Run&, const edm::EventSetup&) override;
0031   void bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) override;
0032   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0033   // analyzeSeparate is templated to accept distinct types of SoAs
0034   // The default use case is to use rechits from Alpaka reconstructed on CPU and GPU;
0035   template <typename U, typename V>
0036   void analyzeSeparate(U tokenRef, V tokenTar, const edm::Event& iEvent);
0037   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0038 
0039 private:
0040   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0041   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0042   // these two are both on Host but originally they have been produced on Host or on Device
0043   const edm::EDGetTokenT<HitsSoA> tokenSoAHitsReference_;
0044   const edm::EDGetTokenT<HitsSoA> tokenSoAHitsTarget_;
0045   const std::string topFolderName_;
0046   const float mind2cut_;
0047   static constexpr uint32_t invalidHit_ = std::numeric_limits<uint32_t>::max();
0048   static constexpr float micron_ = 10000.;
0049   const TrackerGeometry* tkGeom_ = nullptr;
0050   const TrackerTopology* tTopo_ = nullptr;
0051   MonitorElement* hnHits_;
0052   MonitorElement* hBchargeL_[4];  // max 4 barrel hits
0053   MonitorElement* hBsizexL_[4];
0054   MonitorElement* hBsizeyL_[4];
0055   MonitorElement* hBposxL_[4];
0056   MonitorElement* hBposyL_[4];
0057   MonitorElement* hFchargeD_[2][12];  // max 12 endcap disks
0058   MonitorElement* hFsizexD_[2][12];
0059   MonitorElement* hFsizeyD_[2][12];
0060   MonitorElement* hFposxD_[2][12];
0061   MonitorElement* hFposyD_[2][12];
0062   //differences
0063   MonitorElement* hBchargeDiff_;
0064   MonitorElement* hFchargeDiff_;
0065   MonitorElement* hBsizeXDiff_;
0066   MonitorElement* hFsizeXDiff_;
0067   MonitorElement* hBsizeYDiff_;
0068   MonitorElement* hFsizeYDiff_;
0069   MonitorElement* hBposXDiff_;
0070   MonitorElement* hFposXDiff_;
0071   MonitorElement* hBposYDiff_;
0072   MonitorElement* hFposYDiff_;
0073 };
0074 
0075 //
0076 // constructors
0077 //
0078 template <typename T>
0079 SiPixelCompareRecHits<T>::SiPixelCompareRecHits(const edm::ParameterSet& iConfig)
0080     : geomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()),
0081       topoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>()),
0082       tokenSoAHitsReference_(consumes(iConfig.getParameter<edm::InputTag>("pixelHitsReferenceSoA"))),
0083       tokenSoAHitsTarget_(consumes(iConfig.getParameter<edm::InputTag>("pixelHitsTargetSoA"))),
0084       topFolderName_(iConfig.getParameter<std::string>("topFolderName")),
0085       mind2cut_(iConfig.getParameter<double>("minD2cut")) {}
0086 
0087 //
0088 // Begin Run
0089 //
0090 template <typename T>
0091 void SiPixelCompareRecHits<T>::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0092   tkGeom_ = &iSetup.getData(geomToken_);
0093   tTopo_ = &iSetup.getData(topoToken_);
0094 }
0095 
0096 template <typename T>
0097 template <typename U, typename V>
0098 void SiPixelCompareRecHits<T>::analyzeSeparate(U tokenRef, V tokenTar, const edm::Event& iEvent) {
0099   const auto& rhsoaHandleRef = iEvent.getHandle(tokenRef);
0100   const auto& rhsoaHandleTar = iEvent.getHandle(tokenTar);
0101 
0102   if (not rhsoaHandleRef or not rhsoaHandleTar) {
0103     edm::LogWarning out("SiPixelCompareRecHits");
0104     if (not rhsoaHandleRef) {
0105       out << "reference rechits not found; ";
0106     }
0107     if (not rhsoaHandleTar) {
0108       out << "target rechits not found; ";
0109     }
0110     out << "the comparison will not run.";
0111     return;
0112   }
0113 
0114   auto const& rhsoaRef = *rhsoaHandleRef;
0115   auto const& rhsoaTar = *rhsoaHandleTar;
0116 
0117   auto const& soa2dRef = rhsoaRef.const_view();
0118   auto const& soa2dTar = rhsoaTar.const_view();
0119 
0120   uint32_t nHitsRef = soa2dRef.metadata().size();
0121   uint32_t nHitsTar = soa2dTar.metadata().size();
0122 
0123   hnHits_->Fill(nHitsRef, nHitsTar);
0124   auto detIds = tkGeom_->detUnitIds();
0125   for (uint32_t i = 0; i < nHitsRef; i++) {
0126     float minD = mind2cut_;
0127     uint32_t matchedHit = invalidHit_;
0128     uint16_t indRef = soa2dRef[i].detectorIndex();
0129     float xLocalRef = soa2dRef[i].xLocal();
0130     float yLocalRef = soa2dRef[i].yLocal();
0131     for (uint32_t j = 0; j < nHitsTar; j++) {
0132       if (soa2dTar.detectorIndex(j) == indRef) {
0133         float dx = xLocalRef - soa2dTar[j].xLocal();
0134         float dy = yLocalRef - soa2dTar[j].yLocal();
0135         float distance = dx * dx + dy * dy;
0136         if (distance < minD) {
0137           minD = distance;
0138           matchedHit = j;
0139         }
0140       }
0141     }
0142     DetId id = detIds[indRef];
0143     uint32_t chargeRef = soa2dRef[i].chargeAndStatus().charge;
0144     int16_t sizeXRef = std::ceil(float(std::abs(soa2dRef[i].clusterSizeX()) / 8.));
0145     int16_t sizeYRef = std::ceil(float(std::abs(soa2dRef[i].clusterSizeY()) / 8.));
0146     uint32_t chargeTar = 0;
0147     int16_t sizeXTar = -99;
0148     int16_t sizeYTar = -99;
0149     float xLocalTar = -999.;
0150     float yLocalTar = -999.;
0151     if (matchedHit != invalidHit_) {
0152       chargeTar = soa2dTar[matchedHit].chargeAndStatus().charge;
0153       sizeXTar = std::ceil(float(std::abs(soa2dTar[matchedHit].clusterSizeX()) / 8.));
0154       sizeYTar = std::ceil(float(std::abs(soa2dTar[matchedHit].clusterSizeY()) / 8.));
0155       xLocalTar = soa2dTar[matchedHit].xLocal();
0156       yLocalTar = soa2dTar[matchedHit].yLocal();
0157     }
0158     switch (id.subdetId()) {
0159       case PixelSubdetector::PixelBarrel:
0160         hBchargeL_[tTopo_->pxbLayer(id) - 1]->Fill(chargeRef, chargeTar);
0161         hBsizexL_[tTopo_->pxbLayer(id) - 1]->Fill(sizeXRef, sizeXTar);
0162         hBsizeyL_[tTopo_->pxbLayer(id) - 1]->Fill(sizeYRef, sizeYTar);
0163         hBposxL_[tTopo_->pxbLayer(id) - 1]->Fill(xLocalRef, xLocalTar);
0164         hBposyL_[tTopo_->pxbLayer(id) - 1]->Fill(yLocalRef, yLocalTar);
0165         hBchargeDiff_->Fill(chargeRef - chargeTar);
0166         hBsizeXDiff_->Fill(sizeXRef - sizeXTar);
0167         hBsizeYDiff_->Fill(sizeYRef - sizeYTar);
0168         hBposXDiff_->Fill(micron_ * (xLocalRef - xLocalTar));
0169         hBposYDiff_->Fill(micron_ * (yLocalRef - yLocalTar));
0170         break;
0171       case PixelSubdetector::PixelEndcap:
0172         hFchargeD_[tTopo_->pxfSide(id) - 1][tTopo_->pxfDisk(id) - 1]->Fill(chargeRef, chargeTar);
0173         hFsizexD_[tTopo_->pxfSide(id) - 1][tTopo_->pxfDisk(id) - 1]->Fill(sizeXRef, sizeXTar);
0174         hFsizeyD_[tTopo_->pxfSide(id) - 1][tTopo_->pxfDisk(id) - 1]->Fill(sizeYRef, sizeYTar);
0175         hFposxD_[tTopo_->pxfSide(id) - 1][tTopo_->pxfDisk(id) - 1]->Fill(xLocalRef, xLocalTar);
0176         hFposyD_[tTopo_->pxfSide(id) - 1][tTopo_->pxfDisk(id) - 1]->Fill(yLocalRef, yLocalTar);
0177         hFchargeDiff_->Fill(chargeRef - chargeTar);
0178         hFsizeXDiff_->Fill(sizeXRef - sizeXTar);
0179         hFsizeYDiff_->Fill(sizeYRef - sizeYTar);
0180         hFposXDiff_->Fill(micron_ * (xLocalRef - xLocalTar));
0181         hFposYDiff_->Fill(micron_ * (yLocalRef - yLocalTar));
0182         break;
0183     }
0184   }
0185 }
0186 
0187 //
0188 // -- Analyze
0189 //
0190 template <typename T>
0191 void SiPixelCompareRecHits<T>::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0192   // The default use case is to use vertices from Alpaka reconstructed on CPU and GPU;
0193   // The function is left templated if any other cases need to be added
0194   analyzeSeparate(tokenSoAHitsReference_, tokenSoAHitsTarget_, iEvent);
0195 }
0196 
0197 //
0198 // -- Book Histograms
0199 //
0200 template <typename T>
0201 void SiPixelCompareRecHits<T>::bookHistograms(DQMStore::IBooker& iBook,
0202                                               edm::Run const& iRun,
0203                                               edm::EventSetup const& iSetup) {
0204   iBook.cd();
0205   iBook.setCurrentFolder(topFolderName_);
0206 
0207   // clang-format off
0208   //Global
0209   hnHits_ = iBook.book2I("nHits", "ReferencevsTarget RecHits per event;#Reference RecHits;#Target RecHits", 200, 0, 5000,200, 0, 5000);
0210   //Barrel Layer
0211   for(unsigned int il=0;il<tkGeom_->numberOfLayers(PixelSubdetector::PixelBarrel);il++){
0212     hBchargeL_[il] = iBook.book2I(Form("recHitsBLay%dCharge",il+1), Form("ReferencevsTarget RecHits Charge Barrel Layer%d;Reference Charge;Target Charge",il+1), 250, 0, 100000, 250, 0, 100000);
0213     hBsizexL_[il] = iBook.book2I(Form("recHitsBLay%dSizex",il+1), Form("ReferencevsTarget RecHits SizeX Barrel Layer%d;Reference SizeX;Target SizeX",il+1), 30, 0, 30, 30, 0, 30);
0214     hBsizeyL_[il] = iBook.book2I(Form("recHitsBLay%dSizey",il+1), Form("ReferencevsTarget RecHits SizeY Barrel Layer%d;Reference SizeY;Target SizeY",il+1), 30, 0, 30, 30, 0, 30);
0215     hBposxL_[il] = iBook.book2D(Form("recHitsBLay%dPosx",il+1), Form("ReferencevsTarget RecHits x-pos in Barrel Layer%d;Reference pos x;Target pos x",il+1), 200, -5, 5, 200,-5,5);
0216     hBposyL_[il] = iBook.book2D(Form("recHitsBLay%dPosy",il+1), Form("ReferencevsTarget RecHits y-pos in Barrel Layer%d;Reference pos y;Target pos y",il+1), 200, -5, 5, 200,-5,5);
0217   }
0218   //Endcaps
0219   //Endcaps Disk
0220   for(int is=0;is<2;is++){
0221     int sign=is==0? -1:1;
0222     for(unsigned int id=0;id<tkGeom_->numberOfLayers(PixelSubdetector::PixelEndcap);id++){
0223       hFchargeD_[is][id] = iBook.book2I(Form("recHitsFDisk%+dCharge",id*sign+sign), Form("ReferencevsTarget RecHits Charge Endcaps Disk%+d;Reference Charge;Target Charge",id*sign+sign), 250, 0, 100000, 250, 0, 100000);
0224       hFsizexD_[is][id] = iBook.book2I(Form("recHitsFDisk%+dSizex",id*sign+sign), Form("ReferencevsTarget RecHits SizeX Endcaps Disk%+d;Reference SizeX;Target SizeX",id*sign+sign), 30, 0, 30, 30, 0, 30);
0225       hFsizeyD_[is][id] = iBook.book2I(Form("recHitsFDisk%+dSizey",id*sign+sign), Form("ReferencevsTarget RecHits SizeY Endcaps Disk%+d;Reference SizeY;Target SizeY",id*sign+sign), 30, 0, 30, 30, 0, 30);
0226       hFposxD_[is][id] = iBook.book2D(Form("recHitsFDisk%+dPosx",id*sign+sign), Form("ReferencevsTarget RecHits x-pos Endcaps Disk%+d;Reference pos x;Target pos x",id*sign+sign), 200, -5, 5, 200, -5, 5);
0227       hFposyD_[is][id] = iBook.book2D(Form("recHitsFDisk%+dPosy",id*sign+sign), Form("ReferencevsTarget RecHits y-pos Endcaps Disk%+d;Reference pos y;Target pos y",id*sign+sign), 200, -5, 5, 200, -5, 5);
0228     }
0229   }
0230   //1D differences
0231   hBchargeDiff_ = iBook.book1D("rechitChargeDiffBpix","Charge differnce of rechits in BPix; rechit charge difference (Reference - Target)", 101, -50.5, 50.5);
0232   hFchargeDiff_ = iBook.book1D("rechitChargeDiffFpix","Charge differnce of rechits in FPix; rechit charge difference (Reference - Target)", 101, -50.5, 50.5);
0233   hBsizeXDiff_ = iBook.book1D("rechitsizeXDiffBpix","SizeX difference of rechits in BPix; rechit sizex difference (Reference - Target)", 21, -10.5, 10.5);
0234   hFsizeXDiff_ = iBook.book1D("rechitsizeXDiffFpix","SizeX difference of rechits in FPix; rechit sizex difference (Reference - Target)", 21, -10.5, 10.5);
0235   hBsizeYDiff_ = iBook.book1D("rechitsizeYDiffBpix","SizeY difference of rechits in BPix; rechit sizey difference (Reference - Target)", 21, -10.5, 10.5);
0236   hFsizeYDiff_ = iBook.book1D("rechitsizeYDiffFpix","SizeY difference of rechits in FPix; rechit sizey difference (Reference - Target)", 21, -10.5, 10.5);
0237   hBposXDiff_ = iBook.book1D("rechitsposXDiffBpix","x-position difference of rechits in BPix; rechit x-pos difference (Reference - Target)", 1000, -10, 10);
0238   hFposXDiff_ = iBook.book1D("rechitsposXDiffFpix","x-position difference of rechits in FPix; rechit x-pos difference (Reference - Target)", 1000, -10, 10);
0239   hBposYDiff_ = iBook.book1D("rechitsposYDiffBpix","y-position difference of rechits in BPix; rechit y-pos difference (Reference - Target)", 1000, -10, 10);
0240   hFposYDiff_ = iBook.book1D("rechitsposYDiffFpix","y-position difference of rechits in FPix; rechit y-pos difference (Reference - Target)", 1000, -10, 10);
0241 }
0242 
0243 template<typename T>
0244 void SiPixelCompareRecHits<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0245   // monitorpixelRecHitsSoAAlpaka
0246   edm::ParameterSetDescription desc;
0247   desc.add<edm::InputTag>("pixelHitsReferenceSoA", edm::InputTag("siPixelRecHitsPreSplittingAlpakaSerial"));
0248   desc.add<edm::InputTag>("pixelHitsTargetSoA", edm::InputTag("siPixelRecHitsPreSplittingAlpaka"));
0249   desc.add<std::string>("topFolderName", "SiPixelHeterogeneous/PixelRecHitsCompareDeviceVSHost");
0250   desc.add<double>("minD2cut", 0.0001);
0251   descriptions.addWithDefaultLabel(desc);
0252 }
0253 
0254 using SiPixelPhase1CompareRecHits = SiPixelCompareRecHits<pixelTopology::Phase1>;
0255 using SiPixelPhase2CompareRecHits = SiPixelCompareRecHits<pixelTopology::Phase2>;
0256 using SiPixelHIonPhase1CompareRecHits = SiPixelCompareRecHits<pixelTopology::HIonPhase1>;
0257 
0258 #include "FWCore/Framework/interface/MakerMacros.h"
0259 // TODO: change module names to SiPixel*CompareRecHitsSoA when CUDA code is removed
0260 DEFINE_FWK_MODULE(SiPixelPhase1CompareRecHits);
0261 DEFINE_FWK_MODULE(SiPixelPhase2CompareRecHits);
0262 DEFINE_FWK_MODULE(SiPixelHIonPhase1CompareRecHits);
0263