FitFunctions.h

Go to the documentation of this file.

 
class Fitf {
public:
    double operator() (double* x, double* p) {
        double sum = 0;
        int kmax = static_cast<int>(ceil(p[1])) * 10;
 
        double mu = p[1];
        double lambda = p[2];
 
        double sigma0 = p[3];
        double sigma1 = p[4];
 
        double G = p[5];
        double B = p[6];
 
        for (int kint = 0; kint <= kmax; kint++) {
            double k = static_cast<double>(kint);
            double sigmaK = sqrt(sigma0 * sigma0 + k * sigma1 * sigma1);
            //generalized poisson envelope
            double gp = mu * TMath::Power((mu + k * lambda), k - 1) * TMath::Exp(-(mu + k * lambda)) / TMath::Factorial(kint);
 
            sum += p[0] * gp * (1. / sqrt(2. * TMath::Pi()) / sigmaK) * TMath::Exp(-TMath::Power(((x[0] - (k * G + B)) / sqrt(2) / sigmaK), 2.));
        }
        return sum;
    };
};
 
class Fitf_full {
public:
    double operator() (double* x, double* p) {
        double sum = 0;
        int kmax = static_cast<int>(ceil(p[1])) * 10;
 
        double mu = p[1];
        double lambda = p[2];
 
        double sigma0 = p[3];
        double sigma1 = p[4];
 
        double G = p[5];
        double B = p[6];
 
        double alpha = p[7];
        double beta = p[8];
 
        for (int kint = 0; kint <= kmax; kint++) {
            //numbe of fired cells
            double k = static_cast<double>(kint);
            //pulse width
            double sigmaK = sqrt(sigma0 * sigma0 + k * sigma1 * sigma1);
            double gausnormsigmak = 1. / (sqrt(2. * TMath::Pi()) * sigmaK);
            //gauss peak
            double gauss = gausnormsigmak * TMath::Exp(-1. * TMath::Power(x[0] - (k * G + B), 2.) / (sigmaK * sigmaK * 2.));
 
            //generalized poisson envelope
            double gp = p[0] * mu * TMath::Power((mu + k * lambda), k - 1) * TMath::Exp(-(mu + k * lambda)) / TMath::Factorial(kint);
 
            if (kint == 0) {
                sum += (gp * gauss);
            }
            else {
                //after-pulse + delayed cross talk
                double bk0 = TMath::Power(1. - alpha, k);
                double bk1 = TMath::Factorial(kint) / TMath::Factorial(kint - 1) * alpha * TMath::Power(1. - alpha, k - 1);
                double pk1 = TMath::Exp(-1. * (x[0] - (k * G + B)) / beta) * gausnormsigmak / beta * sigmaK * sqrt(TMath::Pi() / 2) * (TMath::Erf((x[0] - (k * G + B)) / (sqrt(2) * sigmaK)) + 1.);
 
 
                if (kint == 1) sum += gp * (bk0 * gauss + bk1 * pk1);
                else {
                    double api2k = 0.;
                    for (int ii = 2; ii < kint; ii++) {
                        double iid = static_cast<double>(ii);
                        double bkialpha = TMath::Factorial(kint) / (TMath::Factorial(ii) * TMath::Factorial(kint - ii)) * TMath::Power(alpha, iid) * TMath::Power((1. - alpha), (k - iid));
                        double dpkidph = 0;
                        if (x[0] > (k * G + B)) dpkidph = TMath::Power((x[0] - (k * G + B)), (iid - 1.)) / (TMath::Factorial(ii - 1) * TMath::Power(beta, iid)) * TMath::Exp(-1. * (x[0] - (k * G + B)) / beta);
                        api2k += bkialpha * dpkidph;
                    }
                    sum += gp * (bk0 * gauss + bk1 * pk1 + api2k);
                }
            }
        }
        return sum;
    };
};
 
class Fitf_biased {
public:
    double operator() (double* x, double* p) {
        double sum = 0;
        int kmax = static_cast<int>(ceil(p[1])) * 10;
 
        double mu = p[1];
        double lambda = p[2];
 
        double sigma0 = p[3];
        double sigma1 = p[4];
 
        double G = p[5];
        double B = p[6];
 
        double a_ped = p[7];
        double x_ped = p[8];
 
        for (int kint = 0; kint <= kmax; kint++) {
            double k = static_cast<double>(kint);
            double sigmaK = sqrt(sigma0 * sigma0 + k * sigma1 * sigma1);
            //generalized poisson envelope
            double gp = mu * TMath::Power((mu + k * lambda), k - 1) * TMath::Exp(-(mu + k * lambda)) / TMath::Factorial(kint);
 
            if (kint == 0) {
                sum += a_ped * p[0] * gp * (1. / sqrt(2. * TMath::Pi()) / sigmaK) * TMath::Exp(-TMath::Power(((x[0] - x_ped) / sqrt(2) / sigmaK), 2));
            }
            else {
                sum += p[0] * gp * (1. / sqrt(2. * TMath::Pi()) / sigmaK) * TMath::Exp(-TMath::Power(((x[0] - (k * G + B)) / sqrt(2) / sigmaK), 2));
            }
        }
        return sum;
    };
};
 
class Fitf_PMT {
public:
    double operator() (double* x, double* p) {
        double pmt_charge_spectrum = 0.;
        int kmax = static_cast<int>(ceil(p[5])) * 10;
 
        for (int kint = 0; kint <= kmax; kint++) {
            double k = static_cast<double>(kint);
 
            double poiss = TMath::Power(p[5], k) * TMath::Exp(-p[5]) / TMath::Factorial(kint);
 
            double normgn_term = 1.;
            if (kint != 0) normgn_term /= (p[6] * TMath::Sqrt(2. * TMath::Pi() * k));
            double gn_term = (1. - p[1]) * normgn_term;
            if (kint != 0) gn_term *= TMath::Exp(-1. * TMath::Power(x[0] - k * p[7] - p[4], 2.) / (2. * k * p[6] * p[6]));
            else if (x[0] != p[4]) gn_term *= 0.; // delta function
 
            double Qn = p[4] + k * p[7];
            double sigman = TMath::Sqrt(p[3] * p[3] + k * p[6] * p[6]);
            double ignxe_term_exp = p[2] / 2. * TMath::Exp(-1. * p[2] * (x[0] - Qn - p[2] * sigman * sigman));
 
            double sgn_arg = x[0] - Qn - sigman * sigman * p[2];
            double arg_sgn = 1.;
            if (sgn_arg < 0.) arg_sgn = -1.;
            double ignxe_term_erf = TMath::Erf(fabs(p[4] - Qn - sigman * sigman * p[2]) / (sigman * TMath::Sqrt(2.))) + arg_sgn * TMath::Erf(fabs(sgn_arg) / (sigman * TMath::Sqrt(2.)));
 
            pmt_charge_spectrum += p[0] * poiss * (gn_term + p[1] * ignxe_term_exp * ignxe_term_erf);
        }
        if (pmt_charge_spectrum < 0.) pmt_charge_spectrum = 0.;
        return pmt_charge_spectrum;
    };
};
 
class Fitf_PMT_pedestal {
public:
    double operator() (double* x, double* p) {
        double pmt_charge_spectrum = 0.;
        int kmax = static_cast<int>(ceil(p[5])) * 10;
 
        for (int kint = 0; kint <= kmax; kint++) {
            double k = static_cast<double>(kint);
 
            double poiss = TMath::Power(p[5], k) * TMath::Exp(-p[5]) / TMath::Factorial(kint);
 
            double normgn_term = 1.;
            if (kint != 0) normgn_term /= (p[6] * TMath::Sqrt(2. * TMath::Pi() * k));
            else normgn_term *= p[8] / (p[3] * TMath::Sqrt(2. * TMath::Pi()));
            double gn_term = (1. - p[1]) * normgn_term;
            if (kint != 0) gn_term *= TMath::Exp(-1. * TMath::Power(x[0] - k * p[7] - p[4], 2.) / (2. * k * p[6] * p[6]));
            else  gn_term *= TMath::Exp(-1. * TMath::Power(x[0] - p[4], 2.) / (2. * p[3] * p[3]));
 
            double Qn = p[4] + k * p[7];
            double sigman = TMath::Sqrt(p[3] * p[3] + k * p[6] * p[6]);
            double ignxe_term_exp = p[2] / 2. * TMath::Exp(-1. * p[2] * (x[0] - Qn - p[2] * sigman * sigman));
 
            double sgn_arg = x[0] - Qn - sigman * sigman * p[2];
            double arg_sgn = 1.;
            if (sgn_arg < 0.) arg_sgn = -1.;
            double ignxe_term_erf = TMath::Erf(fabs(p[4] - Qn - sigman * sigman * p[2]) / (sigman * TMath::Sqrt(2.))) + arg_sgn * TMath::Erf(fabs(sgn_arg) / (sigman * TMath::Sqrt(2.)));
 
            pmt_charge_spectrum += p[0] * poiss * (gn_term + p[1] * ignxe_term_exp * ignxe_term_erf);
        }
        if (pmt_charge_spectrum < 0.) pmt_charge_spectrum = 0.;
        return pmt_charge_spectrum;
    };
};
 
class Fitf_PMT_ideal {
public:
    double operator() (double* x, double* p) {
        double pmt_charge_spectrum = 0.;
 
        int kmax = static_cast<int>(ceil(p[1])) * 10;
 
        for (int kint = 1; kint <= kmax; kint++) {
            double k = static_cast<double>(kint);
 
            double poiss = TMath::Power(p[1], k) * TMath::Exp(-p[1]) / TMath::Factorial(kint);
 
            double norm = 1. / (p[2] * TMath::Sqrt(2. * TMath::Pi() * k));
 
            double gauss = TMath::Exp(-1. * TMath::Power(x[0] - k * p[3], 2) / (2. * k * p[2] * p[2]));
 
            pmt_charge_spectrum += p[0] * poiss * norm * gauss;
        }
        return pmt_charge_spectrum;
    };
};
 
class Fitf_langaus {
public:
    double operator() (double* x, double* par) {
        // Numeric constants
        double invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
        double mpshift = -0.22278298;       // Landau maximum location
 
        // Control constants
        double np = 100.0;      // number of convolution steps
        double sc = 5.0;      // convolution extends to +-sc Gaussian sigmas
 
        // Variables
        double xx;
        double mpc;
        double fland;
        double sum = 0.0;
        double xlow, xupp;
        double step;
        double i;
 
 
        // MP shift correction
        mpc = par[1] - mpshift * par[0];
 
        // Range of convolution integral
        xlow = x[0] - sc * par[3];
        xupp = x[0] + sc * par[3];
 
        step = (xupp - xlow) / np;
 
        // Convolution integral of Landau and Gaussian by sum
        for (i = 1.0; i <= np / 2; i++) {
            xx = xlow + (i - .5) * step;
            fland = TMath::Landau(xx, mpc, par[0]) / par[0];
            sum += fland * TMath::Gaus(x[0], xx, par[3]);
 
            xx = xupp - (i - .5) * step;
            fland = TMath::Landau(xx, mpc, par[0]) / par[0];
            sum += fland * TMath::Gaus(x[0], xx, par[3]);
        }
 
        return (par[2] * step * sum * invsq2pi / par[3]);
    };
};
 
class Fitf_plus_DC {
public:
    double operator() (double* x, double* p) {
        double sum = 0;
        int kmax = static_cast<int>(ceil(p[1])) * 10;
 
        double mu = p[1];
        double lambda = p[2];
 
        double sigma0 = p[3];
        double sigma1 = p[4];
 
        double G = p[5];
        double B = p[6];
 
        double mu_dc = p[7];
        double N0_dc = p[8];
 
        for (int kint = 0; kint <= kmax; kint++) {
            double k = static_cast<double>(kint);
            double sigmaK = sqrt(sigma0 * sigma0 + k * sigma1 * sigma1);
            //generalized poisson envelope
            double gp = ((1. - N0_dc) * (mu * TMath::Power((mu + k * lambda), k - 1) * TMath::Exp(-(mu + k * lambda))) + N0_dc * (mu_dc * TMath::Power((mu_dc + k * lambda), k - 1) * TMath::Exp(-(mu_dc + k * lambda)))) / TMath::Factorial(kint);
 
            sum += p[0] * gp * (1. / sqrt(2. * TMath::Pi()) / sigmaK) * TMath::Exp(-TMath::Power(((x[0] - (k * G + B)) / sqrt(2) / sigmaK), 2));
        }
        return sum;
    };
};
 
class Fitf_periodic_gauss {
public:
    double operator() (double* x, double* p) {
        double k = 0;
        int k_pm = 2; // how many Gaussians to consider -> k_pm=2 means 2+1+2=5 Gaussians
        double A = p[0];
        double mu = p[1];
        double sigma = p[2];
        double B = p[3];
 
        double gauss = 0.;
 
        for (int kint = -k_pm; kint <= k_pm; kint++) {
            k = static_cast<double>(kint);
            // periodic gauss
            gauss += A * TMath::Exp(-TMath::Power(x[0] - mu + k * 360., 2.) / (2. * sigma * sigma));
        }
        return gauss + B;
    };
};
 
class Fitf_langaus_int {
public:
    double operator() (double* x, double* par) {
        // Control constants
        Double_t sc = 5.;        // integral extends to MPV-sc*(Gaussian sigma + Landau width)
 
        // Variables
        Double_t start = par[1] - sc * (par[0] + par[3]);
 
        Fitf_langaus langaus;
        TF1* langaus_fun = new TF1("fun", langaus, start - 1, x[0] + 1, 4);
        langaus_fun->SetParameter(0, par[0]);
        langaus_fun->SetParameter(1, par[1]);
        langaus_fun->SetParameter(2, par[2]);
        langaus_fun->SetParameter(3, par[3]);
 
        return (par[4] * (1. - langaus_fun->Integral(start, x[0])));
    };
};
 
class Fitf_langaus_poisson {
public:
 
    double operator() (double* x, double* par) {
        // Range of the sum (should be significantly larger than the mean number of MIPs)
        int kmax = static_cast<int>(ceil(par[5])) * 10;
 
        Fitf_langaus langaus;
        TF1* langaus_fun = new TF1("fun", langaus, 0, 1, 4);
        langaus_fun->SetParameter(0, par[0]);
        langaus_fun->SetParameter(2, par[2]);
 
        double sum = 0.;
        for (int kint = 1; kint <= kmax; kint++) {
            double k = static_cast<double>(kint);
 
            double peak_shift = (k - 1) * par[6];
            langaus_fun->SetParameter(1, par[1] + peak_shift);
 
            double sigmaK = sqrt(par[3] * par[3] + (k - 1) * par[4] * par[4]);
            langaus_fun->SetParameter(3, sigmaK);
 
            //poisson envelope
            double poiss = TMath::Power(par[5], k) / (TMath::Factorial(kint) * (TMath::Exp(par[5]) - 1.));
 
            sum += langaus_fun->Eval(x[0]) * poiss;
        }
        return sum;
    };
};
 
class Fitf_exp_gauss {
public:
    double operator() (double* x, double* par) {
        double tau_eff = par[0];
        double sigma_gauss = par[1];
        double t_0 = par[2];
        double norm = par[3];
 
        double e_g = norm / (2 * TMath::Abs(tau_eff)) *
            TMath::Exp((sigma_gauss * sigma_gauss + 2 * t_0 * tau_eff - 2 * tau_eff * x[0]) /
                (2 * tau_eff * tau_eff)) *
            TMath::Erfc((sigma_gauss * sigma_gauss + tau_eff * (t_0 - x[0])) /
                (TMath::Sqrt2() * TMath::Abs(tau_eff) * sigma_gauss));
 
        return e_g;
    };
};