102 auto const*
geo = lar::providerFrom<geo::Geometry>();
109 std::vector<art::Ptr<recob::Wire>> wires;
118 for (
auto const& channel : (*simChannelHandle)) {
124 if (ch1 % 1000 == 0)
mf::LogInfo(
"EvaluateROIEFF") << ch1;
125 int view =
geo->View(ch1);
126 auto const& timeSlices = channel.TDCIDEMap();
130 int tdctick_previous = -999;
134 std::vector<int> signal_starttick;
135 std::vector<int> signal_endtick;
136 std::vector<double> signal_energydeposits;
137 std::vector<double> signal_energy_max;
138 std::vector<double> signal_max_tdctick;
140 for (
auto const& timeSlice : timeSlices) {
141 auto const tpctime = timeSlice.first;
142 auto const& energyDeposits = timeSlice.second;
143 int tdctick =
static_cast<int>(clockData.TPCTDC2Tick(
double(tpctime)));
144 if (tdctick < 0 || tdctick >
int(detProp.ReadOutWindowSize()) - 1)
continue;
147 double e_deposit = 0;
148 for (
auto const& energyDeposit : energyDeposits) {
149 e_deposit += energyDeposit.energy;
152 if (tdctick_previous == -999) {
153 signal_starttick.push_back(tdctick);
158 else if (tdctick - tdctick_previous != 1) {
159 signal_starttick.push_back(tdctick);
160 signal_endtick.push_back(tdctick_previous);
161 signal_energydeposits.push_back(totalE);
162 signal_energy_max.push_back(maxE);
163 signal_max_tdctick.push_back(maxEtick);
168 else if (tdctick - tdctick_previous == 1) {
170 if (maxE < e_deposit) {
176 tdctick_previous = tdctick;
180 signal_endtick.push_back(tdctick_previous);
181 signal_energydeposits.push_back(totalE);
182 signal_energy_max.push_back(maxE);
183 signal_max_tdctick.push_back(maxEtick);
185 if (signal_starttick.size() == 0 ||
186 (signal_endtick.size() == 1 && signal_endtick.back() == -999))
189 for (
auto& wire : wires) {
190 if (wire->Channel() != ch1)
continue;
194 std::vector<float> vecADC =
203 for (
size_t s = 0; s < signal_starttick.size(); s++) {
205 bool IsSignalOutside =
true;
206 for (
const auto& range : signalROI.
get_ranges()) {
213 signal_max_tdctick[s],
216 IsSignalOutside =
false;
221 if (IsSignalOutside) {
223 h_energy[view]->Fill(signal_energy_max[s]);
229 for (
const auto& range : signalROI.
get_ranges()) {
236 signal_max_tdctick[s],
240 double maxadc_sig = 0;
241 int maxadc_tick = -99;
242 for (
int k = roiFirstBinTick; k < roiLastBinTick; k++) {
243 if (vecADC[k] > maxadc_sig) {
244 maxadc_sig = vecADC[k];
249 double maxE_roi = -999.;
250 double maxE_roi_tick = -999.;
252 if (maxE_roi < signal_energy_max[s]) {
253 maxE_roi = signal_energy_max[s];
254 maxE_roi_tick = signal_max_tdctick[s];
258 for (
size_t s2 = s + 1; s2 < signal_starttick.size(); s2++) {
261 signal_max_tdctick[s2],
264 if (maxE_roi < signal_energy_max[s2]) {
265 maxE_roi = signal_energy_max[s2];
266 maxE_roi_tick = signal_max_tdctick[s2];
268 if (s2 == signal_starttick.size() - 1) { s = s2; }
289 double roi_sig[3] = {0., 0., 0.};
290 double roi_total[3] = {0., 0., 0.};
292 for (
auto& wire : wires) {
294 if (
chStatus.IsBad(wirechannel))
continue;
296 int view = wire->View();
302 std::vector<float> vecADC =
305 roi_total[view] += signalROI.
get_ranges().size();
308 for (
const auto& range : signalROI.
get_ranges()) {
309 bool IsSigROI =
false;
314 double maxadc_sig = -99.;
315 for (
int k = roiFirstBinTick; k < roiLastBinTick; k++) {
321 for (
auto const& channel : (*simChannelHandle)) {
322 if (wirechannel != channel.Channel())
continue;
324 int tdctick_previous = -999;
328 std::vector<int> signal_starttick;
329 std::vector<int> signal_endtick;
330 std::vector<double> signal_energydeposits;
331 std::vector<double> signal_energy_max;
332 std::vector<double> signal_max_tdctick;
334 auto const& timeSlices = channel.TDCIDEMap();
336 for (
auto const& timeSlice : timeSlices) {
337 auto const tpctime = timeSlice.first;
338 auto const& energyDeposits = timeSlice.second;
339 int tdctick =
static_cast<int>(clockData.TPCTDC2Tick(
double(tpctime)));
340 if (tdctick < 0 || tdctick >
int(detProp.ReadOutWindowSize()) - 1)
continue;
342 double e_deposit = 0;
343 for (
auto const& energyDeposit : energyDeposits) {
344 e_deposit += energyDeposit.energy;
347 if (tdctick_previous == -999) {
348 signal_starttick.push_back(tdctick);
353 else if (tdctick - tdctick_previous != 1) {
354 signal_starttick.push_back(tdctick);
355 signal_endtick.push_back(tdctick_previous);
356 signal_energydeposits.push_back(totalE);
357 signal_energy_max.push_back(maxE);
358 signal_max_tdctick.push_back(maxEtick);
363 else if (tdctick - tdctick_previous == 1) {
365 if (maxE < e_deposit) {
371 tdctick_previous = tdctick;
375 signal_endtick.push_back(tdctick_previous);
376 signal_energydeposits.push_back(totalE);
377 signal_energy_max.push_back(maxE);
378 signal_max_tdctick.push_back(maxEtick);
380 if (signal_starttick.size() == 0 ||
381 (signal_endtick.size() == 1 && signal_endtick.back() == -999))
385 for (
size_t s = 0; s < signal_starttick.size(); s++) {
388 signal_max_tdctick[s],
397 h_roi[view]->Fill(0);
401 h_roi[view]->Fill(1);
408 for (
int i = 0; i < 3; i++) {
410 double purity = roi_sig[i] / roi_total[i];
411 if (purity == 1) purity = 1. - 1.e-6;
417 if (roi_total[0] + roi_total[1] + roi_total[2]) {
419 (roi_sig[0] + roi_sig[1] + roi_sig[2]) / (roi_total[0] + roi_total[1] + roi_total[2]);
420 if (purity == 1) purity = 1. - 1.e-6;
art::InputTag fSimulationProducerLabel
art::InputTag fWireProducerLabel
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
constexpr auto abs(T v)
Returns the absolute value of the argument.
TH1D * h1_tickdiff_max[3]
const range_list_t & get_ranges() const
Returns the internal list of non-void ranges.
int TDCtick_t
Type representing a TDC tick.
bool isSignalInROI(int starttick, int endtick, int maxtick, int roistart, int roiend)
unsigned int ChannelID_t
Type representing the ID of a readout channel.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Namespace collecting geometry-related classes utilities.