resolution.h File Reference

#include "Core/CommonTools/interface/binnings.h"
#include "Core/CommonTools/interface/toolbox.h"
#include <TH1.h>
#include <iostream>
#include <cmath>

Include dependency graph for resolution.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
double	safelog (double x)
double	id (double x)
TH1 *	Derivative (TH1 *h, double(*f)(double)=id)
unique_ptr< TH1 >	DerivativeFivePointStencil (const unique_ptr< TH1 > &h, double(*f)(double)=id)
TH1 *	SecondDerivative (TH1 *h, double(*f)(double)=id)
double	mNSC (double *x, double *p)
double	GausExp (double *x, double *p)
double	Novosibirsk (double *x, double *p)
double	CrystalBall (double *x, double *p)
double	ModifiedGaussianCore (double *x, double *p)

Variables
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

Function Documentation

CrystalBall()

double CrystalBall	(	double *	x,
		double *	p
	)

Resolution with modified Gaussian with simple power-law tail

Note: not very stable because of n**n behaviour -> the trick when fitting is to find the k by looking at the derivative of the log of the function, since the Gaussian should be a straight line...

NB: code should be carefully checked ...

https://en.wikipedia.org/wiki/Crystal_Ball_function

NOTE: prefer the more general implementation

{
    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);
    }
}

Derivative()

TH1* Derivative	(	TH1 *	h,
		double(*)(double)	f = id
	)

Calculate numerical derivative of a histogram

https://en.wikipedia.org/wiki/Numerical_differentiation

{
    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;
}

DerivativeFivePointStencil()

unique_ptr<TH1> DerivativeFivePointStencil	(	const unique_ptr< TH1 > &	h,
		double(*)(double)	f = id
	)

Calculate numerical derivative of a histogram

https://en.wikipedia.org/wiki/Numerical_differentiation

{
    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;
}

GausExp()

double GausExp	(	double *	x,
		double *	p
	)

Resolution with modified Gaussian with simple exponentional tail

Note: more stable than Crystal-Ball, but only one parameter

https://home.fnal.gov/~souvik/GaussExp/GaussExp.pdf

{
    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));
    }
}

id()

double id ( double  x )

inline

{ return x; };

mNSC()

double mNSC	(	double *	x,
		double *	p
	)

Resolution width as a function of gen pt

• p[0] = Noise
• p[1] = Stochastic
• p[2] = Constant
• p[3] = modification parameter

Note that the freedom for the sign of the first term is also a modification of the original NSC function.

{
    return sqrt( p[0]*abs(p[0])/(x[0]*x[0]) + p[1]*p[1]/pow(x[0],p[3]) + p[2]*p[2]);
}

ModifiedGaussianCore()

double ModifiedGaussianCore	(	double *	x,
		double *	p
	)

Advanced function to fit resolution including non-Gaussian deviations and pile-up jets.

{
    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);
}

Novosibirsk()

double Novosibirsk	(	double *	x,
		double *	p
	)

Asymmetric Gaussian, a.k.a. Novosibirsk

Todo:

find reference

{
    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);
}

safelog()

double safelog ( double  x )

inline

{ return x > 0 ? log(x) : 0; };

SecondDerivative()

TH1* SecondDerivative	(	TH1 *	h,
		double(*)(double)	f = id
	)

Calculate numerical second derivative of a histogram

Note: seems to work less good that just applying twice Derivative()...

https://en.wikipedia.org/wiki/Numerical_differentiation

{
    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;
}

Variable Documentation

nkL

const size_t nkL = 2

static

order of Chebyshev polynomial for global LHS tail transition of DCB

nkR

const size_t nkR = 2

static

order of Chebyshev polynomial for global RHS tail transition of DCB

nmu

const size_t nmu = 3

static

order of Chebyshev polynomial for global mean of DCB

nN

const size_t nN = 8

static

order of Chebyshev polynomial for global normalisation of DCB

nnL

const size_t nnL = 5

static

order of Chebyshev polynomial for global LHS tail steepness of DCB

nnR

const size_t nnR = 5

static

order of Chebyshev polynomial for global RHS tail steepness of DCB