Skip to content

Commit ecf9fd7

Browse files
mhkim-anlmhkim-anl
committed
Update for PR
1 parent 4d621f2 commit ecf9fd7

File tree

4 files changed

+230
-95
lines changed

4 files changed

+230
-95
lines changed

‎src/algorithms/digi/PulseDigi.cc‎

Lines changed: 115 additions & 95 deletions
Original file line numberDiff line numberDiff line change
@@ -21,103 +21,123 @@
2121
#include "PulseDigi.h"
2222

2323
class HGCROCRawSample {
24-
public:
25-
HGCROCRawSample(std::size_t n_meas)
26-
: meas_types(n_meas, 0), amplitudes(n_meas, 0), TOAs(n_meas, 0), TOTs(n_meas, 0) {}
27-
28-
void setAmplitude(std::size_t idx, double amp) { amplitudes[idx] = amp; }
29-
void setTOA(std::size_t idx, double toa) {
30-
meas_types[idx] = 1;
31-
TOAs[idx] = toa;
32-
}
33-
void setTOT(std::size_t idx, double tot) {
34-
meas_types[idx] = 2;
35-
TOTs[idx] = tot;
36-
}
37-
38-
const std::vector<uint8_t>& getMeasTypes() const { return meas_types; }
39-
const std::vector<double>& getAmplitudes() const { return amplitudes; }
40-
const std::vector<double>& getTOAs() const { return TOAs; }
41-
const std::vector<double>& getTOTs() const { return TOTs; }
42-
43-
private:
44-
// meas_type 1: Pulse crossed below the threshold.
45-
// meas_type 2: Pulse didn't cross the threshold.
46-
// meas_type 0: Neither of the above cases.
47-
std::vector<uint8_t> meas_types;
48-
std::vector<double> amplitudes;
49-
std::vector<double> TOAs;
50-
std::vector<double> TOTs;
24+
public:
25+
struct RawEntry {
26+
double adc{0};
27+
double toa{0};
28+
double tot{0};
29+
bool totProgress{false};
30+
bool totComplete{false};
31+
};
32+
33+
HGCROCRawSample(std::size_t n_samp) { rawEntries.resize(n_samp); }
34+
35+
void setAmplitude(std::size_t idx, double amp) { rawEntries[idx].adc = amp; }
36+
void setTOA(std::size_t idx, double toa) { rawEntries[idx].toa = toa; }
37+
void setTOT(std::size_t idx, double cross_t) { rawEntries[idx].tot = cross_t - rawEntries[idx].toa; }
38+
void setTotProgress(std::size_t idx) { rawEntries[idx].totProgress = true; }
39+
void setTotComplete(std::size_t idx) { rawEntries[idx].totComplete = true; }
40+
41+
const std::vector<RawEntry>& getEntries() const { return rawEntries; }
42+
43+
private:
44+
std::vector<RawEntry> rawEntries;
5145
};
5246

5347
namespace eicrecon {
5448

55-
void PulseDigi::init() {}
56-
57-
void PulseDigi::process(const PulseDigi::Input& input, const PulseDigi::Output& output) const {
58-
const auto [in_pulses] = input;
59-
auto [out_digi_hits] = output;
60-
61-
for (const auto& pulse : *in_pulses) {
62-
double pulse_t = pulse.getTime();
63-
double pulse_dt = pulse.getInterval();
64-
std::size_t n_amps = pulse.getAmplitude().size();
65-
66-
// Estimate the number of samples.
67-
const std::size_t timeIdx_begin =
68-
static_cast<std::size_t>(std::floor(pulse_t / m_cfg.time_window));
69-
const std::size_t timeIdx_end = static_cast<std::size_t>(
70-
std::floor((pulse_t + (n_amps - 1) * pulse_dt) / m_cfg.time_window));
71-
72-
HGCROCRawSample raw_sample(timeIdx_end - timeIdx_begin + 1);
73-
74-
int sample_tick = std::llround(m_cfg.sample_period / pulse_dt);
75-
int adc_counter =
76-
std::llround((timeIdx_begin * m_cfg.sample_period + m_cfg.adc_phase - pulse_t) / pulse_dt);
77-
78-
bool tot_progress = false;
79-
bool tot_complete = false;
80-
std::size_t toaIdx = 0;
81-
82-
for (std::size_t i = 0; i < n_amps; i++) {
83-
double t = pulse_t + i * pulse_dt;
84-
const std::size_t sampleIdx =
85-
static_cast<std::size_t>(std::floor(t / m_cfg.sample_period)) - timeIdx_begin;
86-
87-
adc_counter++;
88-
// Measure amplitudes
89-
if (adc_counter == sample_tick) {
90-
raw_sample.setAmplitude(sampleIdx, pulse.getAmplitude()[i]);
91-
adc_counter = 0;
92-
}
93-
94-
// Measure crossing point for TOA
95-
if (!tot_progress && pulse.getAmplitude()[i] > m_cfg.toa_thres) {
96-
toaIdx = sampleIdx;
97-
raw_sample.setTOA(sampleIdx,
98-
get_crossing_time(m_cfg.threshold, t, pulse_dt, pulse.getAmplitude()[i],
99-
pulse.getAmplitude()[i - 1]));
100-
tot_progress = true;
101-
tot_complete = false;
102-
}
103-
104-
// Measure crossing point for TOT
105-
if (tot_progress && !tot_complete && pulse.getAmplitude()[i] < m_cfg.threshold) {
106-
raw_sample.setTOT(toaIdx,
107-
get_crossing_time(m_cfg.threshold, t, pulse_dt, pulse.getAmplitude()[i],
108-
pulse.getAmplitude()[i - 1]));
109-
tot_complete = true;
110-
tot_progess = false;
111-
}
112-
}
113-
}
114-
} // PulseDigi:process
115-
116-
double PulseDigi::get_crossing_time(double thres, double t1, double t2, double amp1,
117-
double amp2) const {
118-
double numerator = (thres - amp2) * (t2 - t1);
119-
double denominator = amp2 - amp1;
120-
double added = t2;
121-
return (numerator / denominator) + added;
122-
}
49+
void PulseDigi::init() {}
50+
51+
void PulseDigi::process(const PulseDigi::Input& input, const PulseDigi::Output& output) const {
52+
const auto [in_pulses] = input;
53+
auto [out_digi_hits] = output;
54+
55+
for (const auto& pulse : *in_pulses) {
56+
double pulse_t = pulse.getTime();
57+
double pulse_dt = pulse.getInterval();
58+
std::size_t n_amps = pulse.getAmplitude().size();
59+
60+
// Estimate the number of samples.
61+
const std::size_t timeIdx_begin =
62+
static_cast<std::size_t>(std::floor(pulse_t / m_cfg.time_window));
63+
const std::size_t timeIdx_end = static_cast<std::size_t>(
64+
std::floor((pulse_t + (n_amps - 1) * pulse_dt) / m_cfg.time_window));
65+
66+
HGCROCRawSample raw_sample(timeIdx_end - timeIdx_begin + 1);
67+
68+
// For ADC, amplitude is measured with a fixed phase.
69+
// This was reproduced by sample_tick and adc_counter as follows.
70+
// Amplitude is measured whenever adc_counter reaches sample_tick.
71+
int sample_tick = std::llround(m_cfg.time_window / pulse_dt);
72+
int adc_counter =
73+
std::llround((timeIdx_begin * m_cfg.time_window + m_cfg.adc_phase - pulse_t) / pulse_dt);
74+
75+
bool tot_progress = false;
76+
bool tot_complete = false;
77+
std::size_t toaIdx = 0;
78+
79+
for (std::size_t i = 0; i < n_amps; i++) {
80+
double t = pulse_t + i * pulse_dt;
81+
const std::size_t sampleIdx =
82+
static_cast<std::size_t>(std::floor(t / m_cfg.time_window)) - timeIdx_begin;
83+
84+
adc_counter++;
85+
86+
// Measure amplitudes for ADC
87+
if (adc_counter == sample_tick) {
88+
raw_sample.setAmplitude(sampleIdx, pulse.getAmplitude()[i]);
89+
adc_counter = 0;
90+
if(tot_progress) raw_sample.setTotProgress(sampleIdx);
91+
}
92+
93+
// Measure up-crossing time for TOA
94+
if (!tot_progress && pulse.getAmplitude()[i] > m_cfg.toa_thres) {
95+
toaIdx = sampleIdx;
96+
raw_sample.setTOA(sampleIdx,
97+
get_crossing_time(m_cfg.toa_thres, pulse_dt, t, pulse.getAmplitude()[i],
98+
pulse.getAmplitude()[i - 1]));
99+
tot_progress = true;
100+
tot_complete = false;
101+
raw_sample.setTotProgress(sampleIdx);
102+
}
103+
104+
// Measure down-crossing time for TOT
105+
if (tot_progress && !tot_complete && pulse.getAmplitude()[i] < m_cfg.tot_thres) {
106+
raw_sample.setTOT(toaIdx,
107+
get_crossing_time(m_cfg.tot_thres, pulse_dt, t, pulse.getAmplitude()[i],
108+
pulse.getAmplitude()[i - 1]));
109+
tot_progress = false;
110+
tot_complete = true;
111+
raw_sample.setTotComplete(sampleIdx);
112+
}
113+
}
114+
115+
// Fill HGCROCSamples and RawHGCROCHit
116+
auto out_digi_hit = out_digi_hits->create();
117+
out_digi_hit.setCellID(pulse.getCellID());
118+
119+
const auto& entries = raw_sample.getEntries();
120+
121+
for (const auto& entry : entries) {
122+
edm4eic::HGCROCSample sample;
123+
auto adc = std::max(std::llround(entry.adc / m_cfg.dyRangeADC * m_cfg.capADC), 0LL);
124+
sample.ADC = adc > m_cfg.capADC ? m_cfg.capADC : adc;
125+
auto toa = std::max(std::llround(entry.toa / m_cfg.dyRangeTOA * m_cfg.capTOA), 0LL);
126+
sample.timeOfArrival = toa > m_cfg.capTOA ? m_cfg.capTOA : toa;
127+
auto tot = std::max(std::llround(entry.tot / m_cfg.dyRangeTOT * m_cfg.capTOT), 0LL);
128+
sample.timeOverThreshold = tot > m_cfg.capTOT ? m_cfg.capTOT : tot;
129+
sample.TOTInProgress = entry.totProgress;
130+
sample.TOTComplete = entry.totComplete;
131+
out_digi_hit.addToSamples(sample);
132+
}
133+
}
134+
} // PulseDigi:process
135+
136+
double PulseDigi::get_crossing_time(double thres, double dt, double t, double amp1,
137+
double amp2) const {
138+
double numerator = (thres - amp2) * dt;
139+
double denominator = amp2 - amp1;
140+
double added = t;
141+
return (numerator / denominator) + added;
142+
}
123143
} // namespace eicrecon

‎src/algorithms/digi/PulseDigi.h‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
// SPDX-License-Identifier: LGPL-3.0-or-later
2+
// Copyright (C) 2025 Minho Kim
3+
4+
#pragma once
5+
6+
#include <algorithms/algorithm.h>
7+
#include <edm4eic/SimPulseCollection.h>
8+
#include <edm4eic/RawHGCROCHitCollection.h>
9+
#include <edm4eic/HGCROCSample.h>
10+
#include <stdint.h>
11+
#include <functional>
12+
#include <string>
13+
#include <string_view>
14+
15+
#include "algorithms/digi/PulseDigiConfig.h"
16+
#include "algorithms/interfaces/WithPodConfig.h"
17+
18+
namespace eicrecon {
19+
20+
using PulseDigiAlgorithm = algorithms::Algorithm<
21+
algorithms::Input<edm4eic::SimPulseCollection>,
22+
algorithms::Output<edm4eic::RawHGCROCHitCollection>>;
23+
24+
class PulseDigi : public PulseDigiAlgorithm, public WithPodConfig<PulseDigiConfig> {
25+
26+
public:
27+
PulseDigi(std::string_view name)
28+
: PulseDigiAlgorithm{name, {"InputPulses"}, {"OutputRawHGCROCHits"}, {}} {}
29+
virtual void init() final;
30+
void process(const Input&, const Output&) const;
31+
32+
private:
33+
edm4eic::HGCROCSample sample;
34+
35+
private:
36+
double get_crossing_time(double thres, double dt, double t,
37+
double amp1, double amp2) const;
38+
};
39+
40+
} // namespace eicrecon

‎src/algorithms/digi/PulseDigiConfig.h‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
// SPDX-License-Identifier: LGPL-3.0-or-later
2+
// Copyright (C) 2025 Minho Kim
3+
4+
#pragma once
5+
6+
#include <edm4eic/unit_system.h>
7+
8+
namespace eicrecon {
9+
10+
struct PulseDigiConfig {
11+
12+
// Variables for HGCROC measurement
13+
double time_window{25 * edm4eic::unit::ns};
14+
double adc_phase{10 * edm4eic::unit::ns};
15+
double toa_thres{1};
16+
double tot_thres{1};
17+
18+
// Variables for digitization
19+
unsigned int capADC{1024};
20+
double dyRangeADC{1};
21+
unsigned int capTOA{1024};
22+
double dyRangeTOA{1 * edm4eic::unit::ns};
23+
unsigned int capTOT{1024};
24+
double dyRangeTOT{1 * edm4eic::unit::ns};
25+
};
26+
27+
} // namespace eicrecon

‎src/factories/digi/PulseDigi_factory.h‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
// SPDX-License-Identifier: LGPL-3.0-or-later
2+
// Copyright (C) 2025 Minho Kim
3+
4+
#pragma once
5+
6+
#include "algorithms/digi/PulseDigi.h"
7+
#include "services/algorithms_init/AlgorithmsInit_service.h"
8+
#include "extensions/jana/JOmniFactory.h"
9+
10+
namespace eicrecon {
11+
12+
class PulseDigi_factory : public JOmniFactory<PulseDigi_factory, PulseDigiConfig> {
13+
14+
public:
15+
using AlgoT = eicrecon::PulseDigi;
16+
17+
private:
18+
std::unique_ptr<AlgoT> m_algo;
19+
20+
PodioInput<edm4eic::SimPulse> m_pulse_input{this};
21+
PodioOutput<edm4eic::RawHGCROCHit> m_digi_output{this};
22+
23+
ParameterRef<double> m_time_window{this, "time_window", config().time_window};
24+
ParameterRef<double> m_adc_phase{this, "adc_phase", config().adc_phase};
25+
ParameterRef<double> m_toa_thres{this, "toa_thres", config().toa_thres};
26+
ParameterRef<double> m_tot_thres{this, "tot_thres", config().tot_thres};
27+
ParameterRef<unsigned int> m_capADC{this, "capADC", config().capADC};
28+
ParameterRef<double> m_dyRangeADC{this, "dyRangeADC", config().dyRangeADC};
29+
ParameterRef<unsigned int> m_capTOA{this, "capTOA", config().capTOA};
30+
ParameterRef<double> m_dyRangeTOA{this, "dyRangeTOA", config().dyRangeTOA};
31+
ParameterRef<unsigned int> m_capTOT{this, "capTOT", config().capTOT};
32+
ParameterRef<double> m_dyRangeTOT{this, "dyRangeTOT", config().dyRangeTOT};
33+
34+
Service<AlgorithmsInit_service> m_algorithmsInit{this};
35+
36+
public:
37+
void Configure() {
38+
m_algo = std::make_unique<AlgoT>(GetPrefix());
39+
m_algo->level(static_cast<algorithms::LogLevel>(logger()->level()));
40+
m_algo->applyConfig(config());
41+
m_algo->init();
42+
}
43+
44+
void Process(int32_t /* run_number */, uint64_t /* event_number */) {
45+
m_algo->process({m_pulse_input()}, {m_digi_output().get()});
46+
}
47+
};
48+
} // namespace eicrecon

0 commit comments

Comments
 (0)