Inheritance diagram for reco_tool::CandHitDerivative:

Public Member Functions
	CandHitDerivative (const fhicl::ParameterSet &pset)

void	findHitCandidates (const recob::Wire::RegionsOfInterest_t::datarange_t &, const size_t, const size_t, const size_t, HitCandidateVec &) const override

void	MergeHitCandidates (const recob::Wire::RegionsOfInterest_t::datarange_t &, const HitCandidateVec &, MergeHitCandidateVec &) const override

Private Types
using	HitCandidateVec = std::vector< HitCandidate >

using	MergeHitCandidateVec = std::vector< HitCandidateVec >

using	Waveform = std::vector< float >

Private Member Functions
void	findHitCandidates (Waveform::const_iterator, Waveform::const_iterator, const size_t, int, float, HitCandidateVec &) const

Waveform::const_iterator	findNearestMax (Waveform::const_iterator, Waveform::const_iterator) const

Waveform::const_iterator	findNearestMin (Waveform::const_iterator, Waveform::const_iterator) const

Waveform::const_iterator	findStartTick (Waveform::const_iterator, Waveform::const_iterator) const

Waveform::const_iterator	findStopTick (Waveform::const_iterator, Waveform::const_iterator) const

Private Attributes
size_t	fPlane

int	fMinDeltaTicks

int	fMaxDeltaTicks

float	fMinDeltaPeaks

float	fMinHitHeight

size_t	fNumInterveningTicks

bool	fOutputHistograms

art::TFileDirectory *	fHistDirectory

TH1F *	fDStopStartHist

TH1F *	fDMaxTickMinTickHist

TH1F *	fDMaxDerivMinDerivHist

std::map< size_t, int >	fChannelCntMap

std::unique_ptr< reco_tool::IWaveformTool >	fWaveformTool

const geo::GeometryCore *	fGeometry = lar::providerFrom<geo::Geometry>()

Detailed Description

Definition at line 24 of file CandHitDerivative_tool.cc.

Member Typedef Documentation

using reco_tool::ICandidateHitFinder::HitCandidateVec = std::vector<HitCandidate>

inherited

Definition at line 37 of file ICandidateHitFinder.h.

using reco_tool::ICandidateHitFinder::MergeHitCandidateVec = std::vector<HitCandidateVec>

inherited

Definition at line 38 of file ICandidateHitFinder.h.

using reco_tool::ICandidateHitFinder::Waveform = std::vector<float>

inherited

Definition at line 40 of file ICandidateHitFinder.h.

Constructor & Destructor Documentation

reco_tool::CandHitDerivative::CandHitDerivative ( const fhicl::ParameterSet & pset )

explicit

Definition at line 84 of file CandHitDerivative_tool.cc.

References dir, fDMaxDerivMinDerivHist, fDMaxTickMinTickHist, fDStopStartHist, fHistDirectory, fMaxDeltaTicks, fMinDeltaPeaks, fMinDeltaTicks, fMinHitHeight, fNumInterveningTicks, fOutputHistograms, fPlane, fWaveformTool, art::ServiceHandle< T, SCOPE >::get(), and fhicl::ParameterSet::get().

   {
     fPlane = pset.get<size_t>("Plane", 0);
     fMinDeltaTicks = pset.get<int>("MinDeltaTicks", 0);
     fMaxDeltaTicks = pset.get<int>("MaxDeltaTicks", 30);
     fMinDeltaPeaks = pset.get<float>("MinDeltaPeaks", 0.025);
     fMinHitHeight = pset.get<float>("MinHitHeight", 2.0);
     fNumInterveningTicks = pset.get<size_t>("NumInterveningTicks", 6);
     fOutputHistograms = pset.get<bool>("OutputHistograms", false);
 
     // Recover the baseline tool
     fWaveformTool =
       art::make_tool<reco_tool::IWaveformTool>(pset.get<fhicl::ParameterSet>("WaveformAlgs"));
 
     // If asked, define the global histograms
     if (fOutputHistograms) {
       // Access ART's TFileService, which will handle creating and writing
       // histograms and n-tuples for us.
       art::ServiceHandle<art::TFileService> tfs;
 
       fHistDirectory = tfs.get();
 
       // Make a directory for these histograms
       art::TFileDirectory dir = fHistDirectory->mkdir(Form("HitPlane_%1zu", fPlane));
 
       fDStopStartHist =
         dir.make<TH1F>(Form("DStopStart_%1zu", fPlane), ";Delta Stop/Start;", 200, 0., 200.);
       fDMaxTickMinTickHist =
         dir.make<TH1F>(Form("DMaxTMinT_%1zu", fPlane), ";Delta Max/Min Tick;", 200, 0., 200.);
       fDMaxDerivMinDerivHist =
         dir.make<TH1F>(Form("DMaxDMinD_%1zu", fPlane), ";Delta Max/Min Deriv;", 200, 0., 200.);
     }
 
     return;
   }

Member Function Documentation

void reco_tool::CandHitDerivative::findHitCandidates	(	const recob::Wire::RegionsOfInterest_t::datarange_t &	dataRange,
		const size_t	roiStartTick,
		const size_t	channel,
		const size_t	eventCount,
		HitCandidateVec &	hitCandidateVec
	)		const

overridevirtual

Implements reco_tool::ICandidateHitFinder.

Definition at line 120 of file CandHitDerivative_tool.cc.

References geo::GeometryCore::ChannelToWire(), dir, fChannelCntMap, fDMaxDerivMinDerivHist, fDMaxTickMinTickHist, fDStopStartHist, fGeometry, fHistDirectory, fMinDeltaPeaks, fMinDeltaTicks, fOutputHistograms, and fWaveformTool.

Referenced by findHitCandidates().

   {
     // In this case we want to find hit candidates based on the derivative of of the input waveform
     // We get this from our waveform algs too...
     Waveform rawDerivativeVec;
     Waveform derivativeVec;
 
     // Recover the actual waveform
     const Waveform& waveform = dataRange.data();
 
     fWaveformTool->firstDerivative(waveform, rawDerivativeVec);
     fWaveformTool->triangleSmooth(rawDerivativeVec, derivativeVec);
 
     std::vector<geo::WireID> wids = fGeometry->ChannelToWire(channel);
     size_t plane = wids[0].Plane;
     size_t cryo = wids[0].Cryostat;
     size_t tpc = wids[0].TPC;
     size_t wire = wids[0].Wire;
 
     // Just make sure the input candidate hit vector has been cleared
     hitCandidateVec.clear();
 
     // Now find the hits
     findHitCandidates(derivativeVec.begin(),
                       derivativeVec.end(),
                       roiStartTick,
                       fMinDeltaTicks,
                       fMinDeltaPeaks,
                       hitCandidateVec);
 
     if (hitCandidateVec.empty()) {
       if (plane == 0) {
         std::cout << "** C/T/P: " << cryo << "/" << tpc << "/" << plane << ", wire: " << wire
                   << " has not hits with input size: " << waveform.size() << std::endl;
       }
     }
 
     // Reset the hit height from the input waveform
     for (auto& hitCandidate : hitCandidateVec) {
       size_t centerIdx = hitCandidate.hitCenter;
 
       hitCandidate.hitHeight = waveform.at(centerIdx);
     }
 
     // Keep track of histograms if requested
     if (fOutputHistograms) {
       // Recover the details...
       //        std::vector<geo::WireID> wids  = fGeometry->ChannelToWire(channel);
       //        size_t                   plane = wids[0].Plane;
       //        size_t                   cryo  = wids[0].Cryostat;
       //        size_t                   tpc   = wids[0].TPC;
       //        size_t                   wire  = wids[0].Wire;
 
       size_t channelCnt = fChannelCntMap[channel]++;
 
       // Make a directory for these histograms
       art::TFileDirectory dir = fHistDirectory->mkdir(
         Form("HitPlane_%1zu/ev%04zu/c%1zut%1zuwire_%05zu", plane, eventCount, cryo, tpc, wire));
 
       size_t waveformSize = waveform.size();
       int waveStart = roiStartTick;
       int waveStop = waveStart + waveformSize;
 
       TProfile* waveHist = dir.make<TProfile>(
         Form("HWfm_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
         "Waveform",
         waveformSize,
         waveStart,
         waveStop,
         -500.,
         500.);
       TProfile* derivHist = dir.make<TProfile>(
         Form("HDer_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
         "Derivative",
         waveformSize,
         waveStart,
         waveStop,
         -500.,
         500.);
       TProfile* candHitHist = dir.make<TProfile>(
         Form("HCan_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
         "Cand Hits",
         waveformSize,
         waveStart,
         waveStop,
         -500.,
         500.);
       TProfile* maxDerivHist = dir.make<TProfile>(
         Form("HMax_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
         "Maxima",
         waveformSize,
         waveStart,
         waveStop,
         -500.,
         500.);
 
       // Fill wave/derivative
       for (size_t idx = 0; idx < waveform.size(); idx++) {
         waveHist->Fill(roiStartTick + idx, waveform.at(idx));
         derivHist->Fill(roiStartTick + idx, derivativeVec.at(idx));
       }
 
       // Fill hits
       for (const auto& hitCandidate : hitCandidateVec) {
         candHitHist->Fill(hitCandidate.hitCenter, hitCandidate.hitHeight);
         maxDerivHist->Fill(hitCandidate.maxTick, hitCandidate.maxDerivative);
         maxDerivHist->Fill(hitCandidate.minTick, hitCandidate.minDerivative);
 
         fDStopStartHist->Fill(hitCandidate.stopTick - hitCandidate.startTick, 1.);
         fDMaxTickMinTickHist->Fill(hitCandidate.minTick - hitCandidate.maxTick, 1.);
         fDMaxDerivMinDerivHist->Fill(hitCandidate.maxDerivative - hitCandidate.minDerivative, 1.);
       }
     }
 
     return;
   }

void reco_tool::CandHitDerivative::findHitCandidates	(	Waveform::const_iterator	startItr,
		Waveform::const_iterator	stopItr,
		const size_t	roiStartTick,
		int	dTicksThreshold,
		float	dPeakThreshold,
		HitCandidateVec &	hitCandidateVec
	)		const

private

Definition at line 242 of file CandHitDerivative_tool.cc.

   {
     // Search for candidate hits...
     // The idea will be to find the largest deviation in the input derivative waveform as the starting point. Depending
     // on if a maximum or minimum, we search forward or backward to find the minimum or maximum that our extremum
     // corresponds to.
     std::pair<Waveform::const_iterator, Waveform::const_iterator> minMaxPair =
       std::minmax_element(startItr, stopItr);
 
     Waveform::const_iterator maxItr = minMaxPair.second;
     Waveform::const_iterator minItr = minMaxPair.first;
 
     // Use the larger of the two as the starting point and recover the nearest max or min
     if (std::fabs(*maxItr) > std::fabs(*minItr))
       minItr = findNearestMin(maxItr, stopItr);
     else
       maxItr = findNearestMax(minItr, startItr);
 
     int deltaTicks = std::distance(maxItr, minItr);
     float range = *maxItr - *minItr;
 
     // At some point small rolling oscillations on the waveform need to be ignored...
     if (deltaTicks >= dTicksThreshold && range > dPeakThreshold) {
       // Need to back up to find zero crossing, this will be the starting point of our
       // candidate hit but also the endpoint of the pre sub-waveform we'll search next
       Waveform::const_iterator newEndItr = findStartTick(maxItr, startItr);
 
       int startTick = std::distance(startItr, newEndItr);
 
       // Now need to go forward to again get close to zero, this will then be the end point
       // of our candidate hit and the starting point for the post sub-waveform to search
       Waveform::const_iterator newStartItr = findStopTick(minItr, stopItr);
 
       int stopTick = std::distance(startItr, newStartItr);
 
       // Find hits in the section of the waveform leading up to this candidate hit
       if (startTick > dTicksThreshold) {
         // Special handling for merged hits
         if (*(newEndItr - 1) > 0.) {
           dTicksThreshold = 2;
           dPeakThreshold = 0.;
         }
         else {
           dTicksThreshold = fMinDeltaTicks;
           dPeakThreshold = fMinDeltaPeaks;
         }
 
         findHitCandidates(
           startItr, newEndItr + 1, roiStartTick, dTicksThreshold, dPeakThreshold, hitCandidateVec);
       }
 
       // Create a new hit candidate and store away
       HitCandidate hitCandidate;
 
       Waveform::const_iterator peakItr =
         std::min_element(maxItr, minItr, [](const auto& left, const auto& right) {
           return std::fabs(left) < std::fabs(right);
         });
 
       // Check balance
       if (2 * std::distance(peakItr, minItr) < std::distance(maxItr, peakItr))
         peakItr--;
       else if (2 * std::distance(maxItr, peakItr) < std::distance(peakItr, minItr))
         peakItr++;
 
       hitCandidate.startTick = roiStartTick + startTick;
       hitCandidate.stopTick = roiStartTick + stopTick;
       hitCandidate.maxTick = roiStartTick + std::distance(startItr, maxItr);
       hitCandidate.minTick = roiStartTick + std::distance(startItr, minItr);
       hitCandidate.maxDerivative = maxItr != stopItr ? *maxItr : 0.;
       hitCandidate.minDerivative = minItr != stopItr ? *minItr : 0.;
       hitCandidate.hitCenter = roiStartTick + std::distance(startItr, peakItr) + 0.5;
       hitCandidate.hitSigma = 0.5 * float(hitCandidate.minTick - hitCandidate.maxTick);
       hitCandidate.hitHeight =
         hitCandidate.hitSigma * (hitCandidate.maxDerivative - hitCandidate.minDerivative) / 1.2130;
 
       hitCandidateVec.push_back(hitCandidate);
 
       // Finally, search the section of the waveform following this candidate for more hits
       if (std::distance(newStartItr, stopItr) > dTicksThreshold) {
         // Special handling for merged hits
         if (*(newStartItr + 1) < 0.) {
           dTicksThreshold = 2;
           dPeakThreshold = 0.;
         }
         else {
           dTicksThreshold = fMinDeltaTicks;
           dPeakThreshold = fMinDeltaPeaks;
         }
 
         findHitCandidates(newStartItr,
                           stopItr,
                           roiStartTick + stopTick,
                           dTicksThreshold,
                           dPeakThreshold,
                           hitCandidateVec);
       }
     }
 
     return;
   }

ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitDerivative::findNearestMax	(	Waveform::const_iterator	minItr,
		Waveform::const_iterator	startItr
	)		const

private

Definition at line 408 of file CandHitDerivative_tool.cc.

Referenced by findHitCandidates().

   {
     // Set the internal loop variable...
     Waveform::const_iterator lastItr = minItr;
 
     // One extra condition to watch for here, make sure we can actually back up!
     if (std::distance(startItr, minItr) > 0) {
       // Similar to searching for a maximum, we loop backward over ticks looking for the waveform to start decreasing
       while ((lastItr - 1) != startItr) {
         if (*(lastItr - 1) < *lastItr) break;
 
         lastItr--;
       }
     }
 
     return lastItr;
   }

ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitDerivative::findNearestMin	(	Waveform::const_iterator	maxItr,
		Waveform::const_iterator	stopItr
	)		const

private

Definition at line 389 of file CandHitDerivative_tool.cc.

Referenced by findHitCandidates().

   {
     // reset the min iterator and search forward to find the nearest minimum
     Waveform::const_iterator lastItr = maxItr;
 
     // The strategy is simple...
     // We are at a maximum so we search forward until we find the lowest negative point
     while ((lastItr + 1) != stopItr) {
       if (*(lastItr + 1) > *lastItr) break;
 
       lastItr++;
     }
 
     // The minimum will be the last iterator value...
     return lastItr;
   }

ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitDerivative::findStartTick	(	Waveform::const_iterator	maxItr,
		Waveform::const_iterator	startItr
	)		const

private

Definition at line 428 of file CandHitDerivative_tool.cc.

Referenced by findHitCandidates().

   {
     Waveform::const_iterator lastItr = maxItr;
 
     // If we can't back up then there is nothing to do
     if (std::distance(startItr, lastItr) > 0) {
       // In theory, we are starting at a maximum and want to find the "start" of the candidate peak
       // Ideally we would look to search backward to the point where the (derivative) waveform crosses zero again.
       // However, the complication is that we need to watch for the case where two peaks are merged together and
       // we might run through another peak before crossing zero...
       // So... loop until we hit the startItr...
       Waveform::const_iterator loopItr = lastItr - 1;
 
       while (loopItr != startItr) {
         // Ideal world case, we cross zero... but we might encounter a minimum... or an inflection point
         if (*loopItr < 0. || !(*loopItr < *lastItr)) break;
 
         lastItr = loopItr--;
       }
     }
     else
       lastItr = startItr;
 
     return lastItr;
   }

ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitDerivative::findStopTick	(	Waveform::const_iterator	minItr,
		Waveform::const_iterator	stopItr
	)		const

private

Definition at line 456 of file CandHitDerivative_tool.cc.

References DEFINE_ART_CLASS_TOOL.

Referenced by findHitCandidates().

   {
     Waveform::const_iterator lastItr = minItr;
 
     // If we can't go forward then there is really nothing to do
     if (std::distance(minItr, stopItr) > 1) {
       // Pretty much the same strategy as for finding the start tick...
       Waveform::const_iterator loopItr = lastItr + 1;
 
       while (loopItr != stopItr) {
         // Ideal case that we have crossed zero coming from a minimum... but watch for a maximum as well
         if (*loopItr > 0. || !(*loopItr > *lastItr)) break;
 
         lastItr = loopItr++;
       }
     }
 
     return lastItr;
   }

void reco_tool::CandHitDerivative::MergeHitCandidates	(	const recob::Wire::RegionsOfInterest_t::datarange_t &	,
		const HitCandidateVec &	hitCandidateVec,
		MergeHitCandidateVec &	mergedHitsVec
	)		const

overridevirtual

Implements reco_tool::ICandidateHitFinder.

Definition at line 349 of file CandHitDerivative_tool.cc.

References fMinHitHeight, and fNumInterveningTicks.

   {
     // If nothing on the input end then nothing to do
     if (hitCandidateVec.empty()) return;
 
     // The idea is to group hits that "touch" so they can be part of common fit, those that
     // don't "touch" are fit independently. So here we build the output vector to achieve that
     // Get a container for the hits...
     HitCandidateVec groupedHitVec;
 
     // Initialize the end of the last hit which we'll set to the first input hit's stop
     size_t lastStopTick = hitCandidateVec.front().stopTick;
 
     // Step through the input hit candidates and group them by proximity
     for (const auto& hitCandidate : hitCandidateVec) {
       // Small pulse height hits should not be considered?
       if (hitCandidate.hitHeight > fMinHitHeight) {
         // Check condition that we have a new grouping
         if (hitCandidate.startTick > lastStopTick + fNumInterveningTicks &&
             !groupedHitVec.empty()) {
           mergedHitsVec.emplace_back(groupedHitVec);
 
           groupedHitVec.clear();
         }
 
         // Add the current hit to the current group
         groupedHitVec.emplace_back(hitCandidate);
 
         lastStopTick = hitCandidate.stopTick;
       }
     }
 
     // Check end condition
     if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec);
 
     return;
   }