Greta.h

Go to the documentation of this file.

#ifndef __GRETA__
#define __GRETA__
#include <cstdlib>
#include <vector>
#include <iostream>
#include <iomanip>
#include <cassert>
#include <fstream>
 
#include <TF1.h>
#include <TH1.h>
#include <TString.h>
#include <TFitResult.h>
 
using namespace std;
 
static const float feps = numeric_limits<float>::epsilon();
static const double deps = numeric_limits<double>::epsilon();
static const double dmax = numeric_limits<double>::max();
 
namespace DAS {
 
inline double identity (double x) { return x; };
 
namespace Chebyshev {
 
    template<const size_t n> struct Term {
        static double Eval (const double x)
        {
            const size_t m = n-1, k = n-2;
            return 2*x*Term<m>::Eval(x)-Term<k>::Eval(x);
        }
    };
    template<> struct Term<1> {
        static double Eval (double x) { return x; }
    };
    template<> struct Term<0> {
        static double Eval (double x) { return 1.; }
    };
 
    template<const size_t n> struct Polynom {
        static double Eval (double x, const double p[])
        {
            const size_t k = n-1;
            return p[n]*Term<n>::Eval(x) + Polynom<k>::Eval(x,p);
        }
    };
    template<> struct Polynom<0> {
        static double Eval (double x, const double p[]) { return p[0]; }
    };
 
    template<const size_t maxnPt, const size_t nPt, const size_t nY> struct Polynom2D {
        static double Eval (double xPt, double xY, const double p[])
        {
            const size_t kY = nY-1;
            const double * pPt = p + nY * maxnPt * sizeof(double);
            const double q_nY = Polynom<nPt>::Eval(xPt, pPt);
            return q_nY * Term<nY>::Eval(xY) 
                    + Polynom2D<maxnPt, nPt, kY>::Eval(xPt, xY, pPt);
        }
    };
    template<const size_t maxnPt, const size_t nPt> struct Polynom2D<maxnPt,nPt,0> {
        static double Eval (double xPt, double xY, const double p[])
        {
            return Chebyshev::Polynom<nPt>::Eval(xPt, p);
        }
    };
}
 
template<const size_t d, double (*Log)(double) = log, double (*Exp)(double) = exp>
struct Thunberg {
 
    const size_t maxdegree, npars;
    const double m, M;
 
    Thunberg (double mi, double Ma) :
        maxdegree(d), npars(d+1), m(Log(mi)), M(Log(Ma))
    {
        static_assert(d > 1, "degree must be greater than 1");
    }
 
    double operator() (const double * x, const double * p)
    {
        double nx = -1 + 2* (Log(x[0]) - m) / (M - m);
        double chebyshev = Chebyshev::Polynom<d>::Eval(nx,p);
        double result = Exp(chebyshev);
        return result;
    }
};
 
 
template<double (*Log)(double) = log, double (*Exp)(double) = exp>
struct Greta {
 
    const int degree, npars;
    const double m, M;
 
    [[ deprecated ]]
    Greta (int d, double mi, double Ma) :
        degree(d), npars(d+1), m(mi), M(Ma)
    { }
 
    static double T (double x, int i)
    {
        if (i == 0) return 1.;
        if (i == 1) return x;
        return 2*x*T(x, i-1)-T(x, i-2);
    }
 
    static const char * T (const char * x, int i)
    {
        if (i == 0) return Form("1.");
        if (i == 1) return x;
        const char * t1 = T(x, i-1), * t2 = T(x, i-2);
        return Form("2*%s*%s-%s", x, t1, t2);
    }
 
    inline double normalise (const double *x) const
    {
        return -1 + 2* (Log(*x) - Log(m)) / (Log(M) - Log(m));
    }
 
    double operator() (const double *x, const double *p) const
    {
        double nx = normalise(x);
        double result = 0;
        for (int i = 0; i <= degree; ++i)
            result += p[i]*T(nx,i);
            //result += p[i]*pow(nx,i);
        return Exp(result);
    }
 
    static const char * formula (double mi, double Ma, int deg)
    {
        TString result = "[2]"; 
        TString nx = "(-1+2*(log(x)-log([0]))/(log([1])-log([0])))";
        for (int i = 1; i <= deg; ++i) 
            result += Form("+[%d]*%s", i+2, T(nx, i));
        result = "exp(" + result + ")";
        return result.Data();
    }
 
    const char * formula ()
    {
        return formula(m, M, degree);
    }
 
};
 
template<double (*Log)(double) = log, double (*Exp)(double) = exp>
struct GretaPolynomial {
    double * const p;
    const Greta<Log,Exp> g;
    [[ deprecated ]]
    GretaPolynomial (double * const pars, int Npars, double m, double M) :
        p(pars), g(Npars-1, m, M) { }
    double operator() (double x) const {
        x = max(g.m, x);
        x = min(g.M, x);
        return g(&x,p);
    }
};
 
template<double (*Log)(double) = log, double (*Exp)(double) = exp>
class GretaPolynomialReader {
 
    vector<double> ymins;
    vector<GretaPolynomial<Log,Exp>> pols; 
 
public:
    void read (TString fname)
    {
        cout << "Reading " << fname << endl;
        pols.clear();
        ymins.clear();
        ifstream infile;
        infile.open(fname.Data());
        if (infile.fail()) {
            cerr << "Failure at opening the file " << fname << '\n';
            exit(EXIT_FAILURE);
        }
 
        while (infile.good()) { 
 
            // getting values
            double ymin, ptmin, ptmax;
            int Npars;
            infile >> ymin >> ptmin >> ptmax >> Npars;
            double * pars = new double[Npars];
            for (int i = 0; i < Npars; ++i) infile >> pars[i];
            
            if (infile.eof()) break;
 
            // instantiating
            if (ymins.size() == 0) { assert(abs(ymin) < deps);    }
            else                   { assert(ymin > ymins.back()); }
            ymins.push_back(ymin);
            GretaPolynomial<Log,Exp> gp(pars, Npars, ptmin, ptmax);
            pols.push_back(gp);
            //pols.insert( { ymin, gp } );
 
            // print-out
            cout << std::setw(3) << ymin
                 << std::setw(8) << (int) ptmin
                 << std::setw(8) << (int) ptmax
                 << std::setw(3) << Npars;
            for (int j = 0; j < Npars; ++j) std::cout << std::setw(15) << pars[j];
            cout << '\n';
        } 
        cout << flush;
 
        assert(ymins.size() > 0);
        assert(ymins.size() == pols.size());
 
        infile.close(); 
    }
 
    double operator() (double pt, double y)
    {
        y = abs(y);
        size_t iy = ymins.size();
        do --iy;
        while (ymins.at(iy) > y);
        auto& g = pols.at(iy);
        return g(pt);
    }
 
};
 
 
template<const size_t d>
[[ deprecated ]]
TF1 * GetSmoothFit (TH1 * h, 
        int ifm, 
        int ifM, 
        int nSigmaStop = 0, 
        const char * options = "Q0SNRE", 
        bool autoStop = false)
{
    assert(ifm > 0);
    assert(ifM <= h->GetNbinsX());
    assert(h->GetDimension() == 1);
 
    double contentfm = h->GetBinContent(ifm),
           contentfM = h->GetBinContent(ifM);
 
    double m = h->GetBinLowEdge(ifm), M = h->GetBinLowEdge(ifM+1);
 
    cout << "From bin " << ifm << " (m = " << m << ", content = " << contentfm << ")"
         << " to bin "  << ifM << " (M = " << M << ", content = " << contentfM << ")" << endl;
 
    if (abs(contentfm) <= deps || abs(contentfM) <= deps) {
        cerr << "Fit cannot start with empty bins\n";
        //h->Print("all");
        return nullptr;
    }
 
    double errorfm = h->GetBinError(ifm),
           errorfM = h->GetBinError(ifM);
 
    if (abs(errorfm) <= deps || abs(errorfM) <= deps) {
        cerr << "There must be uncertainties for the fit\n";
        //h->Print("all");
        return nullptr;
    }
 
    double fm = log(contentfm),
           fM = log(contentfM);
 
    double x2n = dmax, error = dmax;
 
    static int iname = 0; ++iname;
    //Greta greta(maxdegree, m, M);
    Thunberg<d, log, exp> greta(m, M);
    const size_t maxdegree = greta.maxdegree;
    TF1 * f = new TF1(Form("greta%d", iname), greta, m, M, greta.npars);
 
    // setting initial parameters
    f->SetParameter(0,(fM+fm)/2);
    f->SetParameter(1,(fM-fm)/2);
    for (size_t degree = 2; degree <= maxdegree; ++degree)
        f->FixParameter(degree,0);
 
    auto dump = [f](int degree) {
        cout << "parameters:\t";
        for (int i = 0; i <= degree; ++i) cout << f->GetParameter(i) << ' ';
        cout << endl;
    };
 
    //running
    for (size_t degree = 2; degree <= maxdegree; ++degree) {
 
        dump(degree-1);
        f->ReleaseParameter(degree);
 
        TFitResultPtr r = h->Fit(f, options);
 
        double chi2 = r->Chi2();
        if (!isnormal(chi2)) {
            h->Print("all");
            cerr << "Fit aborted (chi2 = " << chi2 << ")\n";
            continue;
        }
        double ndf = r->Ndf();
 
        int status = r->Status();
        if (status != 0) {
            h->Print("all");
            cerr << "Fit aborted (status = " << status << ")\n";
            continue;
        }
 
        auto lastx2n = x2n;
        x2n = chi2/ndf;
        error = sqrt(2./ndf);
        cout << "\e[1m" << degree << ' ' << status << '\t' << chi2 << ' ' << ndf << '\t' << x2n << "\u00B1" << error << "\e[0m" << endl;
 
        if (ndf == 2) break;
        if (autoStop && x2n > lastx2n) {
            cerr << "Chi2/ndf has increased -- stepping back and stopping!\n";
            f->FixParameter(degree,0);
            --degree;
            h->Fit(f, options);
            break;
        }
        if (abs(x2n-1) <= nSigmaStop*error) break;
    }
    dump(maxdegree);
 
    return f;
}
 
template<const size_t d>
[[ deprecated ]]
TF1 * GetSmoothFit (TH1 * h, 
        double m, 
        double M, 
        int nSigmaStop = 0, 
        const char * options = "Q0SNRE", 
        bool autoStop = false)
{
    assert(m < M);
 
    int ifm = h->FindBin(m+feps),
        ifM = h->FindBin(M-feps);
    cout << "Axis interval (" << m << ',' << M << ") to bin interval (" << ifm << ',' << ifM << ')' << endl;
 
    assert(ifm < ifM);
 
    return GetSmoothFit<d>(h, ifm, ifM, nSigmaStop, options, autoStop);
}
 
} // end of DAS namespace
#endif