From d6ad5e77598637f67cf23f385d5618e0a52c63cf Mon Sep 17 00:00:00 2001 From: Nima Saed-Samii Date: Tue, 3 Apr 2018 22:02:37 +0200 Subject: [PATCH] Fix TheuerkaufFiter which was broken in 8b5330c880a82b269d343b618872bf82dda3dd13 --- hdtv/rootext/fit/EEFitter.hh | 2 +- hdtv/rootext/fit/Fitter.cc | 8 ++++---- hdtv/rootext/fit/Fitter.hh | 2 +- hdtv/rootext/fit/Param.cc | 16 +++++++-------- hdtv/rootext/fit/Param.hh | 8 +++++++- hdtv/rootext/fit/TheuerkaufFitter.cc | 29 +++++++++++++--------------- 6 files changed, 34 insertions(+), 31 deletions(-) diff --git a/hdtv/rootext/fit/EEFitter.hh b/hdtv/rootext/fit/EEFitter.hh index ca949d7..4c052b4 100644 --- a/hdtv/rootext/fit/EEFitter.hh +++ b/hdtv/rootext/fit/EEFitter.hh @@ -43,10 +43,10 @@ class EEPeak { friend class EEFitter; public: + EEPeak() = default; EEPeak(const Param &pos, const Param &, const Param &sigma1, const Param &sigma2, const Param &eta, const Param &gamma); EEPeak(const EEPeak &src); - EEPeak() = default; EEPeak &operator=(const EEPeak &src); double Eval(const double *x, const double *p) const; diff --git a/hdtv/rootext/fit/Fitter.cc b/hdtv/rootext/fit/Fitter.cc index f3d8e0f..d6514a3 100644 --- a/hdtv/rootext/fit/Fitter.cc +++ b/hdtv/rootext/fit/Fitter.cc @@ -29,10 +29,10 @@ namespace HDTV { namespace Fit { -Fitter::Fitter(double r1, double r2) - : fNumParams(0), fFinal(false), fMin(std::min(r1, r2)), - fMax(std::max(r1, r2)), fNumPeaks(0), fIntBgDeg(-1), - fChisquare(std::numeric_limits::quiet_NaN()) {} +Fitter::Fitter(double r1, double r2) noexcept + : fNumParams{0}, fFinal{false}, fMin{std::min(r1, r2)}, + fMax{std::max(r1, r2)}, fNumPeaks{0}, fIntBgDeg{-1}, + fChisquare{std::numeric_limits::quiet_NaN()} {} Param Fitter::AllocParam() { return Param::Free(fNumParams++); } diff --git a/hdtv/rootext/fit/Fitter.hh b/hdtv/rootext/fit/Fitter.hh index 36ce554..6ef1f15 100644 --- a/hdtv/rootext/fit/Fitter.hh +++ b/hdtv/rootext/fit/Fitter.hh @@ -34,7 +34,7 @@ namespace Fit { //! Common base class for all different (foreground) fitters class Fitter { public: - Fitter(double r1, double r2); + Fitter(double r1, double r2) noexcept; Param AllocParam(); Param AllocParam(double ival); diff --git a/hdtv/rootext/fit/Param.cc b/hdtv/rootext/fit/Param.cc index e42848c..c49b02e 100644 --- a/hdtv/rootext/fit/Param.cc +++ b/hdtv/rootext/fit/Param.cc @@ -22,6 +22,7 @@ #include "Param.hh" +#include #include #include @@ -29,14 +30,6 @@ namespace HDTV { namespace Fit { -Param::Param(int id, double value, bool free, bool hasIVal, bool valid) { - fId = id; - fValue = value; - fFree = free; - fHasIVal = hasIVal; - fValid = valid; -} - double Param::Value(TF1 *func) const { if (fFree) { return func ? func->GetParameter(fId) @@ -56,5 +49,12 @@ double Param::Error(TF1 *func) const { } } +std::ostream &operator <<(std::ostream &lhs, const Param &rhs) { + lhs << "[Id=" << rhs._Id() << ", Free=" << rhs.IsFree() + << ", IVal=" << rhs.HasIVal() << ", Valid=" << static_cast(rhs) + << ", Value=" << rhs._Value() << ']'; + return lhs; +} + } // end namespace Fit } // end namespace HDTV diff --git a/hdtv/rootext/fit/Param.hh b/hdtv/rootext/fit/Param.hh index 4d29a36..8263e8b 100644 --- a/hdtv/rootext/fit/Param.hh +++ b/hdtv/rootext/fit/Param.hh @@ -23,6 +23,8 @@ #ifndef __Param_h__ #define __Param_h__ +#include + class TF1; namespace HDTV { @@ -53,12 +55,16 @@ public: double _Value() const { return fValue; } private: - Param(int id, double value, bool free, bool hasIVal, bool valid); + Param(int id, double value, bool free, bool hasIVal, bool valid) noexcept + : fFree{free}, fHasIVal{hasIVal}, fValid{valid}, fId{id}, fValue{value} {} + bool fFree, fHasIVal, fValid; int fId; double fValue; }; +std::ostream &operator <<(std::ostream &lhs, const Param &rhs); + } // end namespace Fit } // end namespace HDTV diff --git a/hdtv/rootext/fit/TheuerkaufFitter.cc b/hdtv/rootext/fit/TheuerkaufFitter.cc index 6bb9112..575d8ae 100644 --- a/hdtv/rootext/fit/TheuerkaufFitter.cc +++ b/hdtv/rootext/fit/TheuerkaufFitter.cc @@ -38,6 +38,7 @@ namespace HDTV { namespace Fit { // *** TheuerkaufPeak *** + //! Constructor //! Note that no tails correspond to tail parameters tl = tr = \infty. However, //! all member functions are supposed to check fHasLeftTail and fHasRightTail @@ -90,12 +91,11 @@ TheuerkaufPeak &TheuerkaufPeak::operator=(const TheuerkaufPeak &src) { return *this; } +//! Restores parameters and error for fit function. +//! Warnings: Restore function of corresponding fitter has to be called +//! beforehand! void TheuerkaufPeak::RestoreParam(const Param ¶m, double value, double error) { - //! Restores parameters and error for fit function. - //! Warnings: Restore function of corresponding fitter has to be called - //! beforehand! - if (fFunc) { fFunc->SetParameter(param._Id(), value); fFunc->SetParError(param._Id(), error); @@ -230,8 +230,8 @@ double TheuerkaufFitter::Eval(const double *x, const double *p) const { // Evaluate internal background sum += std::accumulate( - std::reverse_iterator(p + fNumParams - 1), - std::reverse_iterator(p + fNumParams - fIntBgDeg), + std::reverse_iterator(p + fNumParams), + std::reverse_iterator(p + fNumParams - fIntBgDeg - 1), 0.0, [x = *x](double bg, double param) { return bg * x + param; }); // Evaluate peaks @@ -249,8 +249,8 @@ double TheuerkaufFitter::EvalBg(const double *x, const double *p) const { // Evaluate internal background sum += std::accumulate( - std::reverse_iterator(p + fNumParams - 1), - std::reverse_iterator(p + fNumParams - fIntBgDeg), + std::reverse_iterator(p + fNumParams), + std::reverse_iterator(p + fNumParams - fIntBgDeg - 1), 0.0, [x = *x](double bg, double param) { return bg * x + param; }); // Evaluate steps in peaks @@ -296,9 +296,8 @@ TF1 *TheuerkaufFitter::GetBgFunc() { return fBgFunc.get(); } +//! Do the fit, using the given background function void TheuerkaufFitter::Fit(TH1 &hist, const Background &bg) { - //! Do the fit, using the given background function - // Refuse to fit twice if (IsFinal()) { return; @@ -309,10 +308,9 @@ void TheuerkaufFitter::Fit(TH1 &hist, const Background &bg) { _Fit(hist); } +//! Do the fit, fitting a polynomial of degree intBgDeg for the background at +//! the same time. Set intBgDeg to -1 to disable background completely. void TheuerkaufFitter::Fit(TH1 &hist, int intBgDeg) { - //! Do the fit, fitting a polynomial of degree intBgDeg for the background - //! at the same time. Set intBgDeg to -1 to disable background completely. - // Refuse to fit twice if (IsFinal()) { return; @@ -323,9 +321,8 @@ void TheuerkaufFitter::Fit(TH1 &hist, int intBgDeg) { _Fit(hist); } +//! Private: worker function to actually do the fit void TheuerkaufFitter::_Fit(TH1 &hist) { - //! Private: worker function to actually do the fit - // Allocate additional parameters for internal polynomial background // Note that a polynomial of degree n has n+1 parameters! if (fIntBgDeg >= 0) { @@ -526,7 +523,7 @@ void TheuerkaufFitter::_Fit(TH1 &hist) { double sumFreeVol = sumVol; auto ampIter = amps.begin(); for (auto &peak : fPeaks) { - if (peak.fVol.IsFree()) { + if (!peak.fVol.IsFree()) { sumFreeAmp -= *(ampIter++); sumFreeVol -= peak.fVol._Value(); } -- 2.24.1