fit.h

Go to the documentation of this file.

#include <cmath>
 
#include <utility>
#include <optional>
#include <memory>
#include <iostream>
#include <iomanip>
 
#include <TF1.h>
#include <TFitResult.h>
#include <TH1.h>
#include <TString.h>
 
#include <darwin.h>
 
namespace DAS::JetEnergy {
 
struct AbstractFit {
 
    const std::unique_ptr<TH1>& h; 
    std::uint32_t status; 
    std::ostream& cout; 
 
    double * p, 
           * e; 
 
    std::unique_ptr<TF1> f; 
    std::optional<double> chi2ndf, chi2ndfErr; 
    std::pair<float, float> interval; 
 
    inline virtual bool good () const { return chi2ndf && status; }
 
protected:
 
    AbstractFit (const unique_ptr<TH1>& h, ostream& cout, int npars) :
        h(h), status(0), cout(cout),
        p(new double[npars]), e(new double[npars])
    { }
 
    virtual ~AbstractFit ()
    {
        for (int i = 0; i < f->GetNpar(); ++i) {
            if (isfinite(p[i])) continue;
            const char * what = Form("Fit parameter %d is not finite", i);
            BOOST_THROW_EXCEPTION( std::runtime_error(what) );
        }
        delete p;
        delete e;
    }
 
    void fit (std::pair<float,float>, const char *);
 
    inline virtual void Write (const char * name)
    {
        if (good()) // successful fit
            f->SetTitle(Form("%.1f", *chi2ndf));
        else {
            f->SetLineStyle(kDashed);
            f->SetTitle("fail");
        }
 
        f->SetParameters(p);
        f->SetParErrors (e);
        f->SetRange(interval.first, interval.second);
        f->Write(name);
    }
};
 
} // end of DAS::JetEnergy namespace
 
inline std::ostream& operator<< (std::ostream& s, const DAS::JetEnergy::AbstractFit& fit)
{
    using namespace std;
    s << fixed << setprecision(3);
    for (int i = 0; i < fit.f->GetNpar(); ++i) {
        double m, M;
        fit.f->GetParLimits(i, m, M);
        if (m < M) s << bold;
        s << setw(10) << fit.f->GetParameter(i);
        if (m < M) s << normal;
    }
    s << defaultfloat;
    return s;
}
 
void DAS::JetEnergy::AbstractFit::fit (std::pair<float,float> range,
                                       const char * fitName)
{
    using namespace std;
 
    // reinitialise the params to the best known fit
    f->SetParameters(p);
    f->SetParErrors (e);
    f->SetRange(range.first, range.second);
    for (int i = 0; i < f->GetNpar(); ++i) {
        double m, M;
        f->GetParLimits(i,m,M);
        if (p[i] >= m && p[i] <= M) continue;
        cerr << orange << setw(2) << i
                       << setw(10) << m
                       << setw(10) << M
                       << setw(10) << p[i]
                       << def << '\n';
    }
 
    // run the fit
    TFitResultPtr result = h->Fit(f.get(), "QS0NRE");
    // Q: quiet
    // S: return smart pointer to fit restult (automatic conversion to integer provides fit status)
    // 0: do not draw
    // N: do not store graphics
    // R: use fit range
 
    // get result
    double currentChi2ndf = f->GetChisquare()/f->GetNDF();
 
    bool success = result >= 0 && result % 10 == 0 && currentChi2ndf < 10000;
    //  status = migradStatus + 10*minosStatus + 100*hesseStatus + 1000*improveStatus.
    //  TMinuit returns 0 (for migrad, minos, hesse or improve) in case of success and 4 in case of error (see the documentation of TMinuit::mnexcm).
    //  For example, for an error only in Minos but not in Migrad a fitStatus of 40 will be returned.
    //  Minuit2 returns 0 in case of success and different values in migrad,minos or hesse depending on the error.
    //  See in this case the documentation of Minuit2Minimizer::Minimize for the migrad return status,
    //                                        Minuit2Minimizer::GetMinosError for the minos return status,
    //                                    and Minuit2Minimizer::Hesse for the hesse return status.
 
    bool improvement = !chi2ndf || abs(currentChi2ndf-1) <= abs(*chi2ndf-1);
 
    cout << h->GetTitle() << " with " << bold << fitName << normal << " in ["
         << fixed << setprecision(2) << range.first << "," << range.second << "]: " << defaultfloat
         << (success ? (improvement ? green : orange) : red) << *this
         << fixed << setprecision(3) << setw(15) << currentChi2ndf << defaultfloat
         << "\t(status = " << (int) result << ")" << def << endl;
 
    // skip if the new attempt isn't improving the performance
    if (!success || !improvement) return;
 
    // at this stage, we can consider that the fit has improved, 
    // therefore we save the current parameter estimates
    chi2ndf = currentChi2ndf;
    chi2ndfErr = sqrt(2./f->GetNDF());
    for (int i = 0; i < f->GetNpar(); ++i) {
        p[i] = f->GetParameter(i);
        e[i] = f->GetParError (i);
    }
    interval = range;
}