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