Inheritance diagram for icarus_tool::ROIFinderDifferential:

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

	~ROIFinderDifferential ()

void	configure (const fhicl::ParameterSet &pset) override

void	initializeHistograms (art::TFileDirectory &) const override

size_t	plane () const override

void	FindROIs (const Waveform &, size_t, size_t, double, CandidateROIVec &) const override

Public Member Functions inherited from icarus_tool::IROIFinder
virtual	~IROIFinder () noexcept=default

Private Member Functions
void	findROICandidates (Waveform::const_iterator startItr, Waveform::const_iterator stopItr, size_t channel, size_t roiStartTick, float roiThreshold, CandidateROIVec &roiCandidateVec) const

void	averageInputWaveform (const Waveform &, Waveform &) const

void	smoothInputWaveform (const Waveform &, Waveform &) const

Private Attributes
size_t	fPlane

float	fNumSigma
	"# sigma" rms noise for ROI threshold More...

int	fNumBinsToAve
	Controls the averaging. More...

int	fMax2MinDistance
	Maxmimum allow peak to peak distance. More...

size_t	fMaxTicksGap
	Maximum gap between ROI's before merging. More...

unsigned short	fPreROIPad
	ROI padding. More...

unsigned short	fPostROIPad
	ROI padding. More...

float	fTruncRMSMinFraction
	or at least this fraction of time bins More...

bool	fOutputHistograms
	Output histograms? More...

std::vector< float >	fAveWeightVec

float	fWeightSum

art::TFileDirectory *	fHistDirectory

TH1F *	fDiffMeanHist

TH1F *	fDiffRmsHist

TH1F *	fMinMaxDPeakHist

TH1F *	fMinMaxDistHist

TH1F *	fMinMaxSigmaHist

TH2F *	fDPeakVDistHist

icarus_signal_processing::WaveformTools < float >	fWaveformTool

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

Additional Inherited Members
Public Types inherited from icarus_tool::IROIFinder
using	Waveform = std::vector< float >

using	CandidateROI = std::pair< size_t, size_t >

using	CandidateROIVec = std::vector< CandidateROI >

Detailed Description

Definition at line 27 of file ROIFinderDifferential_tool.cc.

Constructor & Destructor Documentation

icarus_tool::ROIFinderDifferential::ROIFinderDifferential ( const fhicl::ParameterSet & pset )

explicit

Definition at line 85 of file ROIFinderDifferential_tool.cc.

 {
     configure(pset);
 }

icarus_tool::ROIFinderDifferential::~ROIFinderDifferential ( )

Definition at line 90 of file ROIFinderDifferential_tool.cc.

91 {

92 }

Member Function Documentation

void icarus_tool::ROIFinderDifferential::averageInputWaveform	(	const Waveform &	inputWaveform,
		Waveform &	outputWaveform
	)		const

private

Definition at line 409 of file ROIFinderDifferential_tool.cc.

 {
     // Vector reduction - take the 10 bin average
     float aveSum(0.);
     
     if (outputWaveform.size() != inputWaveform.size()) outputWaveform.resize(inputWaveform.size());
     
     for(size_t idx = 0; idx < inputWaveform.size(); idx++)
     {
         aveSum += fAveWeightVec.at(idx % fNumBinsToAve) * inputWaveform.at(idx);
         
         if ((idx + 1) % fNumBinsToAve == 0)
         {
             outputWaveform[idx/fNumBinsToAve] = aveSum / fWeightSum;
             
             aveSum = 0.;
         }
     }
     
     return;
 }

void icarus_tool::ROIFinderDifferential::configure ( const fhicl::ParameterSet & pset )

overridevirtual

Implements icarus_tool::IROIFinder.

Definition at line 94 of file ROIFinderDifferential_tool.cc.

 {
     // Start by recovering the parameters
     std::vector<unsigned short> zin;
     
     fPlane               = pset.get< size_t >                     ("Plane"                     );
     fNumSigma            = pset.get< float >                      ("NumSigma"                  );
     fNumBinsToAve        = pset.get< int >                        ("NumBinsToAve"              );
     fMax2MinDistance     = pset.get< int >                        ("Max2MinDistance"           );
     fMaxTicksGap         = pset.get< size_t >                     ("MaxTicksGap",            50);
     zin                  = pset.get< std::vector<unsigned short> >("ROILeadTrailPadding"       );
     fTruncRMSMinFraction = pset.get< float >                      ("TruncRMSMinFraction",   0.6);
     fOutputHistograms    = pset.get< bool  >                      ("OutputHistograms",    false);
 
     // put the ROI pad sizes into more convenient vectors
     fPreROIPad  = zin[0];
     fPostROIPad = zin[1];
     
     // The "Max2MinDistance" is input in ticks but needs to be scaled if we are averaging
 //    if (fNumBinsToAve > 1) fMax2MinDistance = std::round(fMax2MinDistance / fNumBinsToAve) + 1;
 
     // If asked, define some 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("ROIPlane_%1zu",fPlane));
         
         fDiffMeanHist    = dir.make<TH1F>("DiffMean",    ";Diff Mean;",    100, -20.,  20.);
         fDiffRmsHist     = dir.make<TH1F>("DiffRms",     ";Diff RMS;",     100,   0.,   5.);
         
         fMinMaxDPeakHist = dir.make<TH1F>("MinMaxDPeak", ";Delta Peak;",   200, 0.,  50.);
         fMinMaxDistHist  = dir.make<TH1F>("MinMaxDist",  ";Peak Sep;",      50, 0.,  50.);
         fMinMaxSigmaHist = dir.make<TH1F>("MinMaxSigma", ";Peak Sep/rms;", 200, 0., 200.);
         fDPeakVDistHist  = dir.make<TH2F>("DPeakVDist",  ";Delta Peak; Peak Sep", 200, 0, 50., 50, 0., 50.);
     }
 
     // precalculate the weight vector to use in the averaging
     fAveWeightVec.resize(fNumBinsToAve);
     
     if (fNumBinsToAve > 1)
     {
         for(int idx = 0; idx < fNumBinsToAve/2; idx++)
         {
             float weight = idx + 1;
         
             fAveWeightVec.at(idx)                     = weight;
             fAveWeightVec.at(fNumBinsToAve - idx - 1) = weight;
         }
         
         // Watch for case of fNumBinsToAve being odd
         if (fNumBinsToAve % 2 > 0) fAveWeightVec.at(fNumBinsToAve/2) = fNumBinsToAve/2 + 1;
     }
     else fAveWeightVec.at(0) = 1.;
     
     fWeightSum = std::accumulate(fAveWeightVec.begin(),fAveWeightVec.end(), 0.);
     
     return;
 }

void icarus_tool::ROIFinderDifferential::findROICandidates	(	Waveform::const_iterator	startItr,
		Waveform::const_iterator	stopItr,
		size_t	channel,
		size_t	roiStartTick,
		float	roiThreshold,
		CandidateROIVec &	roiCandidateVec
	)		const

private

Definition at line 293 of file ROIFinderDifferential_tool.cc.

 {
     // Require a minimum length
     size_t roiLength = std::distance(startItr,stopItr);
     
     if (roiLength > 0)
     {
         // The idea is to recover the two maxima from the input waveform and center our search for the ROI
         // around them.
         std::pair<Waveform::const_iterator,Waveform::const_iterator> minMaxItr = std::minmax_element(startItr,  stopItr);
         
         Waveform::const_iterator maxItr = minMaxItr.second;
         Waveform::const_iterator minItr = minMaxItr.first;
 
         // It can be that what has been returned are the lobes of two completely separate pulses which may be separated
         // by a large number of ticks. So we can't simply assume they belong together... unless they are "close" with
         // the minimum after the maximum
         if (std::distance(maxItr,minItr) < 0 || std::distance(maxItr,minItr) > fMax2MinDistance)
         {
             // Given that, we key off the larger of the lobes and do a forward or backward search to find the minimum/maximum
             // which we think is associated to the lobe we have chosen.
             // First consider that the maximum is larger than the minimum...
             if (*maxItr > std::fabs(*minItr))
             {
                 // Check distance to the minimum we found
                 if (std::distance(maxItr,stopItr) > 0)
                 {
                     // The maximum is larger so search forward from here for the true minimum
                     minItr = maxItr;
                 
                     float prevValue = *minItr++;
                 
                     while(minItr != stopItr)
                     {
                         // Look for the point where it turns back up
                         if (prevValue < 0. && *minItr > prevValue)
                         {
                             // reset to what was the actual minimum value
                             minItr -= 1;
                             break;
                         }
                     
                         prevValue = *minItr++;
                     }
                 }
             }
             else
             {
                 // Check distance to the max
                 if (std::distance(startItr,minItr) > 0)
                 {
                     // Otherwise, we are starting at the minimum and searching backwards to find the max
                     maxItr = minItr;
                 
                     float prevValue = *maxItr--;
                 
                     while(maxItr != startItr)
                     {
                         // Decreasing for two bins
                         if (prevValue > 0. && *maxItr < prevValue)
                         {
                             // reset to what was the actual minimum value
                             maxItr += 1;
                             break;
                         }
                     
                         prevValue = *maxItr--;
                     }
                 }
             }
         }
         
         // Check that the range from maximum to minimum is over threshold
         float maxMinRange    = *maxItr - *minItr;
         int   maxMinDistance = std::distance(maxItr,minItr);
         
         if (fOutputHistograms && maxMinRange > 2. * truncRMS)
         {
             fMinMaxDPeakHist->Fill(maxMinRange, 1.);
             fMinMaxDistHist->Fill(maxMinDistance, 1.);
             fMinMaxSigmaHist->Fill(maxMinRange/truncRMS, 1.);
             fDPeakVDistHist->Fill(maxMinRange, maxMinDistance, 1.);
         }
         
         if (maxMinRange > fNumSigma * truncRMS && maxMinDistance >= 0 && maxMinDistance < fMax2MinDistance)
         {
             // To complete the edges of the ROI, search both sides for the point which is essentially back to zero,
             // or in reality back into the rms level...
             Waveform::const_reverse_iterator revItr = std::find_if(std::make_reverse_iterator(maxItr), std::make_reverse_iterator(startItr), std::bind(std::less<float>(),std::placeholders::_1,truncRMS));
 
             maxItr = revItr.base();
             minItr = std::find_if(minItr,stopItr,std::bind(std::greater<float>(),std::placeholders::_1,-truncRMS));
         
             // Before saving this ROI, look for candidates preceeding this one
             // Note that preceeding snippet will reference to the current roiStartTick
             findROICandidates(startItr, maxItr, channel, roiStartTick, truncRMS, roiCandidateVec);
         
             // Save this ROI
             size_t newStartTick = std::distance(startItr,maxItr) + roiStartTick;
             size_t newStopTick  = std::distance(startItr,minItr) + roiStartTick;
         
             roiCandidateVec.push_back(CandidateROI(newStartTick, newStopTick));
         
             // Now look for candidate ROI's downstream of this one
             findROICandidates(minItr, stopItr, channel, newStopTick, truncRMS, roiCandidateVec);
         }
     }
     
     return;
 }

void icarus_tool::ROIFinderDifferential::FindROIs	(	const Waveform &	waveform,
		size_t	channel,
		size_t	cnt,
		double	rmsNoise,
		CandidateROIVec &	roiVec
	)		const

overridevirtual

Implements icarus_tool::IROIFinder.

Definition at line 159 of file ROIFinderDifferential_tool.cc.

 {
     // The idea here is to consider the input waveform - if an induction plane then it is already in differential form,
     // if a collection plane then we form the differential - then smooth and look for peaks. The technique will be to
     // consider the peak signature as a maximum followed some bins later by a mininum and those whose difference between
     // max and min is more than the threshold are kept.
     
     // First up, determine what kind of wire we have
     std::vector<geo::WireID> wids    = fGeometry->ChannelToWire(channel);
     const geo::PlaneID&      planeID = wids[0].planeID();
     geo::SigType_t           sigType = fGeometry->SignalType(planeID);
     
     // Local copy of the input waveform
     Waveform waveformDeriv(waveform.size());
     
     // If we have a collection plane then take the derivative
     if (sigType == geo::kCollection) fWaveformTool.firstDerivative(waveform, waveformDeriv);
 
     // Otherwise a straight copy since the bipolar pulses are, effectively, derivatives
     else std::copy(waveform.begin(),waveform.end(),waveformDeriv.begin());
     
     // Do the averaging
     Waveform aveWaveformDeriv;
     
 //    averageInputWaveform(waveformDeriv, aveWaveformDeriv);
     smoothInputWaveform(waveformDeriv, aveWaveformDeriv);
 //    fWaveformTool->triangleSmooth(waveform,aveWaveformDeriv);
     
     // Scheme for finding a suitable threshold
     float truncMean(0.);
     float truncRMS(0.);
     float fullRMS(0.);
     float nSig(2.5);
     int   nTrunc(0);
     int   range(0);
     
     fWaveformTool.getTruncatedMean(aveWaveformDeriv, truncMean, nTrunc, range);
     fWaveformTool.getTruncatedRMS(aveWaveformDeriv, nSig, fullRMS, truncRMS, nTrunc);
     
     // Put a floor on the value of the truncated RMS...
     float truncRMSFloor = std::max(truncRMS, float(0.25));
 
     // Now find the ROI's
     findROICandidates(aveWaveformDeriv.begin(),aveWaveformDeriv.end(),channel,0,truncRMSFloor,roiVec);
     
     TProfile* roiHist(0);
 
     if (fOutputHistograms)
     {
         fDiffMeanHist->Fill(truncMean, 1.);
         fDiffRmsHist->Fill(truncRMS, 1.);
         
         // Try to limit to the wire number (since we are already segregated by plane)
         std::vector<geo::WireID> wids = fGeometry->ChannelToWire(channel);
         size_t                   wire = wids[0].Wire;
         
         // Make a directory for these histograms
         art::TFileDirectory dir = fHistDirectory->mkdir(Form("ROIPlane_%1zu/%03zu/wire_%05zu",fPlane,cnt,wire));
         
         // We keep track of four histograms:
         try
         {
             TProfile* waveformHist    = dir.make<TProfile>(Form("Wfm_%03zu_c%05zu",cnt,wire), "Waveform",   waveform.size(),         0, waveform.size(),         -500., 500.);
             TProfile* derivativeHist  = dir.make<TProfile>(Form("Der_%03zu_c%05zu",cnt,wire), "Derivative", waveformDeriv.size(),    0, waveformDeriv.size(),    -500., 500.);
             TProfile* aveDerivHist    = dir.make<TProfile>(Form("Ave_%03zu_c%05zu",cnt,wire), "AveDeriv",   aveWaveformDeriv.size(), 0, aveWaveformDeriv.size(), -500., 500.);
             
             roiHist = dir.make<TProfile>(Form("ROI_%03zu_c%05zu",cnt,wire), "ROI", waveform.size(), 0, waveform.size(), -500., 500.);
 
             for(size_t binIdx = 0; binIdx < waveform.size(); binIdx++)
             {
                 waveformHist->Fill(binIdx, waveform.at(binIdx));
                 derivativeHist->Fill(binIdx, waveformDeriv.at(binIdx));
             }
 
             int binIdx(0);
             
             for(const auto& val : aveWaveformDeriv) aveDerivHist->Fill(binIdx++, val);
             
         } catch(...)
         {
             std::cout << "Caught exception trying to make new hists, plane,cnt,wire: " << fPlane << ", " << cnt << ", " << wire << std::endl;
         }
     }
 
     if (roiVec.empty()) return;
     
     int nMultiplier = 1; // fNumBinsToAve
     
     // pad the ROIs
     for(auto& roi : roiVec)
     {
         if (roiHist)
         {
             roiHist->Fill(int(nMultiplier * roi.first),  std::max(3.*truncRMS,1.));
             roiHist->Fill(int(nMultiplier * roi.second), std::max(3.*truncRMS,1.));
         }
         
         // low ROI end
         roi.first  = std::max(int(nMultiplier * roi.first - fPreROIPad),0);
         // high ROI end
         roi.second = std::min(nMultiplier * roi.second + fNumBinsToAve + fPostROIPad, waveform.size() - 1);
     }
     
     // merge overlapping (or touching) ROI's
     if(roiVec.size() > 1)
     {
         // temporary vector for merged ROIs
         CandidateROIVec tempRoiVec;
         
         // Loop through candidate roi's
         size_t startRoi = roiVec.front().first;
         size_t stopRoi  = startRoi;
         
         for(auto& roi : roiVec)
         {
             if (roi.first <= stopRoi + fMaxTicksGap) stopRoi = roi.second;
             else
             {
                 tempRoiVec.push_back(CandidateROI(startRoi,stopRoi));
                 
                 startRoi = roi.first;
                 stopRoi  = roi.second;
             }
         }
         
         // Make sure to get the last one
         tempRoiVec.push_back(CandidateROI(startRoi,stopRoi));
         
         roiVec = tempRoiVec;
     }
     
     return;
 }

void icarus_tool::ROIFinderDifferential::initializeHistograms ( art::TFileDirectory & histDir ) const

overridevirtual

Implements icarus_tool::IROIFinder.

Definition at line 454 of file ROIFinderDifferential_tool.cc.

 {
     // It is assumed that the input TFileDirectory has been set up to group histograms into a common
     // folder at the calling routine's level. Here we create one more level of indirection to keep
     // histograms made by this tool separate.
 /*
     std::string dirName = "ROIFinderPlane_" + std::to_string(fPlane);
     
     art::TFileDirectory dir = histDir.mkdir(dirName.c_str());
     
     auto const* detprop      = lar::providerFrom<detinfo::DetectorPropertiesService>();
     double      samplingRate = detprop->SamplingRate();
     double      numBins      = fROIFinderVec.size();
     double      maxFreq      = 500. / samplingRate;
     std::string histName     = "ROIFinderPlane_" + std::to_string(fPlane);
     
     TH1D*       hist         = dir.make<TH1D>(histName.c_str(), "ROIFinder;Frequency(MHz)", numBins, 0., maxFreq);
     
     for(int bin = 0; bin < numBins; bin++)
     {
         double freq = maxFreq * double(bin + 0.5) / double(numBins);
         
         hist->Fill(freq, fROIFinderVec.at(bin).Re());
     }
 */
     
     return;
 }

size_t icarus_tool::ROIFinderDifferential::plane ( ) const

inlineoverridevirtual

Implements icarus_tool::IROIFinder.

Definition at line 36 of file ROIFinderDifferential_tool.cc.

36 {return fPlane;}

icarus_tool::ROIFinderDifferential::fPlane

size_t fPlane

Definition: ROIFinderDifferential_tool.cc:53

void icarus_tool::ROIFinderDifferential::smoothInputWaveform	(	const Waveform &	inputWaveform,
		Waveform &	outputWaveform
	)		const

private

Definition at line 431 of file ROIFinderDifferential_tool.cc.

 {
     // Vector smoothing - take the 10 bin average
     int   halfBins = fNumBinsToAve / 2;
     
     outputWaveform.resize(inputWaveform.size());
     
     // Set the edge bins which can't be smoothed
     std::copy(inputWaveform.begin(),inputWaveform.begin()+halfBins,outputWaveform.begin());
     std::copy(inputWaveform.end()-halfBins,inputWaveform.end(),outputWaveform.end()-halfBins);
     
     for(size_t idx = halfBins; idx < inputWaveform.size() - halfBins; idx++)
     {
         float weightedSum(0.);
         
         for(int wIdx = 0; wIdx < fNumBinsToAve; wIdx++) weightedSum += fAveWeightVec.at(wIdx) * inputWaveform.at(idx - wIdx + halfBins);
         
         outputWaveform.at(idx) = weightedSum / fWeightSum;
     }
     
     return;
 }