#ifndef _utl_ExponentialFitter_h_ #define _utl_ExponentialFitter_h_ #include #include namespace utl { /** \brief Fit exponential function \author Darko Veberic \date 01 Sep 2007 \version $Id: ExponentialFitter.h 14717 2009-09-17 20:24:36Z lukas $ */ template class ExponentialFitter { public: ExponentialFitter(const Histogram& histogram, const int startBin, const int stopBin) : fStartBin(startBin), fStopBin(stopBin), fHistogram(histogram) { if (fStartBin < 0 || fStopBin >= histogram.GetNBins()) throw OutOfBoundException("start, stop bins out of bound"); if (fStopBin - fStartBin < 2) throw OutOfBoundException("too few data points to fit exponential"); } ExponentialFitter(const Histogram& histogram, const double xStart, const double xStop) : fStartBin(histogram.GetBinIndex(xStart)), fStopBin(histogram.GetBinIndex(xStop)), fHistogram(histogram) { if (fStartBin < 0) throw OutOfBoundException("range start error"); if (fStopBin - fStartBin < 2) throw OutOfBoundException("range order mismatch"); } double GetSlope(double& slopeError, const double offset = 0) { // assumes sigma(y) = sqrt(y) without offset // z = log(y) fSum1 = 0; fSumx = 0; fSumx2 = 0; fSumxz = 0; fSumz = 0; fN = 0; for (int i = fStartBin; i < fStopBin; ++i) { const double counts = fHistogram.GetBin(i); if (counts > 0) { const double width = fHistogram.GetBinSize(i); const double y = (counts - offset) / width; if (y > 0) { const double invSigmaz2 = counts; const double z = log(y); const double x = fHistogram.GetBinCenter(i); fSum1 += invSigmaz2; const double s2x = invSigmaz2 * x; fSumx += s2x; fSumx2 += s2x * x; fSumxz += s2x * z; fSumz += invSigmaz2 * z; ++fN; } } } if (fN > 1) { fDet = fSum1*fSumx2 - Sqr(fSumx); if (fDet) { slopeError = fSum1 / fDet; return (fSum1*fSumxz - fSumx*fSumz) / fDet; } } fDet = slopeError = 0; return 0; } void GetFit(double& amplitude, double& amplitudeError, double& slope, double& slopeError, const double offset = 0) { slope = GetSlope(slopeError, offset); if (fDet) { amplitude = exp((fSumx2*fSumz - fSumx*fSumxz) / fDet); amplitudeError = amplitude * (fSumx2 / fDet); } else { amplitude = 0; amplitudeError = 0; } } void GetFit(double& amplitude, double& amplitudeError, double& slope, double& slopeError, double& chi2, int& ndof, const double offset = 0) { GetFit(amplitude, amplitudeError, slope, slopeError, offset); if (fDet) { const double linOffset = log(amplitude); chi2 = Sqr(linOffset)*fSum1 + Sqr(slope)*fSumx2 + 2*(slope*linOffset*fSumx - linOffset*fSumz - slope*fSumxz); ndof = fN - 2; } else { chi2 = 0; ndof = 0; } } void GetFit(ExponentialFitData& ef, const double offset = 0) { ef.fStart = fHistogram.GetBinCenter(fStartBin); ef.fStop = fHistogram.GetBinCenter(fStopBin); GetFit(ef.fAmplitude, ef.fAmplitudeError, ef.fSlope, ef.fSlopeError, ef.fChi2, ef.fNdof, offset); ef.fOffset = offset; } private: const int fStartBin; const int fStopBin; const Histogram& fHistogram; double fSum1; double fSumx; double fSumx2; double fSumxz; double fSumz; double fDet; int fN; }; template inline ExponentialFitter MakeExponentialFitter(const Histogram& histogram, const T& start, const T& stop) { return ExponentialFitter(histogram, start, stop); } template inline ExponentialFitter MakeExponentialFitter(const Histogram& histogram, const std::pair& startStop) { return ExponentialFitter(histogram, startStop.first, startStop.second); } } #endif