FitFunctions.cpp

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <map>

#include "FitFunctions.h"

using namespace std;


//ClassImp(FitFunctions);

double FitFunctions::fnc_dscb(double*xx,double*pp)
{
    double x   = xx[0];
    // gaussian core
    double N   = pp[0];//norm
    double mu  = pp[1];//mean
    double sig = pp[2];//variance
    // transition parameters
    double a1  = pp[3];
    double p1  = pp[4];
    double a2  = pp[5];
    double p2  = pp[6];

    double u   = (x-mu)/sig;
    double A1  = TMath::Power(p1/TMath::Abs(a1),p1)*TMath::Exp(-a1*a1/2);
    double A2  = TMath::Power(p2/TMath::Abs(a2),p2)*TMath::Exp(-a2*a2/2);
    double B1  = p1/TMath::Abs(a1) - TMath::Abs(a1);
    double B2  = p2/TMath::Abs(a2) - TMath::Abs(a2);

    double result(N);
    if      (u<-a1) result *= A1*TMath::Power(B1-u,-p1);
    else if (u<a2)  result *= TMath::Exp(-u*u/2);
    else            result *= A2*TMath::Power(B2+u,-p2);
    return result;
}


double FitFunctions::DinoFct(double*xx,double*pp)
{
    double x   = xx[0];

    int    NSubPeaks = (int)pp[0]; //Number of subpeaks in the peak range
    double BgFrom    = pp[1]; //First Channel for the Bg estimation
    double BgTo      = pp[2]; //Last Channel for the Bg estimation
    double BgdOff    = pp[3]; //Bg offset
    double BgdSlope  = pp[4]; //Bg slope

    double f_tot = 0.;

    if(x<BgFrom || x>BgTo) return 0.;
    else
    {
        double BGd = BgdSlope*(x-BgFrom) + BgdOff;
        f_tot += BGd;
    }

//    cout<<NSubPeaks<<" "<<BgFrom<<" "<<BgTo<<" "<<BgdOff<<" "<<BgdSlope<<endl;
    for(int i=0 ; i<NSubPeaks ; i++)
    {
        double Ampli     = pp[5+i*6+0];
        double Mean      = pp[5+i*6+1];
        double Sigma     = pp[5+i*6+2]*1./sqrt(8.*log(2.));;
        double Lambda    = pp[5+i*6+3];
        double Rho       = pp[5+i*6+4];
        double S         = pp[5+i*6+5];

//        cout<<Ampli<<" "<<Mean<<" "<<Sigma<<" "<<Lambda<<" "<<Rho<<" "<<S<<endl;

        double U         = (x-Mean)/Sigma;
        double f_g       = Ampli*TMath::Exp(-U*U*0.5);
        double f_lambda  = Ampli*TMath::Exp(-0.5*Lambda*(2.*U-Lambda));
        double f_rho     = Ampli*TMath::Exp(-0.5*Rho*(2.*U-Rho));
        double f_S       = Ampli*S*1./((1+TMath::Exp(U))*(1+TMath::Exp(U)));

        if(U<Lambda) f_tot += f_lambda;
        else if(U>Rho) f_tot += f_rho;
        else f_tot += f_g;

        f_tot += f_S;
    }

    return f_tot;
}


double FitFunctions::Gaus(double *xx, double *par)
{
    double x = xx[0];
    double N = par[0];
    double mean = par[1];
    double sigma = par[2];

    double f = N*exp(-0.5*((x-mean)/sigma)*((x-mean)/sigma));

    return f;
}


TF1 *FitFunctions::GetDinoFct(TString Name,double min, double max, int Npar)
{
    TF1 *f = new TF1(Name,this,&FitFunctions::DinoFct,min,max,Npar,"FitFunctions","DinoFct");

    return f;
}