#ifndef _utl_QuadraticFitter_h_ #define _utl_QuadraticFitter_h_ #include #include #include #include #include namespace utl { class PoissonianErrors { public: static double GetError(const double& value) { return std::sqrt(value); } static double GetError2(const double& value) { return value; } }; class FlatErrors { public: static double GetError(const double& /*value*/) { return 1; } static double GetError2(const double& /*value*/) { return 1; } }; /** \class QuadraticFitter QuadraticFitter.h "utl/QuadraticFitter.h" \author Darko Veberic \date 27 August 2007 \version $Id: QuadraticFitter.h 14717 2009-09-17 20:24:36Z lukas $ */ template class QuadraticFitter { public: QuadraticFitter(const Histogram& h) : fStartBin(0), fStopBin(h.GetNBins()), fHistogram(h) { } QuadraticFitter(const Histogram& h, const int startBin, const int stopBin) : fStartBin(startBin), fStopBin(stopBin), fHistogram(h) { if (startBin < 0 || stopBin - startBin < 3 || stopBin >= int(h.GetNBins())) throw utl::OutOfBoundException("start, stop bin out of bounds of histogram"); } QuadraticFitter(const Histogram& h, const double xStart, const double xStop) : fHistogram(h) { int start = h.GetBinIndex(xStart); if (start == Histogram::eUnderflow) start = 0; fStartBin = start; int stop = h.GetBinIndex(xStop); if (stop == Histogram::eOverflow) stop = h.GetNBins() - 1; fStopBin = stop; if (start < 0 || stop < 0 || fStopBin - fStartBin < 3) throw utl::OutOfBoundException("start, stop bin out of bounds of histogram"); } double GetExtremePosition() { MakeSums(); MakeCoeficientTimesDeterminant(1); MakeCoeficientTimesDeterminant(2); return - fK[1] / (2 * fK[2]); } double GetExtremePosition(double& posError, int& ndof) { const double pos = GetExtremePosition(); const double sigma2k1 = Sqr(fx2Sum) - fx0Sum * fx4Sum; // /det const double sigma2k2 = Sqr(fx1Sum) - fx0Sum * fx2Sum; // /det const double covk1k2 = fx0Sum * fx3Sum - fx1Sum * fx2Sum; // /det const double k22 = Sqr(fK[2]); // /det^2 MakeDeterminant(); posError = 0.5 * sqrt(fDet*(sigma2k1 + pos*covk1k2 + sigma2k2*Sqr(fK[1])/k22) / k22); ndof = fN - 3; return pos; } double GetExtremePosition(double& posError, double& chi2, int& ndof) { const double pos = GetExtremePosition(posError, ndof); MakeCoeficientTimesDeterminant(0); chi2 = ((Sqr(fK[0]) * fx0Sum + 2*fK[0]*fK[1] * fx1Sum + (Sqr(fK[1]) + 2*fK[0]*fK[2]) * fx2Sum + 2*fK[1]*fK[2] * fx3Sum + Sqr(fK[2]) * fx4Sum) / Sqr(fDet) - 2 * (fK[0] * fx0ySum + fK[1] * fx1ySum + fK[2] * fx2ySum) / fDet + fx0y2Sum); return pos; } double GetExtremePosition(double& posError, double& chi2, int& ndof, double polynomialCoeff[3]) { const double pos = GetExtremePosition(posError, chi2, ndof); for (int i = 0; i < 3; ++i) polynomialCoeff[i] = fK[i] / fDet; return pos; } void GetFitData(QuadraticFitData& fit) { fit.fStart = fHistogram.GetBinCenter(fStartBin); fit.fStop = fHistogram.GetBinCenter(fStopBin); fit.fExtremePosition = GetExtremePosition(fit.fExtremePositionError, fit.fChi2, fit.fNdof, fit.fCoefficients); } private: void MakeSums() { fx0Sum = fx1Sum = fx2Sum = fx3Sum = fx4Sum = fx0ySum = fx1ySum = fx2ySum = fx0y2Sum = 0; for (unsigned int i = fStartBin; i < fStopBin; ++i) { const double x = fHistogram.GetBinCenter(i); const double raw = fHistogram.GetBin(i); const double invSigma2 = Sqr(fHistogram.GetBinSize(i)) / ErrorPolicy::GetError2(raw); const double y = fHistogram.GetBinAverage(i); fx0Sum += invSigma2; fx1Sum += x * invSigma2; fx0ySum += y * invSigma2; fx0y2Sum += Sqr(y) * invSigma2; const double xy = x*y; fx1ySum += xy * invSigma2; const double x2 = Sqr(x); fx2Sum += x2 * invSigma2; fx2ySum += x*xy * invSigma2; const double x3 = x2*x; fx3Sum += x3 * invSigma2; fx4Sum += x*x3 * invSigma2; } fN = fStopBin - fStartBin; } void MakeDeterminant() { fDet = Sqr(fx2Sum) * fx2Sum + fx0Sum * Sqr(fx3Sum) + Sqr(fx1Sum) * fx4Sum - fx2Sum * (2 * fx1Sum * fx3Sum + fx0Sum * fx4Sum); } void MakeCoeficientTimesDeterminant(const int i) { switch (i) { case 0: fK[0] = Sqr(fx3Sum) * fx0ySum + fx4Sum * (fx1Sum * fx1ySum - fx2Sum * fx0ySum) + Sqr(fx2Sum) * fx2ySum - fx3Sum * (fx2Sum * fx1ySum + fx1Sum * fx2ySum); // /det case 1: fK[1] = fx4Sum * (fx1Sum * fx0ySum - fx0Sum * fx1ySum) + Sqr(fx2Sum) * fx1ySum + fx0Sum * fx3Sum * fx2ySum - fx2Sum * (fx3Sum * fx0ySum + fx1Sum * fx2ySum); // /det case 2: fK[2] = Sqr(fx2Sum) * fx0ySum + fx3Sum * (fx0Sum * fx1ySum - fx1Sum * fx0ySum) + Sqr(fx1Sum) * fx2ySum - fx2Sum * (fx1Sum * fx1ySum + fx0Sum * fx2ySum); // /det } } unsigned int fStartBin; unsigned int fStopBin; const Histogram& fHistogram; double fx0Sum; double fx1Sum; double fx2Sum; double fx3Sum; double fx4Sum; double fx0ySum; double fx1ySum; double fx2ySum; double fx0y2Sum; int fN; double fDet; double fK[3]; }; // convinience makers template inline QuadraticFitter MakeQuadraticFitter(const Histogram& h, const ErrorPolicy /*ep*/) { return QuadraticFitter(h); } template inline QuadraticFitter MakeQuadraticFitter(const Histogram& h) { return QuadraticFitter(h); } template inline QuadraticFitter MakeQuadraticFitter(const Histogram& h, const Type xStart, const Type xStop, const ErrorPolicy ep) { return QuadraticFitter(h, xStart, xStop); } template inline QuadraticFitter MakeQuadraticFitter(const Histogram& h, const Type xStart, const Type xStop) { return QuadraticFitter(h, xStart, xStop); } } #endif