Commit d6ad5e77 authored by Nima Saed-Samii's avatar Nima Saed-Samii

Fix TheuerkaufFiter which was broken in 8b5330c8

parent 7b5f1468
...@@ -43,10 +43,10 @@ class EEPeak { ...@@ -43,10 +43,10 @@ class EEPeak {
friend class EEFitter; friend class EEFitter;
public: public:
EEPeak() = default;
EEPeak(const Param &pos, const Param &amp, const Param &sigma1, EEPeak(const Param &pos, const Param &amp, const Param &sigma1,
const Param &sigma2, const Param &eta, const Param &gamma); const Param &sigma2, const Param &eta, const Param &gamma);
EEPeak(const EEPeak &src); EEPeak(const EEPeak &src);
EEPeak() = default;
EEPeak &operator=(const EEPeak &src); EEPeak &operator=(const EEPeak &src);
double Eval(const double *x, const double *p) const; double Eval(const double *x, const double *p) const;
......
...@@ -29,10 +29,10 @@ ...@@ -29,10 +29,10 @@
namespace HDTV { namespace HDTV {
namespace Fit { namespace Fit {
Fitter::Fitter(double r1, double r2) Fitter::Fitter(double r1, double r2) noexcept
: fNumParams(0), fFinal(false), fMin(std::min(r1, r2)), : fNumParams{0}, fFinal{false}, fMin{std::min(r1, r2)},
fMax(std::max(r1, r2)), fNumPeaks(0), fIntBgDeg(-1), fMax{std::max(r1, r2)}, fNumPeaks{0}, fIntBgDeg{-1},
fChisquare(std::numeric_limits<double>::quiet_NaN()) {} fChisquare{std::numeric_limits<double>::quiet_NaN()} {}
Param Fitter::AllocParam() { return Param::Free(fNumParams++); } Param Fitter::AllocParam() { return Param::Free(fNumParams++); }
......
...@@ -34,7 +34,7 @@ namespace Fit { ...@@ -34,7 +34,7 @@ namespace Fit {
//! Common base class for all different (foreground) fitters //! Common base class for all different (foreground) fitters
class Fitter { class Fitter {
public: public:
Fitter(double r1, double r2); Fitter(double r1, double r2) noexcept;
Param AllocParam(); Param AllocParam();
Param AllocParam(double ival); Param AllocParam(double ival);
......
...@@ -22,6 +22,7 @@ ...@@ -22,6 +22,7 @@
#include "Param.hh" #include "Param.hh"
#include <iostream>
#include <limits> #include <limits>
#include <TF1.h> #include <TF1.h>
...@@ -29,14 +30,6 @@ ...@@ -29,14 +30,6 @@
namespace HDTV { namespace HDTV {
namespace Fit { 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 { double Param::Value(TF1 *func) const {
if (fFree) { if (fFree) {
return func ? func->GetParameter(fId) return func ? func->GetParameter(fId)
...@@ -56,5 +49,12 @@ double Param::Error(TF1 *func) const { ...@@ -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<bool>(rhs)
<< ", Value=" << rhs._Value() << ']';
return lhs;
}
} // end namespace Fit } // end namespace Fit
} // end namespace HDTV } // end namespace HDTV
...@@ -23,6 +23,8 @@ ...@@ -23,6 +23,8 @@
#ifndef __Param_h__ #ifndef __Param_h__
#define __Param_h__ #define __Param_h__
#include <iosfwd>
class TF1; class TF1;
namespace HDTV { namespace HDTV {
...@@ -53,12 +55,16 @@ public: ...@@ -53,12 +55,16 @@ public:
double _Value() const { return fValue; } double _Value() const { return fValue; }
private: 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; bool fFree, fHasIVal, fValid;
int fId; int fId;
double fValue; double fValue;
}; };
std::ostream &operator <<(std::ostream &lhs, const Param &rhs);
} // end namespace Fit } // end namespace Fit
} // end namespace HDTV } // end namespace HDTV
......
...@@ -38,6 +38,7 @@ namespace HDTV { ...@@ -38,6 +38,7 @@ namespace HDTV {
namespace Fit { namespace Fit {
// *** TheuerkaufPeak *** // *** TheuerkaufPeak ***
//! Constructor //! Constructor
//! Note that no tails correspond to tail parameters tl = tr = \infty. However, //! Note that no tails correspond to tail parameters tl = tr = \infty. However,
//! all member functions are supposed to check fHasLeftTail and fHasRightTail //! all member functions are supposed to check fHasLeftTail and fHasRightTail
...@@ -90,12 +91,11 @@ TheuerkaufPeak &TheuerkaufPeak::operator=(const TheuerkaufPeak &src) { ...@@ -90,12 +91,11 @@ TheuerkaufPeak &TheuerkaufPeak::operator=(const TheuerkaufPeak &src) {
return *this; 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 &param, double value, void TheuerkaufPeak::RestoreParam(const Param &param, double value,
double error) { double error) {
//! Restores parameters and error for fit function.
//! Warnings: Restore function of corresponding fitter has to be called
//! beforehand!
if (fFunc) { if (fFunc) {
fFunc->SetParameter(param._Id(), value); fFunc->SetParameter(param._Id(), value);
fFunc->SetParError(param._Id(), error); fFunc->SetParError(param._Id(), error);
...@@ -230,8 +230,8 @@ double TheuerkaufFitter::Eval(const double *x, const double *p) const { ...@@ -230,8 +230,8 @@ double TheuerkaufFitter::Eval(const double *x, const double *p) const {
// Evaluate internal background // Evaluate internal background
sum += std::accumulate( sum += std::accumulate(
std::reverse_iterator<const double *>(p + fNumParams - 1), std::reverse_iterator<const double *>(p + fNumParams),
std::reverse_iterator<const double *>(p + fNumParams - fIntBgDeg), std::reverse_iterator<const double *>(p + fNumParams - fIntBgDeg - 1),
0.0, [x = *x](double bg, double param) { return bg * x + param; }); 0.0, [x = *x](double bg, double param) { return bg * x + param; });
// Evaluate peaks // Evaluate peaks
...@@ -249,8 +249,8 @@ double TheuerkaufFitter::EvalBg(const double *x, const double *p) const { ...@@ -249,8 +249,8 @@ double TheuerkaufFitter::EvalBg(const double *x, const double *p) const {
// Evaluate internal background // Evaluate internal background
sum += std::accumulate( sum += std::accumulate(
std::reverse_iterator<const double *>(p + fNumParams - 1), std::reverse_iterator<const double *>(p + fNumParams),
std::reverse_iterator<const double *>(p + fNumParams - fIntBgDeg), std::reverse_iterator<const double *>(p + fNumParams - fIntBgDeg - 1),
0.0, [x = *x](double bg, double param) { return bg * x + param; }); 0.0, [x = *x](double bg, double param) { return bg * x + param; });
// Evaluate steps in peaks // Evaluate steps in peaks
...@@ -296,9 +296,8 @@ TF1 *TheuerkaufFitter::GetBgFunc() { ...@@ -296,9 +296,8 @@ TF1 *TheuerkaufFitter::GetBgFunc() {
return fBgFunc.get(); return fBgFunc.get();
} }
//! Do the fit, using the given background function
void TheuerkaufFitter::Fit(TH1 &hist, const Background &bg) { void TheuerkaufFitter::Fit(TH1 &hist, const Background &bg) {
//! Do the fit, using the given background function
// Refuse to fit twice // Refuse to fit twice
if (IsFinal()) { if (IsFinal()) {
return; return;
...@@ -309,10 +308,9 @@ void TheuerkaufFitter::Fit(TH1 &hist, const Background &bg) { ...@@ -309,10 +308,9 @@ void TheuerkaufFitter::Fit(TH1 &hist, const Background &bg) {
_Fit(hist); _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) { 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 // Refuse to fit twice
if (IsFinal()) { if (IsFinal()) {
return; return;
...@@ -323,9 +321,8 @@ void TheuerkaufFitter::Fit(TH1 &hist, int intBgDeg) { ...@@ -323,9 +321,8 @@ void TheuerkaufFitter::Fit(TH1 &hist, int intBgDeg) {
_Fit(hist); _Fit(hist);
} }
//! Private: worker function to actually do the fit
void TheuerkaufFitter::_Fit(TH1 &hist) { void TheuerkaufFitter::_Fit(TH1 &hist) {
//! Private: worker function to actually do the fit
// Allocate additional parameters for internal polynomial background // Allocate additional parameters for internal polynomial background
// Note that a polynomial of degree n has n+1 parameters! // Note that a polynomial of degree n has n+1 parameters!
if (fIntBgDeg >= 0) { if (fIntBgDeg >= 0) {
...@@ -526,7 +523,7 @@ void TheuerkaufFitter::_Fit(TH1 &hist) { ...@@ -526,7 +523,7 @@ void TheuerkaufFitter::_Fit(TH1 &hist) {
double sumFreeVol = sumVol; double sumFreeVol = sumVol;
auto ampIter = amps.begin(); auto ampIter = amps.begin();
for (auto &peak : fPeaks) { for (auto &peak : fPeaks) {
if (peak.fVol.IsFree()) { if (!peak.fVol.IsFree()) {
sumFreeAmp -= *(ampIter++); sumFreeAmp -= *(ampIter++);
sumFreeVol -= peak.fVol._Value(); sumFreeVol -= peak.fVol._Value();
} }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment