13 #include "cetlib_except/exception.h" 132 if (fOutputHistograms)
163 fWaveformTool->triangleSmooth(rawDerivativeVec, derivativeVec);
178 erosionVec.begin(), erosionVec.end(),
179 dilationVec.begin(), dilationVec.end(),
185 for(
auto& hitCandidate : hitCandidateVec)
187 size_t centerIdx = hitCandidate.hitCenter;
189 hitCandidate.hitHeight = waveform.at(centerIdx);
197 size_t plane = wids[0].Plane;
198 size_t cryo = wids[0].Cryostat;
199 size_t tpc = wids[0].TPC;
200 size_t wire = wids[0].Wire;
207 size_t waveformSize = waveform.size();
208 int waveStart = roiStartTick;
209 int waveStop = waveStart + waveformSize;
211 TProfile* waveHist = dir.
make<TProfile>(Form(
"HWfm_%03zu_ctw%01zu-%01zu-%01zu-%05zu",
fChannelCnt,cryo,tpc,plane,wire),
"Waveform", waveformSize, waveStart, waveStop, -500., 500.);
212 TProfile* derivHist = dir.
make<TProfile>(Form(
"HDer_%03zu_ctw%01zu-%01zu-%01zu-%05zu",
fChannelCnt,cryo,tpc,plane,wire),
"Derivative", waveformSize, waveStart, waveStop, -500., 500.);
213 TProfile* erosionHist = dir.
make<TProfile>(Form(
"HEro_%03zu_ctw%01zu-%01zu-%01zu-%05zu",
fChannelCnt,cryo,tpc,plane,wire),
"Erosion", waveformSize, waveStart, waveStop, -500., 500.);
214 TProfile* dilationHist = dir.
make<TProfile>(Form(
"HDil_%03zu_ctw%01zu-%01zu-%01zu-%05zu",
fChannelCnt,cryo,tpc,plane,wire),
"Dilation", waveformSize, waveStart, waveStop, -500., 500.);
215 TProfile* candHitHist = dir.
make<TProfile>(Form(
"HCan_%03zu_ctw%01zu-%01zu-%01zu-%05zu",
fChannelCnt,cryo,tpc,plane,wire),
"Cand Hits", waveformSize, waveStart, waveStop, -500., 500.);
216 TProfile* maxDerivHist = dir.
make<TProfile>(Form(
"HMax_%03zu_ctw%01zu-%01zu-%01zu-%05zu",
fChannelCnt,cryo,tpc,plane,wire),
"Maxima", waveformSize, waveStart, waveStop, -500., 500.);
217 TProfile* strtStopHist = dir.
make<TProfile>(Form(
"HSSS_%03zu_ctw%01zu-%01zu-%01zu-%05zu",
fChannelCnt,cryo,tpc,plane,wire),
"Start/Stop", waveformSize, waveStart, waveStop, -500., 500.);
220 for(
size_t idx = 0; idx < waveform.size(); idx++)
222 waveHist->Fill(roiStartTick + idx, waveform.at(idx));
223 derivHist->Fill(roiStartTick + idx, derivativeVec.at(idx));
224 erosionHist->Fill(roiStartTick + idx, erosionVec.at(idx));
225 dilationHist->Fill(roiStartTick + idx, dilationVec.at(idx));
229 for(
const auto& hitCandidate : hitCandidateVec)
231 candHitHist->Fill(hitCandidate.hitCenter, hitCandidate.hitHeight);
232 maxDerivHist->Fill(hitCandidate.maxTick, hitCandidate.maxDerivative);
233 maxDerivHist->Fill(hitCandidate.minTick, hitCandidate.minDerivative);
234 strtStopHist->Fill(hitCandidate.startTick, waveform.at(hitCandidate.startTick));
235 strtStopHist->Fill(hitCandidate.stopTick, waveform.at(hitCandidate.stopTick));
237 fDStopStartHist->Fill(hitCandidate.stopTick - hitCandidate.startTick, 1.);
253 float dilationThreshold,
261 float maxVal = *maxItr;
264 if (maxVal < dilationThreshold)
return;
266 int maxBin = std::distance(dilationStartItr,maxItr);
273 float firstDilationValue = *firstDilItr;
278 while(firstDilItr != dilationStartItr)
281 if (*erosionItr-- < erosionCutValue && *firstDilItr < dilationCutValue && *firstDilItr >= firstDilationValue)
break;
283 firstDilationValue = *firstDilItr--;
287 int hitRegionStart = std::distance(dilationStartItr,firstDilItr);
293 float lastDilationValue = *lastDilItr;
295 erosionItr = erosionStartItr + maxBin;
298 while(++lastDilItr != dilationStopItr)
300 if (*++erosionItr <= erosionCutValue && *lastDilItr < dilationCutValue && *lastDilItr >= lastDilationValue)
break;
302 lastDilationValue = *lastDilItr;
306 int hitRegionStop = std::distance(dilationStartItr,lastDilItr);
311 erosionStartItr, erosionStartItr + hitRegionStart,
312 dilationStartItr, dilationStartItr + hitRegionStart,
319 derivStartItr + hitRegionStop,
320 roiStartTick + hitRegionStart,
328 erosionStartItr + hitRegionStop, erosionStopItr,
329 dilationStartItr + hitRegionStop, dilationStopItr,
330 roiStartTick + hitRegionStop,
341 float dPeakThreshold,
348 std::pair<Waveform::const_iterator, Waveform::const_iterator> minMaxPair = std::minmax_element(startItr, stopItr);
354 if (std::fabs(*maxItr) > std::fabs(*minItr)) minItr =
findNearestMin(maxItr, stopItr);
357 int deltaTicks = std::distance(maxItr,minItr);
358 float range = *maxItr - *minItr;
361 if (deltaTicks >= dTicksThreshold && range > dPeakThreshold)
367 int startTick = std::distance(startItr,newEndItr);
373 int stopTick = std::distance(startItr,newStartItr);
379 if (*(newEndItr-1) > 0.) {dTicksThreshold = 2; dPeakThreshold = 0.; }
382 findHitCandidates(startItr,newEndItr+1,roiStartTick,dTicksThreshold,dPeakThreshold,hitCandidateVec);
391 if (2 * std::distance(peakItr,minItr) < std::distance(maxItr,peakItr)) peakItr--;
392 else if (2 * std::distance(maxItr,peakItr) < std::distance(peakItr,minItr)) peakItr++;
395 size_t hitCandidateStartTick = roiStartTick + startTick;
397 if (!hitCandidateVec.empty())
399 int deltaTicks = hitCandidateStartTick - hitCandidateVec.back().stopTick;
403 hitCandidateStartTick -= deltaTicks / 2;
404 hitCandidateVec.back().stopTick += deltaTicks / 2;
408 hitCandidate.startTick = hitCandidateStartTick;
409 hitCandidate.stopTick = roiStartTick +
stopTick;
410 hitCandidate.maxTick = roiStartTick + std::distance(startItr,maxItr);
411 hitCandidate.minTick = roiStartTick + std::distance(startItr,minItr);
412 hitCandidate.maxDerivative = maxItr != stopItr ? *maxItr : 0.;
413 hitCandidate.minDerivative = minItr != stopItr ? *minItr : 0.;
414 hitCandidate.hitCenter = roiStartTick + std::distance(startItr,peakItr) + 0.5;
415 hitCandidate.hitSigma = 0.5 * float(hitCandidate.minTick - hitCandidate.maxTick);
416 hitCandidate.hitHeight = hitCandidate.hitSigma * (hitCandidate.maxDerivative - hitCandidate.minDerivative) / 1.2130;
418 hitCandidateVec.push_back(hitCandidate);
421 if (std::distance(newStartItr,stopItr) > 2)
424 if (*(newStartItr+1) < 0.) {dTicksThreshold = 2; dPeakThreshold = 0.; }
427 findHitCandidates(newStartItr,stopItr,roiStartTick + stopTick,dTicksThreshold,dPeakThreshold,hitCandidateVec);
439 if (hitCandidateVec.empty())
return;
447 size_t lastStopTick = hitCandidateVec.front().stopTick;
450 for(
const auto& hitCandidate : hitCandidateVec)
458 mergedHitsVec.emplace_back(groupedHitVec);
460 groupedHitVec.clear();
464 groupedHitVec.emplace_back(hitCandidate);
466 lastStopTick = hitCandidate.stopTick;
471 if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec);
483 while((lastItr + 1) != stopItr)
485 if (*(lastItr + 1) > *lastItr)
break;
500 if (std::distance(startItr,minItr) > 0)
503 while((lastItr - 1) != startItr)
505 if (*(lastItr - 1) < *lastItr)
break;
519 if (std::distance(startItr,lastItr) > 0)
528 while(loopItr != startItr)
531 if (*loopItr < 0. || !(*loopItr < *lastItr))
break;
536 else lastItr = startItr;
546 if (std::distance(minItr,stopItr) > 1)
551 while(loopItr != stopItr)
554 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)
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
T get(std::string const &key) const
Description of geometry of one entire detector.
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
T * make(ARGS...args) const
This provides an interface for tools which are tasked with finding candidate hits on input waveforms...
art framework interface to geometry description