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;
145 <<
"Cannot fill histograms when multiple threads configured, please set " 146 "fOutputHistograms to false or change number of threads to 1\n";
151 <<
"Cannot write output waveforms when multiple threads configured, please set " 152 "fOutputHistograms to false or change number of threads to 1\n";
175 dir.make<TH1F>(Form(
"DStopStart_%1zu",
fPlane),
";Delta Stop/Start;", 100, 0., 100.);
177 dir.make<TH1F>(Form(
"DMaxTMinT_%1zu",
fPlane),
";Delta Max/Min Tick;", 100, 0., 100.);
179 dir.make<TH1F>(Form(
"DMaxDMinD_%1zu",
fPlane),
";Delta Max/Min Deriv;", 200, 0., 100.);
181 dir.make<TH1F>(Form(
"MaxErosion_%1zu",
fPlane),
";Max Erosion;", 200, -50., 150.);
183 dir.make<TH1F>(Form(
"MaxDilation_%1zu",
fPlane),
";Max Dilation;", 200, -50., 150.);
185 dir.make<TH1F>(Form(
"MaxDilEroRat_%1zu",
fPlane),
";Max Dil/Ero;", 200, -1., 1.);
192 const recob::Wire::RegionsOfInterest_t::datarange_t& dataRange,
193 const size_t roiStartTick,
194 const size_t channel,
195 const size_t eventCount,
204 const Waveform& waveform = dataRange.data();
207 fWaveformTool->triangleSmooth(rawDerivativeVec, derivativeVec);
239 for (
auto& hc : hitCandidateVec) {
241 if (startCand >= 0) hc.startTick = std::max(hc.startTick,
size_t(startCand));
248 for (
auto& hitCandidate : hitCandidateVec) {
249 size_t centerIdx = hitCandidate.hitCenter;
251 hitCandidate.hitHeight = waveform.at(centerIdx);
258 size_t plane = wids[0].Plane;
259 size_t cryo = wids[0].Cryostat;
260 size_t tpc = wids[0].TPC;
261 size_t wire = wids[0].Wire;
267 Form(
"Event%04zu/c%1zuT%1zuP%1zu/Wire_%05zu", eventCount, cryo, tpc, plane, wire));
269 size_t waveformSize = waveform.size();
270 size_t waveStart = dataRange.begin_index();
273 dir.make<TProfile>(Form(
"HWfm_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
280 TProfile* derivHist =
281 dir.make<TProfile>(Form(
"HDer_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
288 TProfile* erosionHist =
289 dir.make<TProfile>(Form(
"HEro_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
296 TProfile* dilationHist =
297 dir.make<TProfile>(Form(
"HDil_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
304 TProfile* candHitHist =
305 dir.make<TProfile>(Form(
"HCan_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
312 TProfile* maxDerivHist =
313 dir.make<TProfile>(Form(
"HMax_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
320 TProfile* strtStopHist =
321 dir.make<TProfile>(Form(
"HSSS_%03zu_roiStart-%05zu",
fChannelCnt, waveStart),
330 for (
size_t idx = 0; idx < waveform.size(); idx++) {
331 waveHist->Fill(roiStartTick + idx, waveform.at(idx));
332 derivHist->Fill(roiStartTick + idx, derivativeVec.at(idx));
333 erosionHist->Fill(roiStartTick + idx, erosionVec.at(idx));
334 dilationHist->Fill(roiStartTick + idx, dilationVec.at(idx));
338 for (
const auto& hitCandidate : hitCandidateVec) {
339 candHitHist->Fill(hitCandidate.hitCenter, hitCandidate.hitHeight);
340 maxDerivHist->Fill(hitCandidate.maxTick, hitCandidate.maxDerivative);
341 maxDerivHist->Fill(hitCandidate.minTick, hitCandidate.minDerivative);
342 strtStopHist->Fill(hitCandidate.startTick, waveform.at(hitCandidate.startTick));
343 strtStopHist->Fill(hitCandidate.stopTick, waveform.at(hitCandidate.stopTick));
352 for (
const auto& hitCandidate : hitCandidateVec) {
353 fDStopStartHist->Fill(hitCandidate.stopTick - hitCandidate.startTick, 1.);
360 std::max_element(dilationVec.begin(), dilationVec.end());
362 std::max_element(erosionVec.begin(), erosionVec.end());
366 if (
std::abs(*maxDilationItr) > 0.) dilEroRat = *maxErosionItr / *maxDilationItr;
382 const size_t roiStartTick,
383 float dilationThreshold,
390 int ticksInInputWaveform = std::distance(derivStartItr, derivStopItr);
396 float maxVal = *maxItr;
399 if (maxVal < dilationThreshold)
return;
401 int maxBin = std::distance(dilationStartItr, maxItr);
408 float firstDerivValue = -1.;
412 while (firstDerItr != derivStartItr) {
414 if (*erosionItr-- < erosionCutValue) {
418 if (*firstDerItr * firstDerivValue <= 0. && firstDerivValue > 0.)
break;
421 firstDerivValue = *firstDerItr--;
425 int hitRegionStart = std::distance(derivStartItr, firstDerItr);
431 float lastDerivValue = 1.;
433 erosionItr = erosionStartItr + maxBin;
436 while (lastDerItr != derivStopItr) {
437 if (*erosionItr++ <= erosionCutValue) {
440 if (*lastDerItr * lastDerivValue <= 0. && lastDerivValue < 0.)
break;
443 lastDerivValue = *lastDerItr++;
447 int hitRegionStop = std::distance(derivStartItr, lastDerItr);
452 derivStartItr + hitRegionStart,
454 erosionStartItr + hitRegionStart,
456 dilationStartItr + hitRegionStart,
463 derivStartItr + hitRegionStop,
464 roiStartTick + hitRegionStart,
473 erosionStartItr + hitRegionStop,
475 dilationStartItr + hitRegionStop,
477 roiStartTick + hitRegionStop,
486 const size_t roiStartTick,
488 float dPeakThreshold,
497 startItr, stopItr, dTicksThreshold, dPeakThreshold, candHitParamsVec)) {
500 for (
const auto& tuple : candHitParamsVec) {
510 std::min_element(maxItr, minItr, [](
const auto&
left,
const auto&
right) {
511 return std::fabs(
left) < std::fabs(
right);
515 if (2 * std::distance(peakItr, minItr) < std::distance(maxItr, peakItr))
517 else if (2 * std::distance(maxItr, peakItr) < std::distance(peakItr, minItr))
521 size_t hitCandidateStartTick = roiStartTick + std::distance(startItr, candStartItr);
523 if (!hitCandidateVec.empty()) {
524 int deltaTicks = hitCandidateStartTick - hitCandidateVec.back().stopTick;
526 if (deltaTicks > 0) {
527 hitCandidateStartTick -= deltaTicks / 2;
528 hitCandidateVec.back().stopTick += deltaTicks / 2;
532 hitCandidate.
startTick = hitCandidateStartTick;
533 hitCandidate.
stopTick = roiStartTick + std::distance(startItr, candStopItr);
534 hitCandidate.
maxTick = roiStartTick + std::distance(startItr, maxItr);
535 hitCandidate.
minTick = roiStartTick + std::distance(startItr, minItr);
536 hitCandidate.
maxDerivative = maxItr != stopItr ? *maxItr : 0.;
537 hitCandidate.
minDerivative = minItr != stopItr ? *minItr : 0.;
538 hitCandidate.
hitCenter = roiStartTick + std::distance(startItr, peakItr) + 0.5;
543 hitCandidateVec.push_back(hitCandidate);
639 float dPeakThreshold,
643 bool foundCandidate(
false);
645 int dTicks = std::distance(startItr, stopItr);
654 MaxMinPair minMaxPair = std::minmax_element(startItr, stopItr);
660 if (std::fabs(*maxItr) > std::fabs(*minItr))
665 int deltaTicks = std::distance(maxItr, minItr);
666 float range = *maxItr - *minItr;
668 if (deltaTicks < 2)
return foundCandidate;
671 if (deltaTicks >= dTicksThreshold && range > dPeakThreshold) foundCandidate =
true;
683 startItr, candStartItr, dTicksThreshold, dPeakThreshold, candHitParamsVec);
686 candHitParamsVec.emplace_back(candStartItr, maxItr, minItr, candStopItr);
690 candStopItr, stopItr, dTicksThreshold, dPeakThreshold, candHitParamsVec);
692 return foundCandidate || prevTicks || postTicks;
696 const recob::Wire::RegionsOfInterest_t::datarange_t&,
701 if (hitCandidateVec.empty())
return;
709 size_t lastStopTick = hitCandidateVec.front().stopTick;
712 for (
const auto& hitCandidate : hitCandidateVec) {
717 !groupedHitVec.empty()) {
718 mergedHitsVec.emplace_back(groupedHitVec);
720 groupedHitVec.clear();
724 groupedHitVec.emplace_back(hitCandidate);
726 lastStopTick = hitCandidate.stopTick;
731 if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec);
745 while ((lastItr + 1) != stopItr) {
746 if (*(lastItr + 1) > *lastItr)
break;
763 if (std::distance(startItr, minItr) > 0) {
765 while ((lastItr - 1) != startItr) {
766 if (*(lastItr - 1) < *lastItr)
break;
782 if (std::distance(startItr, lastItr) > 0) {
790 while (loopItr != startItr) {
792 if (*loopItr < 0. || !(*loopItr < *lastItr))
break;
810 if (std::distance(minItr, stopItr) > 1) {
814 while (loopItr != stopItr) {
816 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)
constexpr auto abs(T v)
Returns the absolute value of the argument.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Interface for a class providing readout channel mapping to geometry.
T get(std::string const &key) const
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.
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
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