CrystalBall.hh
00001
00002
00003
00004
00005
00006 #include <math.h>
00007 #include "TMath.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 Double_t CB(Double_t *x, Double_t *par);
00020 Double_t CBIntegral(Double_t mean, Double_t sigma, Double_t alpha, Double_t n, Double_t min, Double_t max);
00021 Double_t ApproxErf(Double_t arg);
00022
00023 Double_t CB(Double_t *x, Double_t *par)
00024 {
00025 Double_t cb_result;
00026 Double_t t = (x[0] - par[0])/par[1];
00027 if(par[2] < 0) t = -t;
00028
00029 Double_t absAlpha = fabs((Double_t)par[2]);
00030
00031 if (t >= -absAlpha)
00032 {
00033 cb_result = exp(-0.5*t*t);
00034 }
00035 else
00036 {
00037 Double_t a = TMath::Power(par[3]/absAlpha,par[3]);
00038 a *= exp(-0.5*absAlpha*absAlpha);
00039 Double_t b= par[3]/absAlpha - absAlpha;
00040
00041 cb_result = a/TMath::Power(b - t, par[3]);
00042 }
00043
00044 Double_t N = 1./CBIntegral(par[0], par[1], par[2], par[3], par[4], par[5]);
00045 cb_result *= N;
00046
00047 cb_result = par[6]/(par[5]-par[4]) + (1-par[6])*cb_result;
00048
00049 return cb_result;
00050 }
00051
00052
00053 Double_t CBIntegral(Double_t mean, Double_t sigma, Double_t alpha, Double_t n, Double_t min, Double_t max)
00054 {
00055 static const double sqrtPiOver2 = 1.2533141373;
00056 static const double sqrt2 = 1.4142135624;
00057
00058 double result = 0.0;
00059 bool useLog = false;
00060
00061 if( fabs(n-1.0) < 1.0e-05 )
00062 useLog = true;
00063
00064 double sig = fabs((Double_t)sigma);
00065
00066 double tmin = (min-mean)/sig;
00067 double tmax = (max-mean)/sig;
00068
00069 if(alpha < 0) {
00070 double tmp = tmin;
00071 tmin = -tmax;
00072 tmax = -tmp;
00073 }
00074
00075 double absAlpha = fabs((Double_t)alpha);
00076
00077 if( tmin >= -absAlpha ) {
00078 result += sig*sqrtPiOver2*( ApproxErf(tmax/sqrt2)
00079 - ApproxErf(tmin/sqrt2) );
00080 }
00081 else if( tmax <= -absAlpha ) {
00082 double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
00083 double b = n/absAlpha - absAlpha;
00084
00085 if(useLog) {
00086 result += a*sig*( log(b-tmin) - log(b-tmax) );
00087 }
00088 else {
00089 result += a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
00090 - 1.0/(TMath::Power(b-tmax,n-1.0)) );
00091 }
00092 }
00093 else {
00094 double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
00095 double b = n/absAlpha - absAlpha;
00096
00097 double term1 = 0.0;
00098 if(useLog) {
00099 term1 = a*sig*( log(b-tmin) - log(n/absAlpha));
00100 }
00101 else {
00102 term1 = a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
00103 - 1.0/(TMath::Power(n/absAlpha,n-1.0)) );
00104 }
00105
00106 double term2 = sig*sqrtPiOver2*( ApproxErf(tmax/sqrt2)
00107 - ApproxErf(-absAlpha/sqrt2) );
00108
00109
00110 result += term1 + term2;
00111 }
00112
00113 return result;
00114 }
00115
00116 Double_t ApproxErf(Double_t arg)
00117 {
00118 static const double erflim = 5.0;
00119 if( arg > erflim )
00120 return 1.0;
00121 if( arg < -erflim )
00122 return -1.0;
00123
00124 return TMath::Erf(arg);
00125 }
00126