CrystalBall.hh

00001 
00002 // Fit function for "crystal ball" lineshapes
00003 // Code adapted from RooFit: RooCBShape.cc, v 1.15
00004 
00005 
00006 #include <math.h>
00007 #include "TMath.h"
00008 
00009 // Here are the parameters for the fit:
00010 // x[0] = x
00011 // par[0] = mean
00012 // par[1] = sigma
00013 // par[2] = alpha
00014 // par[3] = n (order)
00015 // par[4] = range.min
00016 // par[5] = range.max 
00017 // par[6] = f_bkg (constant background term)
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         // Normalization constant
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 

Generated on Tue Nov 16 10:49:57 2010 for CUORE Software by  doxygen 1.5.6