1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
|
#include <cstdio>
#include <chrono>
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/Utilities/interface/StreamID.h"
#include "FWCore/Framework/interface/ESWatcher.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
// #include "HeterogeneousCore/AlpakaCore/interface/alpaka/EDPutToken.h"
// #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESGetToken.h"
// #include "HeterogeneousCore/AlpakaCore/interface/alpaka/Event.h"
// #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
// #include "HeterogeneousCore/AlpakaInterface/interface/CopyToDevice.h"
#include "CondFormats/DataRecord/interface/HGCalElectronicsMappingRcd.h"
#include "CondFormats/DataRecord/interface/HGCalDenseIndexInfoRcd.h"
#include "CondFormats/HGCalObjects/interface/HGCalMappingModuleIndexer.h"
#include "CondFormats/HGCalObjects/interface/HGCalMappingCellIndexer.h"
#include "CondFormats/HGCalObjects/interface/HGCalMappingParameterHost.h"
#include "Geometry/HGCalMapping/interface/HGCalMappingTools.h"
class HGCalMappingESSourceTester : public edm::one::EDAnalyzer<> {
public:
explicit HGCalMappingESSourceTester(const edm::ParameterSet&);
static void fillDescriptions(edm::ConfigurationDescriptions&);
std::map<uint32_t, uint32_t> mapGeoToElectronics(const hgcal::HGCalMappingModuleParamHost& modules,
const hgcal::HGCalMappingCellParamHost& cells,
bool geo2ele,
bool sipm);
private:
void analyze(const edm::Event&, const edm::EventSetup&) override;
edm::ESWatcher<HGCalElectronicsMappingRcd> cfgWatcher_;
edm::ESGetToken<HGCalMappingCellIndexer, HGCalElectronicsMappingRcd> cellIndexTkn_;
edm::ESGetToken<hgcal::HGCalMappingCellParamHost, HGCalElectronicsMappingRcd> cellTkn_;
edm::ESGetToken<HGCalMappingModuleIndexer, HGCalElectronicsMappingRcd> moduleIndexTkn_;
edm::ESGetToken<hgcal::HGCalMappingModuleParamHost, HGCalElectronicsMappingRcd> moduleTkn_;
edm::ESGetToken<hgcal::HGCalDenseIndexInfoHost, HGCalDenseIndexInfoRcd> denseIndexTkn_;
};
//
HGCalMappingESSourceTester::HGCalMappingESSourceTester(const edm::ParameterSet& iConfig)
: cellIndexTkn_(esConsumes()),
cellTkn_(esConsumes()),
moduleIndexTkn_(esConsumes()),
moduleTkn_(esConsumes()),
denseIndexTkn_(esConsumes()) {}
//
void HGCalMappingESSourceTester::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
// if the cfg didn't change there's nothing else to do
if (!cfgWatcher_.check(iSetup))
return;
//get cell indexers and SoA
auto const& cellIdx = iSetup.getData(cellIndexTkn_);
auto const& cells = iSetup.getData(cellTkn_);
printf("[HGCalMappingIndexESSourceTester][analyze] Cell dense indexers and associated SoA retrieved for HGCAL\n");
int nmodtypes = cellIdx.typeCodeIndexer_.size();
printf("[HGCalMappingIndexESSourceTester][analyze] module cell indexer has %d module types\n", nmodtypes);
//printout and test the indexer contents for cells
printf("[HGCalMappingIndexESSourceTester][analyze] Test indexer\n");
uint32_t totOffset(0);
for (size_t idx = 0; idx < cellIdx.di_.size(); idx++) {
//check typecode exists
auto typecode = cellIdx.getTypecodeFromEnum(idx);
assert(cellIdx.typeCodeIndexer_.count(typecode) == 1);
//check that the current offset is consistent with the increment from cells from the previous module
if (idx > 0) {
uint32_t nch_prev = cellIdx.maxErx_[idx - 1] * cellIdx.maxChPerErx_;
uint32_t delta_offset = cellIdx.offsets_[idx] - cellIdx.offsets_[idx - 1];
assert(delta_offset == nch_prev);
}
//assert offset is consistent with the accumulation
auto off = cellIdx.offsets_[idx];
assert(off == totOffset);
totOffset += cellIdx.maxErx_[idx] * cellIdx.maxChPerErx_;
//print
printf("[HGCalMappingIndexESSourceTester][analyze][%s] has index(internal)=%ld #eRx=%d #cells=%d offset=%d\n",
typecode.c_str(),
idx,
cellIdx.maxErx_[idx],
cellIdx.di_[idx].getMaxIndex(),
cellIdx.offsets_[idx]);
}
assert(totOffset == cellIdx.maxDenseIndex());
printf("[HGCalMappingIndexESSourceTester][analyze] SoA size for module cell mapping will be %d\n", totOffset);
//printout and test module cells SoA contents
uint32_t ncells = cells.view().metadata().size();
uint32_t validCells = 0;
assert(ncells == cellIdx.maxDenseIndex()); //check for consistent size
printf("[HGCalMappingIndexESSourceTester][analyze] Module cell mapping contents\n");
for (uint32_t i = 0; i < ncells; i++) {
auto icell = cells.view()[i];
if (!cells.view()[i].valid())
continue;
validCells++;
printf(
"\t idx=%d isHD=%d iscalib=%d isSiPM=%d typeidx=%d chip=%d half=%d seq=%d rocpin=%d sensorcell=%d triglink=%d "
"trigcell=%d i1=%d i2=%d t=%d trace=%f eleid=0x%x detid=0x%x\n",
i,
icell.isHD(),
icell.iscalib(),
icell.isSiPM(),
icell.typeidx(),
icell.chip(),
icell.half(),
icell.seq(),
icell.rocpin(),
icell.sensorcell(),
icell.triglink(),
icell.trigcell(),
icell.i1(),
icell.i2(),
icell.t(),
icell.trace(),
icell.eleid(),
icell.detid());
}
//module mapping
auto const& modulesIdx = iSetup.getData(moduleIndexTkn_);
printf("[HGCalMappingIndexESSourceTester][analyze] Module indexer has FEDs=%d Types in sequences=%ld max idx=%d\n",
modulesIdx.fedCount(),
modulesIdx.getGlobalTypesCounter().size(),
modulesIdx.maxModulesIndex());
printf("[HGCalMappingIndexESSourceTester][analyze] FED Readout sequence\n");
std::unordered_set<uint32_t> unique_modOffsets, unique_erxOffsets, unique_chDataOffsets;
uint32_t totalmods(0);
for (const auto& frs : modulesIdx.getFEDReadoutSequences()) {
std::copy(
frs.modOffsets_.begin(), frs.modOffsets_.end(), std::inserter(unique_modOffsets, unique_modOffsets.end()));
std::copy(
frs.erxOffsets_.begin(), frs.erxOffsets_.end(), std::inserter(unique_erxOffsets, unique_erxOffsets.end()));
std::copy(frs.chDataOffsets_.begin(),
frs.chDataOffsets_.end(),
std::inserter(unique_chDataOffsets, unique_chDataOffsets.end()));
size_t nmods = frs.readoutTypes_.size();
totalmods += nmods;
printf("\t[FED %d] packs data from %ld ECON-Ds - readout types -> (offsets) :", frs.id, nmods);
for (size_t i = 0; i < nmods; i++) {
printf("\t%d -> (%d;%d;%d)", frs.readoutTypes_[i], frs.modOffsets_[i], frs.erxOffsets_[i], frs.chDataOffsets_[i]);
}
printf("\n");
}
//check that there are unique offsets per modules in the full system
assert(unique_modOffsets.size() == totalmods);
assert(unique_erxOffsets.size() == totalmods);
assert(unique_chDataOffsets.size() == totalmods);
//get the module mapper SoA
auto const& modules = iSetup.getData(moduleTkn_);
int nmodules = modulesIdx.maxModulesIndex();
int validModules = 0;
assert(nmodules == modules.view().metadata().size()); //check for consistent size
printf("[HGCalMappingIndexESSourceTester][analyze] Module mapping contents\n");
for (int i = 0; i < nmodules; i++) {
auto imod = modules.view()[i];
if (!modules.view()[i].valid())
continue;
validModules++;
printf(
"\t idx=%d zside=%d isSiPM=%d plane=%d i1=%d i2=%d irot=%d celltype=%d typeidx=%d fedid=%d localfedid=%d "
"captureblock=%d capturesblockidx=%d econdidx=%d eleid=0x%x detid=0x%d\n",
i,
imod.zside(),
imod.isSiPM(),
imod.plane(),
imod.i1(),
imod.i2(),
imod.irot(),
imod.celltype(),
imod.typeidx(),
imod.fedid(),
imod.slinkidx(),
imod.captureblock(),
imod.captureblockidx(),
imod.econdidx(),
imod.eleid(),
imod.detid());
}
printf(
"[HGCalMappingIndexESSourceTester][analyze] Module and cells maps retrieved for HGCAL %d/%d valid modules "
"%d/%d valid cells",
validModules,
modules.view().metadata().size(),
validCells,
cells.view().metadata().size());
//test DetId to ElectronicsId mapping
auto tmap = [](auto geo2ele, auto ele2geo) {
//sizes must match
assert(geo2ele.size() == ele2geo.size());
//test uniqueness of keys
for (auto it : geo2ele) {
assert(ele2geo.count(it.second) == 1);
assert(ele2geo[it.second] == it.first);
}
for (auto it : ele2geo) {
assert(geo2ele.count(it.second) == 1);
assert(geo2ele[it.second] == it.first);
}
return true;
};
//apply test to Si
std::map<uint32_t, uint32_t> sigeo2ele = this->mapGeoToElectronics(modules, cells, true, false);
std::map<uint32_t, uint32_t> siele2geo = this->mapGeoToElectronics(modules, cells, false, false);
printf("[HGCalMappingIndexESSourceTester][produce] Silicon electronics<->geometry map\n");
printf("\tID maps ele2geo=%ld ID maps geo2ele=%ld\n", siele2geo.size(), sigeo2ele.size());
tmap(sigeo2ele, siele2geo);
//apply test to SiPMs
std::map<uint32_t, uint32_t> sipmgeo2ele = this->mapGeoToElectronics(modules, cells, true, true);
std::map<uint32_t, uint32_t> sipmele2geo = this->mapGeoToElectronics(modules, cells, false, true);
printf("[HGCalMappingIndexESSourceTester][produce] SiPM-on-tile electronics<->geometry map\n");
printf("\tID maps ele2geo=%ld ID maps geo2ele=%ld\n", sipmele2geo.size(), sipmgeo2ele.size());
tmap(sipmgeo2ele, sipmele2geo);
// Timing studies
uint16_t fedid(175), econdidx(2), captureblockidx(0), chip(0), half(1), seq(12), rocpin(48);
int zside(0), n_i(6000000);
printf("[HGCalMappingIndexESSourceTester][produce] Creating %d number of raw ElectronicsIds\n", n_i);
uint32_t elecid(0);
auto start = std::chrono::steady_clock::now();
for (int i = 0; i < n_i; i++) {
elecid = ::hgcal::mappingtools::getElectronicsId(zside, fedid, captureblockidx, econdidx, chip, half, seq);
}
auto stop = std::chrono::steady_clock::now();
std::chrono::duration<float> elapsed = stop - start;
printf("\tTime: %f seconds\n", elapsed.count());
HGCalElectronicsId eid(elecid);
assert(eid.localFEDId() == fedid);
assert((uint32_t)eid.captureBlock() == captureblockidx);
assert((uint32_t)eid.econdIdx() == econdidx);
assert((uint32_t)eid.econdeRx() == (uint32_t)(2 * chip + half));
assert((uint32_t)eid.halfrocChannel() == seq);
float eidrocpin = (uint32_t)eid.halfrocChannel() + 36 * ((uint32_t)eid.econdeRx() % 2) -
((uint32_t)eid.halfrocChannel() > 17) * ((uint32_t)eid.econdeRx() % 2 + 1);
assert(eidrocpin == rocpin);
int plane(1), u(-9), v(-6), celltype(2), celliu(3), celliv(7);
printf("[HGCalMappingIndexESSourceTester][produce] Creating %d number of raw HGCSiliconDetIds\n", n_i);
uint32_t geoid(0);
start = std::chrono::steady_clock::now();
for (int i = 0; i < n_i; i++) {
geoid = ::hgcal::mappingtools::getSiDetId(zside, plane, u, v, celltype, celliu, celliv);
}
stop = std::chrono::steady_clock::now();
elapsed = stop - start;
printf("\tTime: %f seconds\n", elapsed.count());
HGCSiliconDetId gid(geoid);
assert(gid.type() == celltype);
assert(gid.layer() == plane);
assert(gid.cellU() == celliu);
assert(gid.cellV() == celliv);
int modidx(0), cellidx(1);
printf(
"[HGCalMappingIndexESSourceTester][produce] Creating %d number of raw ElectronicsIds from SoAs module idx: "
"%d, cell idx: %d\n",
n_i,
modidx,
cellidx);
elecid = 0;
start = std::chrono::steady_clock::now();
for (int i = 0; i < n_i; i++) {
elecid = modules.view()[modidx].eleid() + cells.view()[cellidx].eleid();
}
stop = std::chrono::steady_clock::now();
elapsed = stop - start;
printf("\tTime: %f seconds\n", elapsed.count());
eid = HGCalElectronicsId(elecid);
assert(eid.localFEDId() == modules.view()[modidx].fedid());
assert((uint32_t)eid.captureBlock() == modules.view()[modidx].captureblockidx());
assert((uint32_t)eid.econdIdx() == modules.view()[modidx].econdidx());
assert((uint32_t)eid.halfrocChannel() == cells.view()[cellidx].seq());
uint16_t econderx = cells.view()[cellidx].chip() * 2 + cells.view()[cellidx].half();
assert((uint32_t)eid.econdeRx() == econderx);
printf(
"[HGCalMappingIndexESSourceTester][produce] Creating %d number of raw HGCSiliconDetId from SoAs module idx: "
"%d, cell idx: %d\n",
n_i,
modidx,
cellidx);
uint32_t detid(0);
start = std::chrono::steady_clock::now();
for (int i = 0; i < n_i; i++) {
detid = modules.view()[modidx].detid() + cells.view()[cellidx].detid();
}
stop = std::chrono::steady_clock::now();
elapsed = stop - start;
printf("\tTime: %f seconds\n", elapsed.count());
HGCSiliconDetId did(detid);
assert(did.type() == modules.view()[modidx].celltype());
assert(did.layer() == modules.view()[modidx].plane());
assert(did.cellU() == cells.view()[cellidx].i1());
assert(did.cellV() == cells.view()[cellidx].i2());
//test dense index token
auto const& denseIndexInfo = iSetup.getData(denseIndexTkn_);
printf("Retrieved %d dense index info\n", denseIndexInfo.view().metadata().size());
int nindices = denseIndexInfo.view().metadata().size();
printf("fedId fedReadoutSeq detId eleid modix cellidx channel x y z");
for (int i = 0; i < nindices; i++) {
auto row = denseIndexInfo.view()[i];
printf("%d %d 0x%x 0x%x %d %d %d %f %f %f\n",
row.fedId(),
row.fedReadoutSeq(),
row.detid(),
row.eleid(),
row.modInfoIdx(),
row.cellInfoIdx(),
row.chNumber(),
row.x(),
row.y(),
row.z());
}
}
//
void HGCalMappingESSourceTester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
descriptions.addWithDefaultLabel(desc);
}
//
std::map<uint32_t, uint32_t> HGCalMappingESSourceTester::mapGeoToElectronics(
const hgcal::HGCalMappingModuleParamHost& modules,
const hgcal::HGCalMappingCellParamHost& cells,
bool geo2ele,
bool sipm) {
//loop over different modules
std::map<uint32_t, uint32_t> idmap;
uint32_t ndups(0);
printf("\n");
for (int i = 0; i < modules.view().metadata().size(); i++) {
auto imod = modules.view()[i];
if (!imod.valid())
continue;
//require match to si or SiPM
if (sipm != imod.isSiPM())
continue;
if (imod.plane() == 0) {
printf("WARNING: found plane=0 for i1=%d i2=%d siPM=%d @ index=%i\n", imod.i1(), imod.i2(), imod.isSiPM(), i);
continue;
}
//loop over cells in the module
for (int j = 0; j < cells.view().metadata().size(); j++) {
auto jcell = cells.view()[j];
//use only the information for cells which match the module type index
if (jcell.typeidx() != imod.typeidx())
continue;
//require that it's a valid cell
if (!jcell.valid())
continue;
//assert type of sensor
assert(imod.isSiPM() == jcell.isSiPM());
// make sure the cell is part of the module and it's not a calibration cell
if (jcell.t() != 1)
continue;
uint32_t elecid = imod.eleid() + jcell.eleid();
uint32_t geoid(0);
if (sipm) {
geoid = ::hgcal::mappingtools::getSiPMDetId(
imod.zside(), imod.plane(), imod.i2(), imod.celltype(), jcell.i1(), jcell.i2());
} else {
geoid = imod.detid() + jcell.detid();
}
if (geo2ele) {
auto it = idmap.find(geoid);
ndups += (it != idmap.end());
if (!sipm && it != idmap.end() && imod.plane() <= 26) {
HGCSiliconDetId detid(geoid);
printf("WARNING duplicate found for plane=%d u=%d v=%d cellU=%d cellV=%d valid=%d -> detid=0x%x\n",
imod.plane(),
imod.i1(),
imod.i2(),
jcell.i1(),
jcell.i2(),
jcell.valid(),
detid.rawId());
}
}
if (!geo2ele) {
auto it = idmap.find(elecid);
ndups += (it != idmap.end());
}
//map
idmap[geo2ele ? geoid : elecid] = geo2ele ? elecid : geoid;
}
}
if (ndups > 0) {
printf("[HGCalMappingESSourceTester][mapGeoToElectronics] found %d duplicates with geo2ele=%d for sipm=%d\n",
ndups,
geo2ele,
sipm);
}
return idmap;
}
// define this as a plug-in
DEFINE_FWK_MODULE(HGCalMappingESSourceTester);
|