FitFunctions.cpp

#include "FitFunctions.h"
#include "TMath.h"
#include <cmath>
#include "CachedFunction1D.h"

using namespace std;

namespace SctData {
    
    //Private function to hold the cache and to make sure 
    //the cache is only created if necessary
    CachedFunction1D& getErfCache() throw() {
        static CachedFunction1D cachedErf(-4., 4., 0.0001, &::erf );
        return cachedErf;
    }
    
    Double_t erf(double *x, double *par) throw() {
    double xx, y;
    
    // if width==0, return zero:
    if(par[2] == 0){
        xx = 0.;
    }else{
        xx = ( (par[1]-x[0]) / par[2]) / 1.414213562; 
    }
    
    //y  = (par[0]/2.) * ( 1 - TMath::Erf(xx) );
    //y  = (par[0]/2.) * ( 1 - ::erf(xx) );
    y = (par[0]/2.) * ( 1 - getErfCache().eval(xx)) ;
    return y;
    }
    
    Double_t erfc(double *x, double *par) throw() {
        //static CachedFunction1D cachedErf(-4., 4., 0.0001, &::erf );
    double xx, y;
    
    // if width==0, return zero:
    if(par[2] == 0.){
        xx = 0.;
    }else{
        xx = ( (par[1]-x[0]) / par[2]) / 1.414213562; 
    }
    
    //y  = (par[0]/2.) * ( 1 + TMath::Erf(xx) );
    //y  = (par[0]/2.) * ( 1 + ::erf(xx) );
    y  = (par[0]/2.) * ( 1 + getErfCache().eval(xx));
    return y;
    }
    
    Double_t tophat(double *x, double *par) throw() {
        //static CachedFunction1D cachedErf(-4., 4., 0.0001, &::erf );
    double trailPart = 2;
    double leadPart = 2;
    
    if (par[2]<0) par[2]=-par[2];
    if (par[4]<0) par[4]=-par[4];
    
    if (par[2]>1e-100) {
        double rise = ( (par[1]-x[0]) / par[2]) / 1.414213562;
        leadPart = 1 - getErfCache().eval(rise);
    } else {
        if (x[0] < par[1]) leadPart = 0;
    }
    if (par[4]>1e-100) {
        double fall = ( (par[3]-x[0]) / par[4]) / 1.414213562;
        trailPart = 1 + getErfCache().eval(fall);
    } else {
        if (x[0] > par[3]) trailPart = 0;
    }   
        
    //   double y=(par[0]/4.)*( 1 - TMath::Erf(rise) ) * ( 1 + TMath::Erf(fall) );
    //   double y=(par[0]/4.)*( 1 - erf(rise) ) * ( 1 + erf(fall) );
    double y=(par[0]/4.)*leadPart*trailPart;
    return y;
    }
}

---