14 #include "art_root_io/TFileService.h" 15 #include "cetlib_except/exception.h" 99 if (fOutputHistograms) {
110 dir.make<TH1F>(Form(
"DStopStart_%1zu",
fPlane),
";Delta Stop/Start;", 200, 0., 200.);
112 dir.make<TH1F>(Form(
"DMaxTMinT_%1zu",
fPlane),
";Delta Max/Min Tick;", 200, 0., 200.);
114 dir.make<TH1F>(Form(
"DMaxDMinD_%1zu",
fPlane),
";Delta Max/Min Deriv;", 200, 0., 200.);
121 const recob::Wire::RegionsOfInterest_t::datarange_t& dataRange,
122 const size_t roiStartTick,
123 const size_t channel,
124 const size_t eventCount,
133 const Waveform& waveform = dataRange.data();
136 fWaveformTool->triangleSmooth(rawDerivativeVec, derivativeVec);
139 size_t plane = wids[0].Plane;
140 size_t cryo = wids[0].Cryostat;
141 size_t tpc = wids[0].TPC;
142 size_t wire = wids[0].Wire;
145 hitCandidateVec.clear();
155 if (hitCandidateVec.empty()) {
157 std::cout <<
"** C/T/P: " << cryo <<
"/" << tpc <<
"/" << plane <<
", wire: " << wire
158 <<
" has not hits with input size: " << waveform.size() << std::endl;
163 for (
auto& hitCandidate : hitCandidateVec) {
164 size_t centerIdx = hitCandidate.hitCenter;
166 hitCandidate.hitHeight = waveform.at(centerIdx);
182 Form(
"HitPlane_%1zu/ev%04zu/c%1zut%1zuwire_%05zu", plane, eventCount, cryo, tpc, wire));
184 size_t waveformSize = waveform.size();
185 int waveStart = roiStartTick;
186 int waveStop = waveStart + waveformSize;
188 TProfile* waveHist = dir.make<TProfile>(
189 Form(
"HWfm_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
196 TProfile* derivHist = dir.make<TProfile>(
197 Form(
"HDer_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
204 TProfile* candHitHist = dir.make<TProfile>(
205 Form(
"HCan_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
212 TProfile* maxDerivHist = dir.make<TProfile>(
213 Form(
"HMax_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
222 for (
size_t idx = 0; idx < waveform.size(); idx++) {
223 waveHist->Fill(roiStartTick + idx, waveform.at(idx));
224 derivHist->Fill(roiStartTick + idx, derivativeVec.at(idx));
228 for (
const auto& hitCandidate : hitCandidateVec) {
229 candHitHist->Fill(hitCandidate.hitCenter, hitCandidate.hitHeight);
230 maxDerivHist->Fill(hitCandidate.maxTick, hitCandidate.maxDerivative);
231 maxDerivHist->Fill(hitCandidate.minTick, hitCandidate.minDerivative);
233 fDStopStartHist->Fill(hitCandidate.stopTick - hitCandidate.startTick, 1.);
244 const size_t roiStartTick,
246 float dPeakThreshold,
253 std::pair<Waveform::const_iterator, Waveform::const_iterator> minMaxPair =
254 std::minmax_element(startItr, stopItr);
260 if (std::fabs(*maxItr) > std::fabs(*minItr))
265 int deltaTicks = std::distance(maxItr, minItr);
266 float range = *maxItr - *minItr;
269 if (deltaTicks >= dTicksThreshold && range > dPeakThreshold) {
274 int startTick = std::distance(startItr, newEndItr);
280 int stopTick = std::distance(startItr, newStartItr);
283 if (startTick > dTicksThreshold) {
285 if (*(newEndItr - 1) > 0.) {
295 startItr, newEndItr + 1, roiStartTick, dTicksThreshold, dPeakThreshold, hitCandidateVec);
302 std::min_element(maxItr, minItr, [](
const auto&
left,
const auto&
right) {
303 return std::fabs(
left) < std::fabs(
right);
307 if (2 * std::distance(peakItr, minItr) < std::distance(maxItr, peakItr))
309 else if (2 * std::distance(maxItr, peakItr) < std::distance(peakItr, minItr))
312 hitCandidate.
startTick = roiStartTick + startTick;
313 hitCandidate.
stopTick = roiStartTick + stopTick;
314 hitCandidate.
maxTick = roiStartTick + std::distance(startItr, maxItr);
315 hitCandidate.
minTick = roiStartTick + std::distance(startItr, minItr);
316 hitCandidate.
maxDerivative = maxItr != stopItr ? *maxItr : 0.;
317 hitCandidate.
minDerivative = minItr != stopItr ? *minItr : 0.;
318 hitCandidate.
hitCenter = roiStartTick + std::distance(startItr, peakItr) + 0.5;
323 hitCandidateVec.push_back(hitCandidate);
326 if (std::distance(newStartItr, stopItr) > dTicksThreshold) {
328 if (*(newStartItr + 1) < 0.) {
339 roiStartTick + stopTick,
354 if (hitCandidateVec.empty())
return;
362 size_t lastStopTick = hitCandidateVec.front().stopTick;
365 for (
const auto& hitCandidate : hitCandidateVec) {
370 !groupedHitVec.empty()) {
371 mergedHitsVec.emplace_back(groupedHitVec);
373 groupedHitVec.clear();
377 groupedHitVec.emplace_back(hitCandidate);
379 lastStopTick = hitCandidate.stopTick;
384 if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec);
398 while ((lastItr + 1) != stopItr) {
399 if (*(lastItr + 1) > *lastItr)
break;
416 if (std::distance(startItr, minItr) > 0) {
418 while ((lastItr - 1) != startItr) {
419 if (*(lastItr - 1) < *lastItr)
break;
435 if (std::distance(startItr, lastItr) > 0) {
443 while (loopItr != startItr) {
445 if (*loopItr < 0. || !(*loopItr < *lastItr))
break;
463 if (std::distance(minItr, stopItr) > 1) {
467 while (loopItr != stopItr) {
469 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.
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)
This provides an interface for tools which are tasked with finding candidate hits on input waveforms...
art framework interface to geometry description