reco_tool::PeakFitterGaussian Class Reference

Inheritance diagram for reco_tool::PeakFitterGaussian:

Public Member Functions
	PeakFitterGaussian (const fhicl::ParameterSet &pset)
void	findPeakParameters (const std::vector< float > &, const ICandidateHitFinder::HitCandidateVec &, PeakParamsVec &, double &, int &) const override

Private Attributes
const double	fMinWidth
	minimum initial width for gaussian fit More...
const double	fMaxWidthMult
	multiplier for max width for gaussian fit More...
const double	fPeakRange
	set range limits for peak center More...
const double	fAmpRange
	set range limit for peak amplitude More...
const bool	fFloatBaseline
	Allow baseline to "float" away from zero. More...
const bool	fOutputHistograms
	If true will generate summary style histograms. More...
TH1F *	fNumCandHitsHist
TH1F *	fROISizeHist
TH1F *	fCandPeakPositionHist
TH1F *	fCandPeakWidHist
TH1F *	fCandPeakAmpitudeHist
TH1F *	fCandBaselineHist
TH1F *	fFitPeakPositionHist
TH1F *	fFitPeakWidHist
TH1F *	fFitPeakAmpitudeHist
TH1F *	fFitBaselineHist
BaselinedGausFitCache	fFitCache
	Preallocated ROOT functions for the fits. More...
TH1F	fHistogram
const geo::GeometryCore *	fGeometry = lar::providerFrom<geo::Geometry>()

Additional Inherited Members
Private Types inherited from reco_tool::IPeakFitter
using	PeakParamsVec = std::vector< PeakFitParams_t >
Private Member Functions inherited from reco_tool::IPeakFitter
virtual	~IPeakFitter ()=default

Detailed Description

Definition at line 56 of file PeakFitterGaussian_tool.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 96 of file PeakFitterGaussian_tool.cc.

                                                                   :
    fMinWidth(pset.get<double>("MinWidth",         0.5)),
    fMaxWidthMult (pset.get<double>("MaxWidthMult",     3.)),
    fPeakRange(pset.get<double>("PeakRangeFact",    2.)),
    fAmpRange(pset.get<double>("PeakAmpRange",     2.)),
    fFloatBaseline(pset.get< bool >("FloatBaseline",    false)),
    fOutputHistograms(pset.get< bool >("OutputHistograms", false))
{
    fHistogram = TH1F("PeakFitterHitSignal","",500,0.,500.);

    fHistogram.Sumw2();

    std::string function = "Gaus(0)";

    // 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;

        // Make a directory for these histograms
        art::TFileDirectory dir = tfs->mkdir("PeakFit");

        fNumCandHitsHist      = dir.make<TH1F>("NumCandHits",    "# Candidate Hits", 100,   0., 100.);
        fROISizeHist          = dir.make<TH1F>("ROISize",        "ROI Size",         400,   0., 400.);
        fCandPeakPositionHist = dir.make<TH1F>("CPeakPosition",  "Peak Position",    200,   0., 400.);
        fCandPeakWidHist      = dir.make<TH1F>("CPeadWidth",     "Peak Width",       100,   0.,  25.);
        fCandPeakAmpitudeHist = dir.make<TH1F>("CPeakAmplitude", "Peak Amplitude",   100,   0., 200.);
        fCandBaselineHist     = dir.make<TH1F>("CBaseline",      "Baseline",         200, -25.,  25.);
        fFitPeakPositionHist  = dir.make<TH1F>("FPeakPosition",  "Peak Position",    200,   0., 400.);
        fFitPeakWidHist       = dir.make<TH1F>("FPeadWidth",     "Peak Width",       100,   0.,  25.);
        fFitPeakAmpitudeHist  = dir.make<TH1F>("FPeakAmplitude", "Peak Amplitude",   100,   0., 200.);
        fFitBaselineHist      = dir.make<TH1F>("FBaseline",      "Baseline",         200, -25.,  25.);
    }

    return;
}

Member Function Documentation

void reco_tool::PeakFitterGaussian::findPeakParameters ( const std::vector< float > &  roiSignalVec, const ICandidateHitFinder::HitCandidateVec &  hitCandidateVec, PeakParamsVec &  peakParamsVec, double &  chi2PerNDF, int &  NDF  ) const

overridevirtual

Implements reco_tool::IPeakFitter.

Definition at line 136 of file PeakFitterGaussian_tool.cc.

{
    // The following is a translation of the original FitGaussians function in the original
    // GausHitFinder module originally authored by Jonathan Asaadi
    //
    // *** NOTE: this algorithm assumes the reference time for input hit candidates is to
    //           the first tick of the input waveform (ie 0)
    //
    if (hitCandidateVec.empty()) return;

    // in case of a fit failure, set the chi-square to infinity
    chi2PerNDF = std::numeric_limits<double>::infinity();

    int startTime = hitCandidateVec.front().startTick;
    int endTime   = hitCandidateVec.back().stopTick;
    int roiSize   = endTime - startTime;

    // Check to see if we need a bigger histogram for fitting
    if (roiSize > fHistogram.GetNbinsX())
    {
        std::string histName = "PeakFitterHitSignal_" + std::to_string(roiSize);
        fHistogram = TH1F(histName.c_str(),"",roiSize,0.,roiSize);
        fHistogram.Sumw2();
    }

    for(int idx = 0; idx < roiSize; idx++) fHistogram.SetBinContent(idx+1,roiSignalVec[startTime+idx]);

    // Build the string to describe the fit formula
#if 0
    std::string equation = "gaus(0)";

    for(size_t idx = 1; idx < hitCandidateVec.size(); idx++) equation += "+gaus(" + std::to_string(3*idx) + ")";

    // Set the baseline if so desired
    float baseline(0.);

    if (fFloatBaseline)
    {
        baseline = roiSignalVec[startTime];

        equation += "+" + std::to_string(baseline);
    }

    // Now define the complete function to fit
    TF1 Gaus("Gaus",equation.c_str(),0,roiSize,TF1::EAddToList::kNo);
#else
    unsigned int const nGaus = hitCandidateVec.size();
    assert(fFitCache.Get(nGaus));
    TF1& Gaus = *(fFitCache.Get(nGaus));

    // Set the baseline if so desired
    float baseline(0.);

    if (fFloatBaseline)
    {
        baseline = roiSignalVec[startTime];
        Gaus.SetParameter(nGaus * 3, baseline);
        Gaus.SetParLimits( nGaus * 3, baseline - 12., baseline + 12.);
    }
    else Gaus.FixParameter(nGaus * 3, baseline); // last parameter is the baseline

#endif // 0

    if (fOutputHistograms)
    {
        fNumCandHitsHist->Fill(hitCandidateVec.size(), 1.);
        fROISizeHist->Fill(roiSize, 1.);
        fCandBaselineHist->Fill(baseline, 1.);
    }

    // ### Setting the parameters for the Gaussian Fit ###
    int parIdx{0};
    for(auto const& candidateHit : hitCandidateVec)
    {
        double const peakMean   = candidateHit.hitCenter - float(startTime);
        double const peakWidth  = candidateHit.hitSigma;
        double const amplitude  = candidateHit.hitHeight - baseline;
        double const meanLowLim = std::max(peakMean - fPeakRange * peakWidth,              0.);
        double const meanHiLim  = std::min(peakMean + fPeakRange * peakWidth, double(roiSize));

        if (fOutputHistograms)
        {
            fCandPeakPositionHist->Fill(peakMean, 1.);
            fCandPeakWidHist->Fill(peakWidth, 1.);
            fCandPeakAmpitudeHist->Fill(amplitude, 1.);
        }

        Gaus.SetParameter(  parIdx, amplitude);
        Gaus.SetParameter(1+parIdx, peakMean);
        Gaus.SetParameter(2+parIdx, peakWidth);
        Gaus.SetParLimits(  parIdx, 0.1 * amplitude,  fAmpRange * amplitude);
        Gaus.SetParLimits(1+parIdx, meanLowLim,       meanHiLim);
        Gaus.SetParLimits(2+parIdx, std::max(fMinWidth, 0.1 * peakWidth), fMaxWidthMult * peakWidth);

        parIdx += 3;
    }

    int fitResult{-1};

    try
    { fitResult = fHistogram.Fit(&Gaus,"QNWB","", 0., roiSize);}
    catch(...)
    {mf::LogWarning("GausHitFinder") << "Fitter failed finding a hit";}

    // If the fit result is not zero there was an error
    if (!fitResult)
    {
        // ##################################################
        // ### Getting the fitted parameters from the fit ###
        // ##################################################
        chi2PerNDF = (Gaus.GetChisquare() / Gaus.GetNDF());
        NDF        = Gaus.GetNDF();

        int parIdx { 0 };
        for(size_t idx = 0; idx < hitCandidateVec.size(); idx++)
        {
            PeakFitParams_t peakParams;

            peakParams.peakAmplitude      = Gaus.GetParameter(parIdx);
            peakParams.peakAmplitudeError = Gaus.GetParError( parIdx);
            peakParams.peakCenter         = Gaus.GetParameter(parIdx + 1) + float(startTime);
            peakParams.peakCenterError    = Gaus.GetParError( parIdx + 1);
            peakParams.peakSigma          = Gaus.GetParameter(parIdx + 2);
            peakParams.peakSigmaError     = Gaus.GetParError( parIdx + 2);

            if (fOutputHistograms)
            {
                fFitPeakPositionHist->Fill(peakParams.peakCenter, 1.);
                fFitPeakWidHist->Fill(peakParams.peakSigma, 1.);
                fFitPeakAmpitudeHist->Fill(peakParams.peakAmplitude, 1.);
            }

            peakParamsVec.emplace_back(peakParams);

            parIdx += 3;
        }

        if (fOutputHistograms) fFitBaselineHist->Fill(Gaus.GetParameter(3*nGaus), 1.);
    }
#if 0
    Gaus.Delete();
#endif // 0
    return;
}

Member Data Documentation

const double reco_tool::PeakFitterGaussian::fAmpRange

private

set range limit for peak amplitude

Definition at line 72 of file PeakFitterGaussian_tool.cc.

TH1F* reco_tool::PeakFitterGaussian::fCandBaselineHist

private

Definition at line 81 of file PeakFitterGaussian_tool.cc.

TH1F* reco_tool::PeakFitterGaussian::fCandPeakAmpitudeHist

private

Definition at line 80 of file PeakFitterGaussian_tool.cc.

TH1F* reco_tool::PeakFitterGaussian::fCandPeakPositionHist

private

Definition at line 78 of file PeakFitterGaussian_tool.cc.

TH1F* reco_tool::PeakFitterGaussian::fCandPeakWidHist

private

Definition at line 79 of file PeakFitterGaussian_tool.cc.

TH1F* reco_tool::PeakFitterGaussian::fFitBaselineHist

private

Definition at line 85 of file PeakFitterGaussian_tool.cc.

BaselinedGausFitCache reco_tool::PeakFitterGaussian::fFitCache

mutableprivate

Preallocated ROOT functions for the fits.

Definition at line 87 of file PeakFitterGaussian_tool.cc.

TH1F* reco_tool::PeakFitterGaussian::fFitPeakAmpitudeHist

private

Definition at line 84 of file PeakFitterGaussian_tool.cc.

TH1F* reco_tool::PeakFitterGaussian::fFitPeakPositionHist

private

Definition at line 82 of file PeakFitterGaussian_tool.cc.

TH1F* reco_tool::PeakFitterGaussian::fFitPeakWidHist

private

Definition at line 83 of file PeakFitterGaussian_tool.cc.

const bool reco_tool::PeakFitterGaussian::fFloatBaseline

private

Allow baseline to "float" away from zero.

Definition at line 73 of file PeakFitterGaussian_tool.cc.

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

private

Definition at line 91 of file PeakFitterGaussian_tool.cc.

TH1F reco_tool::PeakFitterGaussian::fHistogram

mutableprivate

Definition at line 89 of file PeakFitterGaussian_tool.cc.

const double reco_tool::PeakFitterGaussian::fMaxWidthMult

private

multiplier for max width for gaussian fit

Definition at line 70 of file PeakFitterGaussian_tool.cc.

const double reco_tool::PeakFitterGaussian::fMinWidth

private

minimum initial width for gaussian fit

Definition at line 69 of file PeakFitterGaussian_tool.cc.

TH1F* reco_tool::PeakFitterGaussian::fNumCandHitsHist

private

Definition at line 76 of file PeakFitterGaussian_tool.cc.

const bool reco_tool::PeakFitterGaussian::fOutputHistograms

private

If true will generate summary style histograms.

Definition at line 74 of file PeakFitterGaussian_tool.cc.

const double reco_tool::PeakFitterGaussian::fPeakRange

private

set range limits for peak center

Definition at line 71 of file PeakFitterGaussian_tool.cc.

TH1F* reco_tool::PeakFitterGaussian::fROISizeHist

private

Definition at line 77 of file PeakFitterGaussian_tool.cc.

---

The documentation for this class was generated from the following file:

• larreco/larreco/HitFinder/HitFinderTools/PeakFitterGaussian_tool.cc