HGCalMapping/test/HGCalMappingESSourceTester.cc

0001 #include <cstdio>
0002 #include <chrono>
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008 #include "FWCore/Utilities/interface/StreamID.h"
0009 #include "FWCore/Framework/interface/ESWatcher.h"
0010 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012
0013 // #include "HeterogeneousCore/AlpakaCore/interface/alpaka/EDPutToken.h"
0014 // #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESGetToken.h"
0015 // #include "HeterogeneousCore/AlpakaCore/interface/alpaka/Event.h"
0016 // #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
0017 // #include "HeterogeneousCore/AlpakaInterface/interface/CopyToDevice.h"
0018
0019 #include "CondFormats/DataRecord/interface/HGCalElectronicsMappingRcd.h"
0020 #include "CondFormats/DataRecord/interface/HGCalDenseIndexInfoRcd.h"
0021 #include "CondFormats/HGCalObjects/interface/HGCalMappingModuleIndexer.h"
0022 #include "CondFormats/HGCalObjects/interface/HGCalMappingCellIndexer.h"
0023 #include "CondFormats/HGCalObjects/interface/HGCalMappingParameterHost.h"
0024 #include "Geometry/HGCalMapping/interface/HGCalMappingTools.h"
0025
0026 class HGCalMappingESSourceTester : public edm::one::EDAnalyzer<> {
0027 public:
0028   explicit HGCalMappingESSourceTester(const edm::ParameterSet&);
0029   static void fillDescriptions(edm::ConfigurationDescriptions&);
0030   std::map<uint32_t, uint32_t> mapGeoToElectronics(const hgcal::HGCalMappingModuleParamHost& modules,
0031                                                    const hgcal::HGCalMappingCellParamHost& cells,
0032                                                    bool geo2ele,
0033                                                    bool sipm);
0034
0035 private:
0036   void analyze(const edm::Event&, const edm::EventSetup&) override;
0037
0038   edm::ESWatcher<HGCalElectronicsMappingRcd> cfgWatcher_;
0039   edm::ESGetToken<HGCalMappingCellIndexer, HGCalElectronicsMappingRcd> cellIndexTkn_;
0040   edm::ESGetToken<hgcal::HGCalMappingCellParamHost, HGCalElectronicsMappingRcd> cellTkn_;
0041   edm::ESGetToken<HGCalMappingModuleIndexer, HGCalElectronicsMappingRcd> moduleIndexTkn_;
0042   edm::ESGetToken<hgcal::HGCalMappingModuleParamHost, HGCalElectronicsMappingRcd> moduleTkn_;
0043   edm::ESGetToken<hgcal::HGCalDenseIndexInfoHost, HGCalDenseIndexInfoRcd> denseIndexTkn_;
0044 };
0045
0046 //
0047 HGCalMappingESSourceTester::HGCalMappingESSourceTester(const edm::ParameterSet& iConfig)
0048     : cellIndexTkn_(esConsumes()),
0049       cellTkn_(esConsumes()),
0050       moduleIndexTkn_(esConsumes()),
0051       moduleTkn_(esConsumes()),
0052       denseIndexTkn_(esConsumes()) {}
0053
0054 //
0055 void HGCalMappingESSourceTester::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0056   // if the cfg didn't change there's nothing else to do
0057   if (!cfgWatcher_.check(iSetup))
0058     return;
0059
0060   //get cell indexers and SoA
0061   auto const& cellIdx = iSetup.getData(cellIndexTkn_);
0062   auto const& cells = iSetup.getData(cellTkn_);
0063   printf("[HGCalMappingIndexESSourceTester][analyze] Cell dense indexers and associated SoA retrieved for HGCAL\n");
0064   int nmodtypes = cellIdx.typeCodeIndexer_.size();
0065   printf("[HGCalMappingIndexESSourceTester][analyze] module cell indexer has %d module types\n", nmodtypes);
0066
0067   //printout and test the indexer contents for cells
0068   printf("[HGCalMappingIndexESSourceTester][analyze] Test indexer\n");
0069   uint32_t totOffset(0);
0070   for (size_t idx = 0; idx < cellIdx.di_.size(); idx++) {
0071     //check typecode exists
0072     auto typecode = cellIdx.getTypecodeFromEnum(idx);
0073     assert(cellIdx.typeCodeIndexer_.count(typecode) == 1);
0074
0075     //check that the current offset is consistent with the increment from cells from the previous module
0076     if (idx > 0) {
0077       uint32_t nch_prev = cellIdx.maxErx_[idx - 1] * cellIdx.maxChPerErx_;
0078       uint32_t delta_offset = cellIdx.offsets_[idx] - cellIdx.offsets_[idx - 1];
0079       assert(delta_offset == nch_prev);
0080     }
0081
0082     //assert offset is consistent with the accumulation
0083     auto off = cellIdx.offsets_[idx];
0084     assert(off == totOffset);
0085
0086     totOffset += cellIdx.maxErx_[idx] * cellIdx.maxChPerErx_;
0087
0088     //print
0089     printf("[HGCalMappingIndexESSourceTester][analyze][%s] has index(internal)=%ld #eRx=%d #cells=%d offset=%d\n",
0090            typecode.c_str(),
0091            idx,
0092            cellIdx.maxErx_[idx],
0093            cellIdx.di_[idx].getMaxIndex(),
0094            cellIdx.offsets_[idx]);
0095   }
0096
0097   assert(totOffset == cellIdx.maxDenseIndex());
0098   printf("[HGCalMappingIndexESSourceTester][analyze] SoA size for module cell mapping will be %d\n", totOffset);
0099
0100   //printout and test module cells SoA contents
0101   uint32_t ncells = cells.view().metadata().size();
0102   uint32_t validCells = 0;
0103   assert(ncells == cellIdx.maxDenseIndex());  //check for consistent size
0104   printf("[HGCalMappingIndexESSourceTester][analyze] Module cell mapping contents\n");
0105   for (uint32_t i = 0; i < ncells; i++) {
0106     auto icell = cells.view()[i];
0107     if (!cells.view()[i].valid())
0108       continue;
0109     validCells++;
0110     printf(
0111         "\t idx=%d isHD=%d iscalib=%d isSiPM=%d typeidx=%d chip=%d half=%d seq=%d rocpin=%d sensorcell=%d triglink=%d "
0112         "trigcell=%d i1=%d i2=%d t=%d trace=%f eleid=0x%x detid=0x%x\n",
0113         i,
0114         icell.isHD(),
0115         icell.iscalib(),
0116         icell.isSiPM(),
0117         icell.typeidx(),
0118         icell.chip(),
0119         icell.half(),
0120         icell.seq(),
0121         icell.rocpin(),
0122         icell.sensorcell(),
0123         icell.triglink(),
0124         icell.trigcell(),
0125         icell.i1(),
0126         icell.i2(),
0127         icell.t(),
0128         icell.trace(),
0129         icell.eleid(),
0130         icell.detid());
0131   }
0132
0133   //module mapping
0134   auto const& modulesIdx = iSetup.getData(moduleIndexTkn_);
0135   printf("[HGCalMappingIndexESSourceTester][analyze] Module indexer has FEDs=%d Types in sequences=%ld max idx=%d\n",
0136          modulesIdx.fedCount(),
0137          modulesIdx.getGlobalTypesCounter().size(),
0138          modulesIdx.maxModulesIndex());
0139   printf("[HGCalMappingIndexESSourceTester][analyze] FED Readout sequence\n");
0140   std::unordered_set<uint32_t> unique_modOffsets, unique_erxOffsets, unique_chDataOffsets;
0141   uint32_t totalmods(0);
0142   for (const auto& frs : modulesIdx.getFEDReadoutSequences()) {
0143     std::copy(
0144         frs.modOffsets_.begin(), frs.modOffsets_.end(), std::inserter(unique_modOffsets, unique_modOffsets.end()));
0145     std::copy(
0146         frs.erxOffsets_.begin(), frs.erxOffsets_.end(), std::inserter(unique_erxOffsets, unique_erxOffsets.end()));
0147     std::copy(frs.chDataOffsets_.begin(),
0148               frs.chDataOffsets_.end(),
0149               std::inserter(unique_chDataOffsets, unique_chDataOffsets.end()));
0150
0151     size_t nmods = frs.readoutTypes_.size();
0152     totalmods += nmods;
0153     printf("\t[FED %d] packs data from %ld ECON-Ds - readout types -> (offsets) :", frs.id, nmods);
0154     for (size_t i = 0; i < nmods; i++) {
0155       printf("\t%d -> (%d;%d;%d)", frs.readoutTypes_[i], frs.modOffsets_[i], frs.erxOffsets_[i], frs.chDataOffsets_[i]);
0156     }
0157     printf("\n");
0158   }
0159   //check that there are unique offsets per modules in the full system
0160   assert(unique_modOffsets.size() == totalmods);
0161   assert(unique_erxOffsets.size() == totalmods);
0162   assert(unique_chDataOffsets.size() == totalmods);
0163
0164   //get the module mapper SoA
0165   auto const& modules = iSetup.getData(moduleTkn_);
0166   int nmodules = modulesIdx.maxModulesIndex();
0167   int validModules = 0;
0168   assert(nmodules == modules.view().metadata().size());  //check for consistent size
0169   printf("[HGCalMappingIndexESSourceTester][analyze] Module mapping contents\n");
0170   for (int i = 0; i < nmodules; i++) {
0171     auto imod = modules.view()[i];
0172     if (!modules.view()[i].valid())
0173       continue;
0174     validModules++;
0175     printf(
0176         "\t idx=%d zside=%d isSiPM=%d plane=%d i1=%d i2=%d irot=%d celltype=%d typeidx=%d fedid=%d localfedid=%d "
0177         "captureblock=%d capturesblockidx=%d econdidx=%d eleid=0x%x detid=0x%d\n",
0178         i,
0179         imod.zside(),
0180         imod.isSiPM(),
0181         imod.plane(),
0182         imod.i1(),
0183         imod.i2(),
0184         imod.irot(),
0185         imod.celltype(),
0186         imod.typeidx(),
0187         imod.fedid(),
0188         imod.slinkidx(),
0189         imod.captureblock(),
0190         imod.captureblockidx(),
0191         imod.econdidx(),
0192         imod.eleid(),
0193         imod.detid());
0194   }
0195
0196   printf(
0197       "[HGCalMappingIndexESSourceTester][analyze] Module and cells maps retrieved for HGCAL %d/%d valid modules "
0198       "%d/%d valid cells",
0199       validModules,
0200       modules.view().metadata().size(),
0201       validCells,
0202       cells.view().metadata().size());
0203
0204   //test DetId to ElectronicsId mapping
0205   auto tmap = [](auto geo2ele, auto ele2geo) {
0206     //sizes must match
0207     assert(geo2ele.size() == ele2geo.size());
0208
0209     //test uniqueness of keys
0210     for (auto it : geo2ele) {
0211       assert(ele2geo.count(it.second) == 1);
0212       assert(ele2geo[it.second] == it.first);
0213     }
0214     for (auto it : ele2geo) {
0215       assert(geo2ele.count(it.second) == 1);
0216       assert(geo2ele[it.second] == it.first);
0217     }
0218
0219     return true;
0220   };
0221
0222   //apply test to Si
0223   std::map<uint32_t, uint32_t> sigeo2ele = this->mapGeoToElectronics(modules, cells, true, false);
0224   std::map<uint32_t, uint32_t> siele2geo = this->mapGeoToElectronics(modules, cells, false, false);
0225   printf("[HGCalMappingIndexESSourceTester][produce] Silicon electronics<->geometry map\n");
0226   printf("\tID maps ele2geo=%ld ID maps geo2ele=%ld\n", siele2geo.size(), sigeo2ele.size());
0227   tmap(sigeo2ele, siele2geo);
0228
0229   //apply test to SiPMs
0230   std::map<uint32_t, uint32_t> sipmgeo2ele = this->mapGeoToElectronics(modules, cells, true, true);
0231   std::map<uint32_t, uint32_t> sipmele2geo = this->mapGeoToElectronics(modules, cells, false, true);
0232   printf("[HGCalMappingIndexESSourceTester][produce] SiPM-on-tile electronics<->geometry map\n");
0233   printf("\tID maps ele2geo=%ld ID maps geo2ele=%ld\n", sipmele2geo.size(), sipmgeo2ele.size());
0234   tmap(sipmgeo2ele, sipmele2geo);
0235
0236   // Timing studies
0237   uint16_t fedid(175), econdidx(2), captureblockidx(0), chip(0), half(1), seq(12), rocpin(48);
0238   int zside(0), n_i(6000000);
0239   printf("[HGCalMappingIndexESSourceTester][produce]  Creating %d number of raw ElectronicsIds\n", n_i);
0240   uint32_t elecid(0);
0241   auto start = std::chrono::steady_clock::now();
0242   for (int i = 0; i < n_i; i++) {
0243     elecid = ::hgcal::mappingtools::getElectronicsId(zside, fedid, captureblockidx, econdidx, chip, half, seq);
0244   }
0245   auto stop = std::chrono::steady_clock::now();
0246   std::chrono::duration<float> elapsed = stop - start;
0247   printf("\tTime: %f seconds\n", elapsed.count());
0248
0249   HGCalElectronicsId eid(elecid);
0250   assert(eid.localFEDId() == fedid);
0251   assert((uint32_t)eid.captureBlock() == captureblockidx);
0252   assert((uint32_t)eid.econdIdx() == econdidx);
0253   assert((uint32_t)eid.econdeRx() == (uint32_t)(2 * chip + half));
0254   assert((uint32_t)eid.halfrocChannel() == seq);
0255   float eidrocpin = (uint32_t)eid.halfrocChannel() + 36 * ((uint32_t)eid.econdeRx() % 2) -
0256                     ((uint32_t)eid.halfrocChannel() > 17) * ((uint32_t)eid.econdeRx() % 2 + 1);
0257   assert(eidrocpin == rocpin);
0258
0259   int plane(1), u(-9), v(-6), celltype(2), celliu(3), celliv(7);
0260   printf("[HGCalMappingIndexESSourceTester][produce]  Creating %d number of raw  HGCSiliconDetIds\n", n_i);
0261   uint32_t geoid(0);
0262   start = std::chrono::steady_clock::now();
0263   for (int i = 0; i < n_i; i++) {
0264     geoid = ::hgcal::mappingtools::getSiDetId(zside, plane, u, v, celltype, celliu, celliv);
0265   }
0266   stop = std::chrono::steady_clock::now();
0267   elapsed = stop - start;
0268   printf("\tTime: %f seconds\n", elapsed.count());
0269   HGCSiliconDetId gid(geoid);
0270   assert(gid.type() == celltype);
0271   assert(gid.layer() == plane);
0272   assert(gid.cellU() == celliu);
0273   assert(gid.cellV() == celliv);
0274
0275   int modidx(0), cellidx(1);
0276   printf(
0277       "[HGCalMappingIndexESSourceTester][produce]  Creating %d number of raw ElectronicsIds from SoAs module idx: "
0278       "%d, cell idx: %d\n",
0279       n_i,
0280       modidx,
0281       cellidx);
0282   elecid = 0;
0283   start = std::chrono::steady_clock::now();
0284   for (int i = 0; i < n_i; i++) {
0285     elecid = modules.view()[modidx].eleid() + cells.view()[cellidx].eleid();
0286   }
0287   stop = std::chrono::steady_clock::now();
0288   elapsed = stop - start;
0289   printf("\tTime: %f seconds\n", elapsed.count());
0290   eid = HGCalElectronicsId(elecid);
0291   assert(eid.localFEDId() == modules.view()[modidx].fedid());
0292   assert((uint32_t)eid.captureBlock() == modules.view()[modidx].captureblockidx());
0293   assert((uint32_t)eid.econdIdx() == modules.view()[modidx].econdidx());
0294   assert((uint32_t)eid.halfrocChannel() == cells.view()[cellidx].seq());
0295   uint16_t econderx = cells.view()[cellidx].chip() * 2 + cells.view()[cellidx].half();
0296   assert((uint32_t)eid.econdeRx() == econderx);
0297
0298   printf(
0299       "[HGCalMappingIndexESSourceTester][produce] Creating %d number of raw HGCSiliconDetId from SoAs module idx: "
0300       "%d, cell idx: %d\n",
0301       n_i,
0302       modidx,
0303       cellidx);
0304   uint32_t detid(0);
0305   start = std::chrono::steady_clock::now();
0306   for (int i = 0; i < n_i; i++) {
0307     detid = modules.view()[modidx].detid() + cells.view()[cellidx].detid();
0308   }
0309   stop = std::chrono::steady_clock::now();
0310   elapsed = stop - start;
0311   printf("\tTime: %f seconds\n", elapsed.count());
0312   HGCSiliconDetId did(detid);
0313   assert(did.type() == modules.view()[modidx].celltype());
0314   assert(did.layer() == modules.view()[modidx].plane());
0315   assert(did.cellU() == cells.view()[cellidx].i1());
0316   assert(did.cellV() == cells.view()[cellidx].i2());
0317
0318   //test dense index token
0319   auto const& denseIndexInfo = iSetup.getData(denseIndexTkn_);
0320   printf("Retrieved %d dense index info\n", denseIndexInfo.view().metadata().size());
0321   int nindices = denseIndexInfo.view().metadata().size();
0322   printf("fedId fedReadoutSeq detId eleid modix cellidx channel x y z");
0323   for (int i = 0; i < nindices; i++) {
0324     auto row = denseIndexInfo.view()[i];
0325     printf("%d %d 0x%x 0x%x %d %d %d %f %f %f\n",
0326            row.fedId(),
0327            row.fedReadoutSeq(),
0328            row.detid(),
0329            row.eleid(),
0330            row.modInfoIdx(),
0331            row.cellInfoIdx(),
0332            row.chNumber(),
0333            row.x(),
0334            row.y(),
0335            row.z());
0336   }
0337 }
0338
0339 //
0340 void HGCalMappingESSourceTester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0341   edm::ParameterSetDescription desc;
0342   descriptions.addWithDefaultLabel(desc);
0343 }
0344
0345 //
0346 std::map<uint32_t, uint32_t> HGCalMappingESSourceTester::mapGeoToElectronics(
0347     const hgcal::HGCalMappingModuleParamHost& modules,
0348     const hgcal::HGCalMappingCellParamHost& cells,
0349     bool geo2ele,
0350     bool sipm) {
0351   //loop over different modules
0352   std::map<uint32_t, uint32_t> idmap;
0353   uint32_t ndups(0);
0354   printf("\n");
0355   for (int i = 0; i < modules.view().metadata().size(); i++) {
0356     auto imod = modules.view()[i];
0357
0358     if (!imod.valid())
0359       continue;
0360
0361     //require match to si or SiPM
0362     if (sipm != imod.isSiPM())
0363       continue;
0364
0365     if (imod.plane() == 0) {
0366       printf("WARNING: found plane=0 for i1=%d i2=%d siPM=%d @ index=%i\n", imod.i1(), imod.i2(), imod.isSiPM(), i);
0367       continue;
0368     }
0369
0370     //loop over cells in the module
0371     for (int j = 0; j < cells.view().metadata().size(); j++) {
0372       auto jcell = cells.view()[j];
0373
0374       //use only the information for cells which match the module type index
0375       if (jcell.typeidx() != imod.typeidx())
0376         continue;
0377
0378       //require that it's a valid cell
0379       if (!jcell.valid())
0380         continue;
0381
0382       //assert type of sensor
0383       assert(imod.isSiPM() == jcell.isSiPM());
0384
0385       // make sure the cell is part of the module and it's not a calibration cell
0386       if (jcell.t() != 1)
0387         continue;
0388
0389       uint32_t elecid = imod.eleid() + jcell.eleid();
0390
0391       uint32_t geoid(0);
0392
0393       if (sipm) {
0394         geoid = ::hgcal::mappingtools::getSiPMDetId(
0395             imod.zside(), imod.plane(), imod.i2(), imod.celltype(), jcell.i1(), jcell.i2());
0396       } else {
0397         geoid = imod.detid() + jcell.detid();
0398       }
0399
0400       if (geo2ele) {
0401         auto it = idmap.find(geoid);
0402         ndups += (it != idmap.end());
0403         if (!sipm && it != idmap.end() && imod.plane() <= 26) {
0404           HGCSiliconDetId detid(geoid);
0405           printf("WARNING duplicate found for plane=%d u=%d v=%d cellU=%d cellV=%d valid=%d -> detid=0x%x\n",
0406                  imod.plane(),
0407                  imod.i1(),
0408                  imod.i2(),
0409                  jcell.i1(),
0410                  jcell.i2(),
0411                  jcell.valid(),
0412                  detid.rawId());
0413         }
0414       }
0415       if (!geo2ele) {
0416         auto it = idmap.find(elecid);
0417         ndups += (it != idmap.end());
0418       }
0419
0420       //map
0421       idmap[geo2ele ? geoid : elecid] = geo2ele ? elecid : geoid;
0422     }
0423   }
0424
0425   if (ndups > 0) {
0426     printf("[HGCalMappingESSourceTester][mapGeoToElectronics] found %d duplicates with geo2ele=%d for sipm=%d\n",
0427            ndups,
0428            geo2ele,
0429            sipm);
0430   }
0431
0432   return idmap;
0433 }
0434
0435 // define this as a plug-in
0436 DEFINE_FWK_MODULE(HGCalMappingESSourceTester);