resolution.h

Go to the documentation of this file.

#include "Core/CommonTools/interface/binnings.h"
#include "Core/CommonTools/interface/toolbox.h"
 
#include <TH1.h>
 
#include <iostream>
#include <cmath>
 
using namespace DAS;
using namespace std;
 
static const size_t nN  = 8; 
static const size_t nmu = 3; 
static const size_t nkL = 2; 
static const size_t nnL = 5; 
static const size_t nkR = 2; 
static const size_t nnR = 5; 
 
inline double safelog (double x) { return x > 0 ? log(x) : 0; };
 
inline double id (double x) { return x; };
 
TH1 * Derivative (TH1 * h, double (*f)(double) = id)
{
    const char * name = Form("derivative_%s", h->GetName());
    TH1 * d = dynamic_cast<TH1*>(h->Clone(name));
    d->Reset();
    for (int i = 2; i < d->GetNbinsX(); ++i) {
        auto x1 = h->GetBinCenter(i-1),
             x2 = h->GetBinCenter(i+1),
             f1 = f(h->GetBinContent(i-1)),
             f2 = f(h->GetBinContent(i+1));
        auto dx = x2-x1;
        if (abs(dx) < 1e-20) continue;
        auto df = f2-f1;
        d->SetBinContent(i, df/dx);
    }
    return d;
}
 
unique_ptr<TH1> DerivativeFivePointStencil (const unique_ptr<TH1>& h, double (*f)(double) = id)
{
    const char * name = Form("derivative_%s", h->GetName());
    auto d = unique_ptr<TH1>( dynamic_cast<TH1*>(h->Clone(name)) );
    d->Reset();
    for (int i = 3; i < d->GetNbinsX()-1; ++i) {
        auto w = d->GetBinWidth(i);
        //assert(abs(w - d->GetBinWidth(i)) < 1e-4); /// \todo not necessary met in the tails, where bins are merged....
        auto f1 = f(h->GetBinContent(i-2)),
             f2 = f(h->GetBinContent(i-1)),
             f3 = f(h->GetBinContent(i+1)),
             f4 = f(h->GetBinContent(i+2));
        auto dx = 12*w;
        if (abs(dx) < 1e-20) continue;
        auto df = -f4 + 8*f3 - 8*f2 + f1;
        d->SetBinContent(i, df/dx);
    }
    return d;
}
 
TH1 * SecondDerivative (TH1 * h, double (*f)(double) = id)
{
    const char * name = Form("secondderivative_%s", h->GetName());
    TH1 * d2 = dynamic_cast<TH1*>(h->Clone(name));
    d2->Reset();
    for (int i = 2; i < d2->GetNbinsX(); ++i) {
        auto x0 = h->GetBinCenter(i-1),
             x1 = h->GetBinCenter(i  ),
             x2 = h->GetBinCenter(i+1),
             f0 = f(h->GetBinContent(i-1)),
             f1 = f(h->GetBinContent(i  )),
             f2 = f(h->GetBinContent(i+1));
        auto dx2 = (x2-x1)*(x1-x0);
        if (abs(dx2) < 1e-20) continue;
        auto d2f = f2-2*f1+f0;
        d2->SetBinContent(i, d2f/dx2);
    }
    return d2;
}
 
double mNSC (double * x, double * p)
{
    return sqrt( p[0]*abs(p[0])/(x[0]*x[0]) + p[1]*p[1]/pow(x[0],p[3]) + p[2]*p[2]);
}
 
double GausExp (double * x, double * p)
{
    double N = p[0],
           mu = p[1],
           sigma = p[2],
           kR = p[3];
 
    auto z = [&](double x) { return (x - mu)/sigma; };
 
    double zx = z(*x);
 
    if (*x > kR) {
        double zk = z(kR);
        return N*exp(0.5*pow(zk,2)-zk*zx);
    }
    else {
        return N*exp(-0.5*pow(zx,2));
    }
}
 
double Novosibirsk (double * x, double * p)
{
    double N = p[0],
           mu = p[1],
           sigma = p[2],
           tau = p[3];
 
    if (sigma <= 0) return 0;
 
    double y = (*x-mu)/sigma;
 
    static const double C = sqrt(log(4.));
    double coeff = sinh(C*tau) / C;
 
    double z2 = pow( log(1.+coeff*y)/tau, 2) + pow(tau,2);
    return N*exp(-0.5*z2);
}
 
[[ deprecated ]]
double CrystalBall (double * x, double * p)
{
    double N = p[0],
           mu = p[1],
           sigma = p[2],
           k = p[3],
           n = p[4];
 
    double z =    (*x - mu)/sigma,
           a = abs( k - mu)/sigma;
 
    double expa2 = exp(-pow(a,2)/2.);
 
    double C = n * expa2 / (a*(n-1)),
           D = sqrt(M_PI/2.)*(1 + erf(a/sqrt(2)));
 
    N /= sigma * (C + D);
 
    if (*x > k) {
        return N*exp(-pow(z,2)/2.);
    }
    else {
        double A = pow(n/a, n) * expa2,
               B = n/a - a;
        return N*A*pow(B-z,-n);
    }
}
 
double ModifiedGaussianCore (double * x, double * p)
{
    double N = p[0], // Gaussian core
           mu = p[1],
           sigma = max(1e-8,abs(p[2])),
           B = p[3],
           C = p[4];
 
    return N*exp(-  pow((*x-mu)/sigma,2)/2
                 -B*pow((*x-mu)      ,3)/3
                 -C*pow((*x-mu)      ,4)/4);
}