Newer
Older
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
R__LOAD_LIBRARY(libDDRec)
R__LOAD_LIBRARY(libfmt)
// Includes necessary for compilation
#include <DD4hep/Detector.h>
#include <DDRec/CellIDPositionConverter.h>
#include <DDSegmentation/BitFieldCoder.h>
#include <iostream>
#include <vector>
#include <ROOT/RDataFrame.hxx>
#include <dd4pod/CalorimeterHitData.h>
#include <fmt/format.h>
int hcal_segmentation() {
std::cout << "hello world" << std::endl;
// Get a handle to the global detector
dd4hep::Detector* det = &(dd4hep::Detector::getInstance());
// Load the detector from the compact file. Use the correct compact file for
// your simulation here. Also make sure your environment variables are setup
// correctly Example:
// - source /opt/detector/athena-canyonlands-v2.1/setup.sh
// - use "/opt/detector/athena-canyonlands-v2.1/share/athena/athena.xml"
det->fromCompact(
"/opt/detector/athena-canyonlands-v2.1/share/athena/athena.xml");
det->volumeManager();
det->apply("DD4hepVolumeManager", 0, 0);
// Create our cellID converter. This is useful to get cell coordinates. Not
// used below.
// auto cellIDConverter = std::make_shared<const
// dd4hep::rec::CellIDPositionConverter>(*det);
// Get our readout ID spec for HcalBarrelHits
auto idSpecBarrel = det->readout("HcalBarrelHits").idSpec();
auto idSpecEndcapN = det->readout("HcalEndcapNHits").idSpec();
auto idSpecEndcapP = det->readout("HcalEndcapPHits").idSpec();
// Get our cell ID decoder
auto idDecBarrel = idSpecBarrel.decoder();
auto idDecEndcapN = idSpecEndcapN.decoder();
auto idDecEndcapP = idSpecEndcapP.decoder();
// Get a bit mask for module, layer, x and y
// module: Which stave are we on?
// layer: layer within the module (stave)
// (x,): pad coordinates within the layer
//
const int moduleIdxBarrel = idDecBarrel->index("module");
const int layerIdxBarrel = idDecBarrel->index("layer");
const int xIdxBarrel = idDecBarrel->index("x");
const int yIdxBarrel = idDecBarrel->index("y");
// same for negative and positive endcaps (they don't use module)
const int layerIdxEndcapN = idDecEndcapN->index("layer");
const int xIdxEndcapN = idDecEndcapN->index("x");
const int yIdxEndcapN = idDecEndcapN->index("y");
const int layerIdxEndcapP = idDecEndcapP->index("layer");
const int xIdxEndcapP = idDecEndcapP->index("x");
const int yIdxEndcapP = idDecEndcapP->index("y");
// Simple lambda function that prints out the hit segmentation for the barrel
auto HcalBarrelPrinter =
[&](int event, const std::vector<dd4pod::CalorimeterHitData>& h) {
for (size_t i = 0; i < h.size(); ++i) {
// Only do hits above .1 MeV
if (h[i].energyDeposit < 1e-4) {
continue;
}
const int module = idDecBarrel->get(h[i].cellID, moduleIdxBarrel);
const int layer = idDecBarrel->get(h[i].cellID, layerIdxBarrel);
const int x = idDecBarrel->get(h[i].cellID, xIdxBarrel);
const int y = idDecBarrel->get(h[i].cellID, yIdxBarrel);
fmt::print(
"HcalBarrel event: {:3d} hit: {:3d} module: {:2d} layer: {:2d} x: {:4d} "
"y {:4d} energy: {:1.3e}\n",
event, i, module, layer, x, y, h[i].energyDeposit);
}
};
// Open the ROOT file from full simulation, change this by the name of
// your ROOT file
ROOT::RDataFrame d{"events", "dis.root"};
// For this demonstration we only use the first 50 events, not something
// you would want to do in general. For production, remove this line and
// change all occurences of d50 with d
auto d50 = d.Range(50);
// Print out hit segmentation data for the barren HCal
int event = 0;
d50.Define("event", [&]{return event++;}, {})
.Foreach(HcalBarrelPrinter, {"event", "HcalBarrelHits"});
// Let's do the same for the negative endcap
auto HcalEndcapNPrinter =
[&](int event, const std::vector<dd4pod::CalorimeterHitData>& h) {
for (size_t i = 0; i < h.size(); ++i) {
// Only do hits above .1 MeV
if (h[i].energyDeposit < 1e-4) {
continue;
}
const int layer = idDecEndcapN->get(h[i].cellID, layerIdxEndcapN);
const int x = idDecEndcapN->get(h[i].cellID, xIdxEndcapN);
const int y = idDecEndcapN->get(h[i].cellID, yIdxEndcapN);
fmt::print(
"HcalEndcapN event: {:3d} hit: {:3d} layer: {:2d} x: {:4d} "
"y {:4d} energy: {:1.3e}\n",
event, i, layer, x, y, h[i].energyDeposit);
}
};
event = 0;
d50.Define("event", [&]{return event++;}, {})
.Foreach(HcalEndcapNPrinter, {"event", "HcalEndcapNHits"});
return 0;
}