HGCalMapping/test/HGCalMappingESSourceTester.cc

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