15 #include "art_root_io/TFileService.h" 16 #include "cetlib_except/exception.h" 61 using MaxMinPair = std::pair<Waveform::const_iterator, Waveform::const_iterator>;
65 Waveform::const_iterator>;
69 Waveform::const_iterator,
76 Waveform::const_iterator)
const;
78 Waveform::const_iterator)
const;
81 Waveform::const_iterator
findStartTick(Waveform::const_iterator,
82 Waveform::const_iterator)
const;
83 Waveform::const_iterator
findStopTick(Waveform::const_iterator, Waveform::const_iterator)
const;
144 <<
"Cannot fill histograms when multiple threads configured, please set " 145 "fOutputHistograms to false or change number of threads to 1\n";
150 <<
"Cannot write output waveforms when multiple threads configured, please set " 151 "fOutputHistograms to false or change number of threads to 1\n";
174 dir.make<TH1F>(Form(
"DStopStart_%1zu",
fPlane),
";Delta Stop/Start;", 100, 0., 100.);
176 dir.make<TH1F>(Form(
"DMaxTMinT_%1zu",
fPlane),
";Delta Max/Min Tick;", 100, 0., 100.);
178 dir.make<TH1F>(Form(
"DMaxDMinD_%1zu",
fPlane),
";Delta Max/Min Deriv;", 200, 0., 100.);
180 dir.make<TH1F>(Form(
"MaxErosion_%1zu",
fPlane),
";Max Erosion;", 200, -50., 150.);
182 dir.make<TH1F>(Form(
"MaxDilation_%1zu",
fPlane),
";Max Dilation;", 200, -50., 150.);
184 dir.make<TH1F>(Form(
"MaxDilEroRat_%1zu",
fPlane),
";Max Dil/Ero;", 200, -1., 1.);
191 const recob::Wire::RegionsOfInterest_t::datarange_t& dataRange,
192 const size_t roiStartTick,
193 const size_t channel,
194 const size_t eventCount,
203 const Waveform& waveform = dataRange.data();
206 fWaveformTool->triangleSmooth(rawDerivativeVec, derivativeVec);
238 for (
auto& hc : hitCandidateVec) {
240 if (startCand >= 0) hc.startTick = std::max(hc.startTick,
size_t(startCand));
247 for (
auto& hitCandidate : hitCandidateVec) {
248 size_t centerIdx = hitCandidate.hitCenter;
250 hitCandidate.hitHeight = waveform.at(centerIdx);
257 size_t plane = wids[0].Plane;
258 size_t cryo = wids[0].Cryostat;
259 size_t tpc = wids[0].TPC;
260 size_t wire = wids[0].Wire;
266 Form(
"Event%04zu/c%1zuT%1zuP%1zu/Wire_%05zu", eventCount, cryo, tpc, plane, wire));
268 size_t waveformSize = waveform.size();
269 size_t waveStart = dataRange.begin_index();
272 dir.make<TProfile>(Form(
"HWfm_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
279 TProfile* derivHist =
280 dir.make<TProfile>(Form(
"HDer_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
287 TProfile* erosionHist =
288 dir.make<TProfile>(Form(
"HEro_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
295 TProfile* dilationHist =
296 dir.make<TProfile>(Form(
"HDil_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
303 TProfile* candHitHist =
304 dir.make<TProfile>(Form(
"HCan_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
311 TProfile* maxDerivHist =
312 dir.make<TProfile>(Form(
"HMax_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
319 TProfile* strtStopHist =
320 dir.make<TProfile>(Form(
"HSSS_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
329 for (
size_t idx = 0; idx < waveform.size(); idx++) {
330 waveHist->Fill(roiStartTick + idx, waveform.at(idx));
331 derivHist->Fill(roiStartTick + idx, derivativeVec.at(idx));
332 erosionHist->Fill(roiStartTick + idx, erosionVec.at(idx));
333 dilationHist->Fill(roiStartTick + idx, dilationVec.at(idx));
337 for (
const auto& hitCandidate : hitCandidateVec) {
338 candHitHist->Fill(hitCandidate.hitCenter, hitCandidate.hitHeight);
339 maxDerivHist->Fill(hitCandidate.maxTick, hitCandidate.maxDerivative);
340 maxDerivHist->Fill(hitCandidate.minTick, hitCandidate.minDerivative);
341 strtStopHist->Fill(hitCandidate.startTick, waveform.at(hitCandidate.startTick));
342 strtStopHist->Fill(hitCandidate.stopTick, waveform.at(hitCandidate.stopTick));
351 for (
const auto& hitCandidate : hitCandidateVec) {
352 fDStopStartHist->Fill(hitCandidate.stopTick - hitCandidate.startTick, 1.);
359 std::max_element(dilationVec.begin(), dilationVec.end());
361 std::max_element(erosionVec.begin(), erosionVec.end());
365 if (
std::abs(*maxDilationItr) > 0.) dilEroRat = *maxErosionItr / *maxDilationItr;
381 const size_t roiStartTick,
382 float dilationThreshold,
389 int ticksInInputWaveform = std::distance(derivStartItr, derivStopItr);
395 float maxVal = *maxItr;
398 if (maxVal < dilationThreshold)
return;
400 int maxBin = std::distance(dilationStartItr, maxItr);
407 float firstDerivValue = -1.;
411 while (firstDerItr != derivStartItr) {
413 if (*erosionItr-- < erosionCutValue) {
417 if (*firstDerItr * firstDerivValue <= 0. && firstDerivValue > 0.)
break;
420 firstDerivValue = *firstDerItr--;
424 int hitRegionStart = std::distance(derivStartItr, firstDerItr);
430 float lastDerivValue = 1.;
432 erosionItr = erosionStartItr + maxBin;
435 while (lastDerItr != derivStopItr) {
436 if (*erosionItr++ <= erosionCutValue) {
439 if (*lastDerItr * lastDerivValue <= 0. && lastDerivValue < 0.)
break;
442 lastDerivValue = *lastDerItr++;
446 int hitRegionStop = std::distance(derivStartItr, lastDerItr);
451 derivStartItr + hitRegionStart,
453 erosionStartItr + hitRegionStart,
455 dilationStartItr + hitRegionStart,
462 derivStartItr + hitRegionStop,
463 roiStartTick + hitRegionStart,
472 erosionStartItr + hitRegionStop,
474 dilationStartItr + hitRegionStop,
476 roiStartTick + hitRegionStop,
485 const size_t roiStartTick,
487 float dPeakThreshold,
496 startItr, stopItr, dTicksThreshold, dPeakThreshold, candHitParamsVec)) {
499 for (
const auto& tuple : candHitParamsVec) {
509 std::min_element(maxItr, minItr, [](
const auto&
left,
const auto&
right) {
510 return std::fabs(
left) < std::fabs(
right);
514 if (2 * std::distance(peakItr, minItr) < std::distance(maxItr, peakItr))
516 else if (2 * std::distance(maxItr, peakItr) < std::distance(peakItr, minItr))
520 size_t hitCandidateStartTick = roiStartTick + std::distance(startItr, candStartItr);
522 if (!hitCandidateVec.empty()) {
523 int deltaTicks = hitCandidateStartTick - hitCandidateVec.back().stopTick;
525 if (deltaTicks > 0) {
526 hitCandidateStartTick -= deltaTicks / 2;
527 hitCandidateVec.back().stopTick += deltaTicks / 2;
531 hitCandidate.
startTick = hitCandidateStartTick;
532 hitCandidate.
stopTick = roiStartTick + std::distance(startItr, candStopItr);
533 hitCandidate.
maxTick = roiStartTick + std::distance(startItr, maxItr);
534 hitCandidate.
minTick = roiStartTick + std::distance(startItr, minItr);
535 hitCandidate.
maxDerivative = maxItr != stopItr ? *maxItr : 0.;
536 hitCandidate.
minDerivative = minItr != stopItr ? *minItr : 0.;
537 hitCandidate.
hitCenter = roiStartTick + std::distance(startItr, peakItr) + 0.5;
542 hitCandidateVec.push_back(hitCandidate);
638 float dPeakThreshold,
642 bool foundCandidate(
false);
644 int dTicks = std::distance(startItr, stopItr);
653 MaxMinPair minMaxPair = std::minmax_element(startItr, stopItr);
659 if (std::fabs(*maxItr) > std::fabs(*minItr))
664 int deltaTicks = std::distance(maxItr, minItr);
665 float range = *maxItr - *minItr;
667 if (deltaTicks < 2)
return foundCandidate;
670 if (deltaTicks >= dTicksThreshold && range > dPeakThreshold) foundCandidate =
true;
682 startItr, candStartItr, dTicksThreshold, dPeakThreshold, candHitParamsVec);
685 candHitParamsVec.emplace_back(candStartItr, maxItr, minItr, candStopItr);
689 candStopItr, stopItr, dTicksThreshold, dPeakThreshold, candHitParamsVec);
691 return foundCandidate || prevTicks || postTicks;
695 const recob::Wire::RegionsOfInterest_t::datarange_t&,
700 if (hitCandidateVec.empty())
return;
708 size_t lastStopTick = hitCandidateVec.front().stopTick;
711 for (
const auto& hitCandidate : hitCandidateVec) {
716 !groupedHitVec.empty()) {
717 mergedHitsVec.emplace_back(groupedHitVec);
719 groupedHitVec.clear();
723 groupedHitVec.emplace_back(hitCandidate);
725 lastStopTick = hitCandidate.stopTick;
730 if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec);
744 while ((lastItr + 1) != stopItr) {
745 if (*(lastItr + 1) > *lastItr)
break;
762 if (std::distance(startItr, minItr) > 0) {
764 while ((lastItr - 1) != startItr) {
765 if (*(lastItr - 1) < *lastItr)
break;
781 if (std::distance(startItr, lastItr) > 0) {
789 while (loopItr != startItr) {
791 if (*loopItr < 0. || !(*loopItr < *lastItr))
break;
809 if (std::distance(minItr, stopItr) > 1) {
813 while (loopItr != stopItr) {
815 if (*loopItr > 0. || !(*loopItr > *lastItr))
break;
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Utilities related to art service access.
std::vector< WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
constexpr auto abs(T v)
Returns the absolute value of the argument.
T get(std::string const &key) const
Description of geometry of one entire detector.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
This provides an interface for tools which are tasked with finding candidate hits on input waveforms...
static Globals * instance()
art framework interface to geometry description