97 auto const*
geo = lar::providerFrom<geo::Geometry>();
104 std::vector<art::Ptr<recob::Wire>> wires;
113 for (
auto const&
channel : (*simChannelHandle)) {
119 if (ch1 % 1000 == 0)
mf::LogInfo(
"EvaluateROIEFF") << ch1;
120 int view =
geo->View(ch1);
121 auto const& timeSlices =
channel.TDCIDEMap();
125 int tdctick_previous = -999;
129 std::vector<int> signal_starttick;
130 std::vector<int> signal_endtick;
131 std::vector<double> signal_energydeposits;
132 std::vector<double> signal_energy_max;
133 std::vector<double> signal_max_tdctick;
135 for (
auto const& timeSlice : timeSlices) {
136 auto const tpctime = timeSlice.first;
137 auto const& energyDeposits = timeSlice.second;
138 int tdctick =
static_cast<int>(clockData.TPCTDC2Tick(
double(tpctime)));
139 if (tdctick < 0 || tdctick >
int(detProp.ReadOutWindowSize()) - 1)
continue;
142 double e_deposit = 0;
143 for (
auto const& energyDeposit : energyDeposits) {
144 e_deposit += energyDeposit.energy;
147 if (tdctick_previous == -999) {
148 signal_starttick.push_back(tdctick);
153 else if (tdctick - tdctick_previous != 1) {
154 signal_starttick.push_back(tdctick);
155 signal_endtick.push_back(tdctick_previous);
156 signal_energydeposits.push_back(totalE);
157 signal_energy_max.push_back(maxE);
158 signal_max_tdctick.push_back(maxEtick);
163 else if (tdctick - tdctick_previous == 1) {
165 if (maxE < e_deposit) {
171 tdctick_previous = tdctick;
175 signal_endtick.push_back(tdctick_previous);
176 signal_energydeposits.push_back(totalE);
177 signal_energy_max.push_back(maxE);
178 signal_max_tdctick.push_back(maxEtick);
180 if (signal_starttick.size() == 0 ||
181 (signal_endtick.size() == 1 && signal_endtick.back() == -999))
184 for (
auto& wire : wires) {
185 if (wire->Channel() != ch1)
continue;
189 std::vector<float> vecADC =
198 for (
size_t s = 0;
s < signal_starttick.size();
s++) {
200 bool IsSignalOutside =
true;
201 for (
const auto& range : signalROI.
get_ranges()) {
208 signal_max_tdctick[s],
211 IsSignalOutside =
false;
216 if (IsSignalOutside) {
218 h_energy[view]->Fill(signal_energy_max[
s]);
224 for (
const auto& range : signalROI.
get_ranges()) {
231 signal_max_tdctick[s],
235 double maxadc_sig = 0;
236 int maxadc_tick = -99;
237 for (
int k = roiFirstBinTick;
k < roiLastBinTick;
k++) {
238 if (vecADC[
k] > maxadc_sig) {
239 maxadc_sig = vecADC[
k];
244 double maxE_roi = -999.;
245 double maxE_roi_tick = -999.;
247 if (maxE_roi < signal_energy_max[s]) {
248 maxE_roi = signal_energy_max[
s];
249 maxE_roi_tick = signal_max_tdctick[
s];
253 for (
size_t s2 = s + 1; s2 < signal_starttick.size(); s2++) {
256 signal_max_tdctick[s2],
259 if (maxE_roi < signal_energy_max[s2]) {
260 maxE_roi = signal_energy_max[s2];
261 maxE_roi_tick = signal_max_tdctick[s2];
263 if (s2 == signal_starttick.size() - 1) { s = s2; }
284 double roi_sig[3] = {0., 0., 0.};
285 double roi_total[3] = {0., 0., 0.};
287 for (
auto& wire : wires) {
289 if (
chStatus.IsBad(wirechannel))
continue;
291 int view = wire->View();
297 std::vector<float> vecADC =
300 roi_total[view] += signalROI.
get_ranges().size();
303 for (
const auto& range : signalROI.
get_ranges()) {
304 bool IsSigROI =
false;
309 double maxadc_sig = -99.;
310 for (
int k = roiFirstBinTick;
k < roiLastBinTick;
k++) {
316 for (
auto const&
channel : (*simChannelHandle)) {
317 if (wirechannel !=
channel.Channel())
continue;
319 int tdctick_previous = -999;
323 std::vector<int> signal_starttick;
324 std::vector<int> signal_endtick;
325 std::vector<double> signal_energydeposits;
326 std::vector<double> signal_energy_max;
327 std::vector<double> signal_max_tdctick;
329 auto const& timeSlices =
channel.TDCIDEMap();
331 for (
auto const& timeSlice : timeSlices) {
332 auto const tpctime = timeSlice.first;
333 auto const& energyDeposits = timeSlice.second;
334 int tdctick =
static_cast<int>(clockData.TPCTDC2Tick(
double(tpctime)));
335 if (tdctick < 0 || tdctick >
int(detProp.ReadOutWindowSize()) - 1)
continue;
337 double e_deposit = 0;
338 for (
auto const& energyDeposit : energyDeposits) {
339 e_deposit += energyDeposit.energy;
342 if (tdctick_previous == -999) {
343 signal_starttick.push_back(tdctick);
348 else if (tdctick - tdctick_previous != 1) {
349 signal_starttick.push_back(tdctick);
350 signal_endtick.push_back(tdctick_previous);
351 signal_energydeposits.push_back(totalE);
352 signal_energy_max.push_back(maxE);
353 signal_max_tdctick.push_back(maxEtick);
358 else if (tdctick - tdctick_previous == 1) {
360 if (maxE < e_deposit) {
366 tdctick_previous = tdctick;
370 signal_endtick.push_back(tdctick_previous);
371 signal_energydeposits.push_back(totalE);
372 signal_energy_max.push_back(maxE);
373 signal_max_tdctick.push_back(maxEtick);
375 if (signal_starttick.size() == 0 ||
376 (signal_endtick.size() == 1 && signal_endtick.back() == -999))
380 for (
size_t s = 0; s < signal_starttick.size(); s++) {
383 signal_max_tdctick[s],
392 h_roi[view]->Fill(0);
396 h_roi[view]->Fill(1);
403 for (
int i = 0; i < 3; i++) {
405 double purity = roi_sig[i] / roi_total[i];
406 if (purity == 1) purity = 1. - 1.e-6;
412 if (roi_total[0] + roi_total[1] + roi_total[2]) {
414 (roi_sig[0] + roi_sig[1] + roi_sig[2]) / (roi_total[0] + roi_total[1] + roi_total[2]);
415 if (purity == 1) purity = 1. - 1.e-6;
art::InputTag fSimulationProducerLabel
art::InputTag fWireProducerLabel
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
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)
LArSoft geometry interface.