221 TH1::AddDirectory(kFALSE);
251 tbb::concurrent_vector<hitstruct> hitstruct_vec;
252 tbb::concurrent_vector<hitstruct> filthitstruct_vec;
266 std::function<double(double, double, double, double, int, int)> chargeFunc =
272 int ) {
return std::sqrt(2 * TMath::Pi()) * peakAmp * peakWidth / areaNorm; };
278 chargeFunc = [](
double peakMean,
285 for (
int sigPos = low; sigPos < hi; sigPos++)
286 charge += peakAmp * TMath::Gaus(sigPos, peakMean, peakWidth);
296 static_cast<std::size_t>(0),
297 wireVecHandle->size(),
298 [&](
size_t& wireIter) {
309 std::vector<geo::WireID> wids = geom->
ChannelToWire(channel);
327 static_cast<std::size_t>(0),
329 [&](
size_t& rangeIter) {
330 const auto& range = signalROI.
range(rangeIter);
342 range, 0, channel, count, hitCandidateVec);
344 range, hitCandidateVec, mergedCandidateHitVec);
350 for (
auto& mergedCands : mergedCandidateHitVec) {
351 int startT = mergedCands.front().startTick;
352 int endT = mergedCands.back().stopTick;
357 if (endT - startT < 5)
continue;
364 int nGausForFit = mergedCands.size();
369 double chi2PerNDF(0.);
381 range.data(), mergedCands, peakParamsVec, chi2PerNDF, NDF);
384 if (!(chi2PerNDF < std::numeric_limits<double>::infinity())) {
400 int nHitsThisPulse = (endT - startT) / longPulseWidth;
404 longPulseWidth = (endT - startT) / nHitsThisPulse;
407 if (nHitsThisPulse * longPulseWidth < endT - startT) nHitsThisPulse++;
409 int firstTick = startT;
410 int lastTick = std::min(firstTick + longPulseWidth, endT);
412 peakParamsVec.clear();
413 nGausForFit = nHitsThisPulse;
415 chi2PerNDF = chi2PerNDF >
fChi2NDF ? chi2PerNDF : -1.;
417 for (
int hitIdx = 0; hitIdx < nHitsThisPulse; hitIdx++) {
420 std::accumulate(range.begin() + firstTick, range.begin() + lastTick, 0.);
421 double peakSigma = (lastTick - firstTick) / 3.;
422 double peakAmp = 0.3989 * ROIsumADC / peakSigma;
423 double peakMean = (firstTick + lastTick) / 2.;
435 peakParamsVec.push_back(peakParams);
438 firstTick = lastTick;
439 lastTick = std::min(lastTick + longPulseWidth, endT);
450 std::vector<recob::Hit> filteredHitVec;
454 float nextpeakSig(0);
455 float prevpeakSig(0);
456 float nsigmaADC(2.0);
459 for (
const auto& peakParams : peakParamsVec) {
461 float peakAmp = peakParams.peakAmplitude;
462 float peakMean = peakParams.peakCenter;
463 float peakWidth = peakParams.peakSigma;
476 if (numHits < nGausForFit - 1) {
477 nextpeak = (peakParamsVec.at(numHits + 1)).peakCenter;
478 nextpeakSig = (peakParamsVec.at(numHits + 1)).peakSigma;
482 prevpeak = (peakParamsVec.at(numHits - 1)).peakCenter;
483 prevpeakSig = (peakParamsVec.at(numHits - 1)).peakSigma;
488 if (std::isnan(peakAmp)) {
489 std::cout <<
"**** hit peak amplitude is a nan! Channel: " << channel
490 <<
", start tick: " << startT << std::endl;
495 float peakAmpErr = peakParams.peakAmplitudeError;
496 float peakMeanErr = peakParams.peakCenterError;
497 float peakWidthErr = peakParams.peakSigmaError;
501 chargeFunc(peakMean, peakAmp, peakWidth,
fAreaNormsVec[plane], startT, endT);
504 std::sqrt(TMath::Pi()) * (peakAmpErr * peakWidthErr + peakWidthErr * peakAmpErr);
512 range.begin() + peakMean - nsigmaADC * peakWidth;
514 range.begin() + peakMean + nsigmaADC * peakWidth;
516 if (nGausForFit > 1) {
518 if ((peakMean - nsigmaADC * peakWidth) < (prevpeak + nsigmaADC * prevpeakSig)) {
519 float difPeak = peakMean - prevpeak;
520 float weightpeak = prevpeakSig / (prevpeakSig + peakWidth);
521 HitsumStartItr = range.begin() + prevpeak + difPeak * weightpeak;
522 newleft = prevpeak + difPeak * weightpeak;
526 if (numHits < nGausForFit - 1) {
527 if ((peakMean + nsigmaADC * peakWidth) > (nextpeak - nsigmaADC * nextpeakSig)) {
528 float difPeak = nextpeak - peakMean;
529 float weightpeak = peakWidth / (nextpeakSig + peakWidth);
530 HitsumEndItr = range.begin() + peakMean + difPeak * weightpeak;
531 newright = peakMean + difPeak * weightpeak;
537 if (newright - newleft < 0)
continue;
540 if (HitsumStartItr < sumStartItr) HitsumStartItr = sumStartItr;
542 if (HitsumEndItr > sumEndItr) HitsumEndItr = sumEndItr;
545 double ROIsumADC = std::accumulate(sumStartItr, sumEndItr, 0.);
546 double HitsumADC = std::accumulate(HitsumStartItr, HitsumEndItr, 0.);
552 startT + roiFirstBinTick,
553 endT + roiFirstBinTick,
555 peakMean + roiFirstBinTick,
569 if (filteredHitCol) filteredHitVec.push_back(hitcreator.copy());
574 hitstruct
tmp{std::move(
hit), wire};
575 hitstruct_vec.push_back(std::move(
tmp));
581 if (filteredHitCol && !filteredHitVec.empty()) {
590 std::sort(filteredHitVec.begin(),
591 filteredHitVec.end(),
592 [](
const auto&
left,
const auto&
right) {
593 return left.PeakAmplitude() >
right.PeakAmplitude();
599 filteredHitVec.clear();
602 if (filteredHitVec.size() > 1) {
604 float largestPH = filteredHitVec.front().PeakAmplitude();
610 std::find_if(filteredHitVec.begin(),
611 filteredHitVec.end(),
612 [largestPH, threshold](
const auto&
hit) {
613 return hit.PeakAmplitude() < 8. &&
614 hit.PeakAmplitude() / largestPH < threshold;
618 if (smallHitItr != filteredHitVec.end())
619 filteredHitVec.resize(std::distance(filteredHitVec.begin(), smallHitItr));
622 std::sort(filteredHitVec.begin(),
623 filteredHitVec.end(),
624 [](
const auto&
left,
const auto&
right) {
625 return left.PeakTime() <
right.PeakTime();
630 for (
const auto& filteredHit : filteredHitVec)
632 hitstruct
tmp{std::move(filteredHit), wire};
633 filthitstruct_vec.push_back(std::move(
tmp));
644 for (
size_t i = 0; i < hitstruct_vec.size(); i++) {
645 allHitCol.emplace_back(hitstruct_vec[i].hit_tbb, hitstruct_vec[i].wire_tbb);
648 for (
size_t j = 0; j < filthitstruct_vec.size(); j++) {
649 filteredHitCol->
emplace_back(filthitstruct_vec[j].hit_tbb, filthitstruct_vec[j].wire_tbb);
655 if (filteredHitCol) {
667 allHitCol.put_into(evt);
671 allHitCol.put_into(evt);
const datarange_t & range(size_t i) const
Returns the i-th non-void range (zero-based)
size_type n_ranges() const
Returns the internal list of non-void ranges.
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
std::vector< WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
std::unique_ptr< HitFilterAlg > fHitFilterAlg
algorithm used to filter out noise hits
unsigned int PlaneID_t
Type for the ID number.
const std::vector< int > fLongPulseWidthVec
Sets width of hits used to describe long pulses.
std::unique_ptr< reco_tool::IPeakFitter > fPeakFitterTool
Perform fit to candidate peaks.
const int fAreaMethod
Type of area calculation.
const std::string fCalDataModuleLabel
const std::vector< float > fPulseWidthCuts
int TDCtick_t
Type representing a TDC tick.
Class managing the creation of a new recob::Hit object.
A class handling a collection of hits and its associations.
std::atomic< size_t > fEventCount
const size_t fMaxMultiHit
maximum hits for multi fit
const std::vector< int > fLongMaxHitsVec
Maximum number hits on a really long pulse train.
void emplace_back(recob::Hit &&hit, art::Ptr< recob::Wire > const &wire=art::Ptr< recob::Wire >(), art::Ptr< raw::RawDigit > const &digits=art::Ptr< raw::RawDigit >())
Adds the specified hit to the data collection.
PlaneID_t Plane
Index of the plane within its TPC.
void put_into(art::Event &)
Moves the data into an event.
Detector simulation of raw signals on wires.
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
const std::vector< float > fPulseRatioCuts
const std::vector< float > fPulseHeightCuts
maximum Chisquared / NDF allowed for a hit to be saved
const std::string fAllHitsInstanceName
2D representation of charge deposited in the TDC/wire plane
const std::vector< double > fAreaNormsVec
factors for converting area to same units as peak height
unsigned int ChannelID_t
Type representing the ID of a readout channel.
std::vector< std::unique_ptr< reco_tool::ICandidateHitFinder > > fHitFinderToolVec
For finding candidate hits.
1.8.11 
