Skip to content

Commit d6330be

Browse files
gertvanhoeygertvanhoey
committed
Implement pile-up
1 parent d6c81f1 commit d6330be

File tree

4 files changed

+244
-7
lines changed

4 files changed

+244
-7
lines changed

‎core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp‎

Lines changed: 185 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -10,34 +10,214 @@
1010

1111
GateDigitizerPileupActor::GateDigitizerPileupActor(py::dict &user_info)
1212
: GateVDigitizerWithOutputActor(user_info, false) {
13-
14-
// Actions
1513
fActions.insert("EndOfEventAction");
14+
fActions.insert("EndOfRunAction");
1615
}
1716

1817
GateDigitizerPileupActor::~GateDigitizerPileupActor() = default;
1918

2019
void GateDigitizerPileupActor::InitializeUserInfo(py::dict &user_info) {
2120
GateVDigitizerWithOutputActor::InitializeUserInfo(user_info);
21+
22+
// Get time window parameter in ns.
23+
fTimeWindow = 1000.0; // default value
24+
if (py::len(user_info) > 0 && user_info.contains("time_window")) {
25+
fTimeWindow = DictGetDouble(user_info, "time_window");
26+
}
27+
fGroupVolumeDepth = -1;
28+
}
29+
30+
void GateDigitizerPileupActor::SetGroupVolumeDepth(const int depth) {
31+
fGroupVolumeDepth = depth;
2232
}
2333

2434
void GateDigitizerPileupActor::DigitInitialize(
2535
const std::vector<std::string> &attributes_not_in_filler) {
36+
// We handle all attributes explicitly by storing their values,
37+
// because we need to keep singles across consecutive events.
2638
auto a = attributes_not_in_filler;
39+
for (const auto &att_name : fInputDigiCollection->GetDigiAttributeNames()) {
40+
a.push_back(att_name);
41+
}
42+
2743
GateVDigitizerWithOutputActor::DigitInitialize(a);
44+
45+
// Get output attribute pointers
46+
fOutputEdepAttribute =
47+
fOutputDigiCollection->GetDigiAttribute("TotalEnergyDeposit");
48+
fOutputGlobalTimeAttribute =
49+
fOutputDigiCollection->GetDigiAttribute("GlobalTime");
50+
fOutputVolumeIDAttribute =
51+
fOutputDigiCollection->GetDigiAttribute("PreStepUniqueVolumeID");
52+
53+
// Set up pointers to track specific attributes
54+
auto &lr = fThreadLocalVDigitizerData.Get();
55+
auto &l = fThreadLocalData.Get();
56+
57+
lr.fInputIter = fInputDigiCollection->NewIterator();
58+
lr.fInputIter.TrackAttribute("TotalEnergyDeposit", &l.edep);
59+
lr.fInputIter.TrackAttribute("GlobalTime", &l.time);
60+
lr.fInputIter.TrackAttribute("PreStepUniqueVolumeID", &l.volID);
2861
}
2962

3063
void GateDigitizerPileupActor::BeginOfRunAction(const G4Run *run) {
3164
GateVDigitizerWithOutputActor::BeginOfRunAction(run);
3265
}
3366

34-
void GateDigitizerPileupActor::EndOfEventAction(const G4Event * /*unused*/) {
67+
void GateDigitizerPileupActor::StoreAttributeValues(
68+
threadLocalT::PileupGroup &group, size_t index) {
69+
// Clear previous values
70+
group.stored_attributes.clear();
71+
72+
// Store values from all attributes in the input collection
73+
for (auto *att : fInputDigiCollection->GetDigiAttributes()) {
74+
auto name = att->GetDigiAttributeName();
75+
auto type = att->GetDigiAttributeType();
76+
77+
// Skip the attributes we handle explicitly
78+
if (name == "TotalEnergyDeposit" || name == "GlobalTime" ||
79+
name == "PreStepUniqueVolumeID") {
80+
continue;
81+
}
82+
83+
// Store value based on type
84+
switch (type) {
85+
case 'D': // double
86+
group.stored_attributes[name] = att->GetDValues()[index];
87+
break;
88+
case 'I': // int
89+
group.stored_attributes[name] = att->GetIValues()[index];
90+
break;
91+
case 'L': // int64_t
92+
group.stored_attributes[name] = att->GetLValues()[index];
93+
break;
94+
case 'S': // string
95+
group.stored_attributes[name] = att->GetSValues()[index];
96+
break;
97+
case '3': // G4ThreeVector
98+
group.stored_attributes[name] = att->Get3Values()[index];
99+
break;
100+
case 'U': // GateUniqueVolumeID::Pointer
101+
group.stored_attributes[name] = att->GetUValues()[index];
102+
break;
103+
}
104+
}
105+
}
106+
107+
void GateDigitizerPileupActor::EndOfEventAction(const G4Event *) {
108+
ProcessPileup();
109+
110+
// Create output singles from piled-up groups of input singles.
111+
112+
auto &l = fThreadLocalData.Get();
113+
114+
for (auto &[volume_id, groups] : l.volume_groups) {
115+
if (!groups.empty()) {
116+
// Handle all groups, except the last one, because the last group may
117+
// still get contributions from upcoming events.
118+
for (auto it = groups.begin(); it != std::prev(groups.end()); ++it) {
119+
FillAttributeValues(*it);
120+
}
121+
groups.erase(groups.begin(), std::prev(groups.end()));
122+
}
123+
}
124+
}
125+
126+
void GateDigitizerPileupActor::ProcessPileup() {
127+
35128
auto &lr = fThreadLocalVDigitizerData.Get();
129+
auto &l = fThreadLocalData.Get();
36130
auto &iter = lr.fInputIter;
131+
37132
iter.GoToBegin();
38133
while (!iter.IsAtEnd()) {
39-
auto &i = lr.fInputIter.fIndex;
40-
lr.fDigiAttributeFiller->Fill(i);
134+
135+
const auto current_time = *l.time;
136+
const auto current_edep = *l.edep;
137+
const auto current_vol =
138+
l.volID->get()->GetIdUpToDepthAsHash(fGroupVolumeDepth);
139+
const auto current_index = iter.fIndex;
140+
141+
bool added_to_existing_group = false;
142+
auto &groups = l.volume_groups[current_vol];
143+
if (groups.size() > 0) {
144+
auto &group = groups.back();
145+
if (std::abs(current_time - group.first_time) <= fTimeWindow) {
146+
// Accumulate deposited energy of all singles in the same time window.
147+
group.total_edep += current_edep;
148+
// Keep the attributes of the highest energy single.
149+
if (current_edep > group.highest_edep) {
150+
group.highest_edep = current_edep;
151+
group.time = current_time;
152+
// Store all other attribute values from this single.
153+
StoreAttributeValues(group, current_index);
154+
}
155+
added_to_existing_group = true;
156+
}
157+
}
158+
159+
// If not added to an existing group, create a new group.
160+
if (!added_to_existing_group) {
161+
typename threadLocalT::PileupGroup new_group;
162+
new_group.total_edep = current_edep;
163+
new_group.highest_edep = current_edep;
164+
new_group.first_time = current_time;
165+
new_group.time = current_time;
166+
new_group.volume_id = *l.volID;
167+
// Store all other attribute values from this single
168+
StoreAttributeValues(new_group, current_index);
169+
groups.push_back(new_group);
170+
}
171+
41172
iter++;
42173
}
43174
}
175+
176+
void GateDigitizerPileupActor::EndOfRunAction(const G4Run *) {
177+
178+
auto &l = fThreadLocalData.Get();
179+
180+
// Output any unfinished groups that are still present.
181+
for (auto &[volume_id, groups] : l.volume_groups) {
182+
if (!groups.empty()) {
183+
for (auto &group : groups) {
184+
FillAttributeValues(group);
185+
}
186+
groups.erase(groups.begin(), groups.end());
187+
}
188+
}
189+
190+
// Make sure everything is output into the root file.
191+
fOutputDigiCollection->FillToRootIfNeeded(true);
192+
}
193+
194+
void GateDigitizerPileupActor::FillAttributeValues(
195+
const threadLocalT::PileupGroup &group) {
196+
197+
fOutputEdepAttribute->FillDValue(group.total_edep);
198+
fOutputGlobalTimeAttribute->FillDValue(group.time);
199+
fOutputVolumeIDAttribute->FillUValue(group.volume_id);
200+
201+
// Fill all other stored attributes
202+
for (const auto &[name, value] : group.stored_attributes) {
203+
auto *att = fOutputDigiCollection->GetDigiAttribute(name);
204+
std::visit(
205+
[att](auto &&arg) {
206+
using T = std::decay_t<decltype(arg)>;
207+
if constexpr (std::is_same_v<T, double>) {
208+
att->FillDValue(arg);
209+
} else if constexpr (std::is_same_v<T, int>) {
210+
att->FillIValue(arg);
211+
} else if constexpr (std::is_same_v<T, int64_t>) {
212+
att->FillLValue(arg);
213+
} else if constexpr (std::is_same_v<T, std::string>) {
214+
att->FillSValue(arg);
215+
} else if constexpr (std::is_same_v<T, G4ThreeVector>) {
216+
att->Fill3Value(arg);
217+
} else if constexpr (std::is_same_v<T, GateUniqueVolumeID::Pointer>) {
218+
att->FillUValue(arg);
219+
}
220+
},
221+
value);
222+
}
223+
}

‎core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h‎

Lines changed: 42 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,9 @@
1111
#include "GateVDigitizerWithOutputActor.h"
1212
#include <G4Cache.hh>
1313
#include <G4Navigator.hh>
14+
#include <map>
1415
#include <pybind11/stl.h>
16+
#include <variant>
1517

1618
namespace py = pybind11;
1719

@@ -36,13 +38,52 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor {
3638
// Called every time an Event ends (all threads)
3739
void EndOfEventAction(const G4Event *event) override;
3840

41+
// Called every time a Run ends (all threads)
42+
void EndOfRunAction(const G4Run *run) override;
43+
44+
void SetGroupVolumeDepth(int depth);
45+
3946
protected:
4047
void DigitInitialize(
4148
const std::vector<std::string> &attributes_not_in_filler) override;
4249

50+
// User parameters
51+
int fGroupVolumeDepth;
52+
double fTimeWindow;
53+
54+
// Output attribute pointers
55+
GateVDigiAttribute *fOutputEdepAttribute{};
56+
GateVDigiAttribute *fOutputGlobalTimeAttribute{};
57+
GateVDigiAttribute *fOutputVolumeIDAttribute{};
58+
59+
void ProcessPileup();
60+
4361
// During computation (thread local)
44-
struct threadLocalT {};
62+
struct threadLocalT {
63+
double *edep;
64+
double *time;
65+
GateUniqueVolumeID::Pointer *volID;
66+
67+
// Storage for piled up singles
68+
struct PileupGroup {
69+
double highest_edep;
70+
double total_edep;
71+
double first_time;
72+
double time;
73+
GateUniqueVolumeID::Pointer volume_id;
74+
75+
using AttributeValue =
76+
std::variant<double, int, int64_t, std::string, G4ThreeVector,
77+
GateUniqueVolumeID::Pointer>;
78+
std::map<std::string, AttributeValue> stored_attributes;
79+
};
80+
81+
std::map<uint64_t, std::vector<PileupGroup>> volume_groups;
82+
};
4583
G4Cache<threadLocalT> fThreadLocalData;
84+
85+
void StoreAttributeValues(threadLocalT::PileupGroup &group, size_t index);
86+
void FillAttributeValues(const threadLocalT::PileupGroup &group);
4687
};
4788

4889
#endif // GateDigitizerPileupActor_h

‎core/opengate_core/opengate_lib/digitizer/pyGateDigitizerPileupActor.cpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,5 +17,7 @@ void init_GateDigitizerPileupActor(py::module &m) {
1717
py::class_<GateDigitizerPileupActor,
1818
std::unique_ptr<GateDigitizerPileupActor, py::nodelete>,
1919
GateVDigitizerWithOutputActor>(m, "GateDigitizerPileupActor")
20-
.def(py::init<py::dict &>());
20+
.def(py::init<py::dict &>())
21+
.def("SetGroupVolumeDepth",
22+
&GateDigitizerPileupActor::SetGroupVolumeDepth);
2123
}

‎opengate/actors/digitizers.py‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -564,6 +564,12 @@ class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor)
564564
"doc": "FIXME",
565565
},
566566
),
567+
"group_volume": (
568+
None,
569+
{
570+
"doc": "FIXME",
571+
},
572+
),
567573
}
568574

569575
def __init__(self, *args, **kwargs):
@@ -579,6 +585,14 @@ def initialize(self):
579585
self.InitializeUserInfo(self.user_info)
580586
self.InitializeCpp()
581587

588+
def set_group_by_depth(self):
589+
depth = -1
590+
if self.user_info.group_volume is not None:
591+
depth = self.simulation.volume_manager.get_volume(
592+
self.user_info.group_volume
593+
).volume_depth_in_tree
594+
self.SetGroupVolumeDepth(depth)
595+
582596
def StartSimulationAction(self):
583597
DigitizerBase.StartSimulationAction(self)
584598
g4.GateDigitizerPileupActor.StartSimulationAction(self)

0 commit comments

Comments
 (0)