Fixes for report #63583: High CPU time spent in TMath::Erf

diff --git a/ITS/AliITSsimulationSPD.cxx b/ITS/AliITSsimulationSPD.cxx
index 428a2a0fd249fdaf5f7c3cee0655a4b3bddbf260..076f11b2b43227cb33022fd98ce8e193b06022d0 100644 (file)
--- a/ITS/AliITSsimulationSPD.cxx
+++ b/ITS/AliITSsimulationSPD.cxx
@@ -32,6 +32,7 @@ $Id$
 #include "AliLog.h"
 #include "AliRun.h"
 #include "AliMagF.h"
+#include "AliMathBase.h"
 
 //#define DEBUG
 
@@ -627,13 +628,13 @@ void AliITSsimulationSPD::SpreadCharge(Double_t x0,Double_t z0,
        z1 -= z0; // Distance from where track traveled
        z2 -= z0; // Distance from where track traveled
        s   = 0.25; // Correction based on definision of Erfc
-       s  *= TMath::Erfc(sp*x1) - TMath::Erfc(sp*x2);
+       s  *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2);
        if(GetDebug(3)){
            cout <<"el="<<el<<" ix0="<<ix0<<" ix="<<ix<<" x0="<<x<<
                " iz0="<<iz0<<" iz="<<iz<<" z0="<<z<< 
                " sp*x1="<<sp*x1<<" sp*x2="<<sp*x2<<" s="<<s;
        } // end if GetDebug
-       s  *= TMath::Erfc(sp*z1) - TMath::Erfc(sp*z2);
+       s  *= AliMathBase::ErfcFast(sp*z1) - AliMathBase::ErfcFast(sp*z2);
        if(GetDebug(3)){
            cout<<" sp*z1="<<sp*z1<<" sp*z2="<<sp*z2<<" s="<<s<< endl;
        } // end if GetDebug
@@ -707,13 +708,13 @@ void AliITSsimulationSPD::SpreadChargeAsym(Double_t x0,Double_t z0,
        z1 -= z0; // Distance from where track traveled
        z2 -= z0; // Distance from where track traveled
        s   = 0.25; // Correction based on definision of Erfc
-       s  *= TMath::Erfc(spx*x1) - TMath::Erfc(spx*x2);
+       s  *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2);
        if(GetDebug(3)){
            cout <<"el="<<el<<" ix0="<<ix0<<" ix="<<ix<<" x0="<<x<<
                " iz0="<<iz0<<" iz="<<iz<<" z0="<<z<< 
                " spx*x1="<<spx*x1<<" spx*x2="<<spx*x2<<" s="<<s;
        } // end if GetDebug
-       s  *= TMath::Erfc(spz*z1) - TMath::Erfc(spz*z2);
+       s  *= AliMathBase::ErfcFast(spz*z1) - AliMathBase::ErfcFast(spz*z2);
        if(GetDebug(3)){
            cout<<" spz*z1="<<spz*z1<<" spz*z2="<<spz*z2<<" s="<<s<< endl;
        } // end if GetDebug

diff --git a/ITS/AliITSsimulationSSD.cxx b/ITS/AliITSsimulationSSD.cxx
index e4293e2ed887065233ca0bf282094f9a4607e0f7..ae473468b46aa202dde3bc4f231e472ff77d7b88 100644 (file)
--- a/ITS/AliITSsimulationSSD.cxx
+++ b/ITS/AliITSsimulationSSD.cxx
@@ -35,7 +35,7 @@
 #include "AliITSsimulationSSD.h"
 #include "AliITSTableSSD.h"
 #include <TF1.h>
-
+#include "AliMathBase.h"
 
 ClassImp(AliITSsimulationSSD)
 ////////////////////////////////////////////////////////////////////////
@@ -518,7 +518,7 @@ Float_t AliITSsimulationSSD::F(Float_t av, Float_t x, Float_t s) {
     Float_t sigm2 = sqrt2*s;
     Float_t integral;
 
-    integral = 0.5 * TMath::Erf( (x - av) / sigm2);
+    integral = 0.5 * AliMathBase::ErfFast( (x - av) / sigm2);
     return integral;
 }
 //______________________________________________________________________

diff --git a/ITS/AliITStrackerMI.cxx b/ITS/AliITStrackerMI.cxx
index aa82b3d6f9ad8c2d6ef1f18e1e8a0aa1a69040d9..6fc40d5a116f8c6bdb3753ae44951ef8f2bc2288 100644 (file)
--- a/ITS/AliITStrackerMI.cxx
+++ b/ITS/AliITStrackerMI.cxx
@@ -61,6 +61,7 @@
 #include "AliITSPlaneEffSSD.h"
 #include "AliITSV0Finder.h"
 #include "AliITStrackerMI.h"
+#include "AliMathBase.h"
 
 ClassImp(AliITStrackerMI)
 
@@ -2593,8 +2594,8 @@ Double_t  AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zer
   }
   // probability that the true z is in the range [zmin,zmax] (i.e. inside 
   // dead zone)
-  probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) - 
-                     TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
+  probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) - 
+                     AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
   return probability;
 }
 //------------------------------------------------------------------------

diff --git a/STEER/AliMathBase.cxx b/STEER/AliMathBase.cxx
index 21bee1a46d9548c1e36a771841154a1927ca7501..cdae2369d06d0477342b06bbae8e26c2f9a8c7d2 100644 (file)
--- a/STEER/AliMathBase.cxx
+++ b/STEER/AliMathBase.cxx
@@ -554,6 +554,19 @@ Float_t AliMathBase::GetCOG(Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMa
 }
 
 
+Double_t AliMathBase::ErfcFast(Double_t x){
+  // Fast implementation of the complementary error function
+  // The error of the approximation is |eps(x)| < 5E-4
+  // See Abramowitz and Stegun, p.299, 7.1.27
+
+  Double_t z = TMath::Abs(x);
+  Double_t ans = 1+z*(0.278393+z*(0.230389+z*(0.000972+z*0.078108)));
+  ans = 1.0/ans;
+  ans *= ans;
+  ans *= ans;
+
+  return (x>=0.0 ? ans : 2.0 - ans);
+}
 
 ///////////////////////////////////////////////////////////////
 //////////////         TEST functions /////////////////////////

diff --git a/STEER/AliMathBase.h b/STEER/AliMathBase.h
index 3b16daa8991015a0ea5f1bd36ddd07aa3a9a1780..fdbec7d01075fb657985aac11ca2043ed2ba9444 100644 (file)
--- a/STEER/AliMathBase.h
+++ b/STEER/AliMathBase.h
@@ -35,6 +35,9 @@ class AliMathBase : public TObject
   static TGraph2D *  MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type);
   static TGraph *  MakeStat1D(TH3 * his, Int_t delta1, Int_t type);
 
+  static Double_t ErfcFast(Double_t x);                           // Complementary error function erfc(x)
+  static Double_t ErfFast(Double_t x) {return 1-ErfcFast(x);}     // Error function erf(x)
+
   //
   // TestFunctions:
   //

diff --git a/STEER/AliSignalProcesor.cxx b/STEER/AliSignalProcesor.cxx
index 782fc2aa9a36f112ab8ee9cbad8f3ccef6fb90b7..bcc6667af422453fe75a701a325341aee1c070d5 100644 (file)
--- a/STEER/AliSignalProcesor.cxx
+++ b/STEER/AliSignalProcesor.cxx
@@ -19,6 +19,7 @@
 #include <TMath.h>
 
 #include "AliSignalProcesor.h"
+#include "AliMathBase.h"
 
 
 ClassImp(AliSignalProcesor)
@@ -48,10 +49,10 @@ Double_t asymgauss(Double_t* x, Double_t* par)
   if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow 
   if (TMath::Abs(par[4])>1) return 0;
   Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
-    *(1-TMath::Erf((par3save*sigma2-dx)/(sqrt2*par2save)));
+    *(1-AliMathBase::ErfFast((par3save*sigma2-dx)/(sqrt2*par2save)));
 
   Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
-    *(1-TMath::Erf((par5save*sigma2-dx)/(sqrt2*par2save)));
+    *(1-AliMathBase::ErfFast((par5save*sigma2-dx)/(sqrt2*par2save)));
 
 
   return par[0]*(exp1+par[4]*exp2);
@@ -82,10 +83,10 @@ Double_t asymgaussN(Double_t* x, Double_t* par)
   if (TMath::Abs(par[4])>=1) return 0;
   
   Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
-    *0.5*(1-TMath::Erf((par3save*sigma2-dx)/(sqrt2*par2save)));
+    *0.5*(1-AliMathBase::ErfFast((par3save*sigma2-dx)/(sqrt2*par2save)));
 
   Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
-    *0.5*(1-TMath::Erf((par5save*sigma2-dx)/(sqrt2*par2save)));
+    *0.5*(1-AliMathBase::ErfFast((par5save*sigma2-dx)/(sqrt2*par2save)));
 
 
   return par[0]*(1.*exp1+par[4]*exp2)/(1.+par[4]);

diff --git a/TOF/AliTOFtrackerMI.cxx b/TOF/AliTOFtrackerMI.cxx
index bd61bd18ab02cdece379f30172470f98dec92e75..8e9aa1f5c030e9f16d538c42092f0a764ed3a4b7 100644 (file)
--- a/TOF/AliTOFtrackerMI.cxx
+++ b/TOF/AliTOFtrackerMI.cxx
@@ -38,6 +38,8 @@
 #include "AliTOFtrackerMI.h"
 #include "AliTOFtrack.h"
 
+#include "AliMathBase.h"
+
 class TGeoManager;
 
 extern TGeoManager *gGeoManager;
@@ -806,28 +808,28 @@ void AliTOFtrackerMI::GetLikelihood(Float_t dy, Float_t dz, const Double_t *cov,
   //
   normwidth = fDy/sigmay;
   normd     = dy/sigmay;
-  p0 = 0.5*(1+TMath::Erf(normd-normwidth*0.5));
-  p1 = 0.5*(1+TMath::Erf(normd+normwidth*0.5));  
+  p0 = 0.5*(1+AliMathBase::ErfFast(normd-normwidth*0.5));
+  p1 = 0.5*(1+AliMathBase::ErfFast(normd+normwidth*0.5));  
   py+= 0.75*(p1-p0);
   //
   normwidth = fDy/(3.*sigmay);
   normd     = dy/(3.*sigmay);
-  p0 = 0.5*(1+TMath::Erf(normd-normwidth*0.5));
-  p1 = 0.5*(1+TMath::Erf(normd+normwidth*0.5));  
+  p0 = 0.5*(1+AliMathBase::ErfFast(normd-normwidth*0.5));
+  p1 = 0.5*(1+AliMathBase::ErfFast(normd+normwidth*0.5));  
   py+= 0.25*(p1-p0);
   // 
   // pz calculation   - 75% admixture of original sigma - 25% tails 
   //
   normwidth = fDz/sigmaz;
   normd     = dz/sigmaz;
-  p0 = 0.5*(1+TMath::Erf(normd-normwidth*0.5));
-  p1 = 0.5*(1+TMath::Erf(normd+normwidth*0.5));  
+  p0 = 0.5*(1+AliMathBase::ErfFast(normd-normwidth*0.5));
+  p1 = 0.5*(1+AliMathBase::ErfFast(normd+normwidth*0.5));  
   pz+= 0.75*(p1-p0);
   //
   normwidth = fDz/(3.*sigmaz);
   normd     = dz/(3.*sigmaz);
-  p0 = 0.5*(1+TMath::Erf(normd-normwidth*0.5));
-  p1 = 0.5*(1+TMath::Erf(normd+normwidth*0.5));  
+  p0 = 0.5*(1+AliMathBase::ErfFast(normd-normwidth*0.5));
+  p1 = 0.5*(1+AliMathBase::ErfFast(normd+normwidth*0.5));  
   pz+= 0.25*(p1-p0);
 }
 //_________________________________________________________________________

diff --git a/TPC/AliTPCClusterParam.cxx b/TPC/AliTPCClusterParam.cxx
index 61bd5911bb447460f54717d91a3f1f3c2709a100..08dd92591b224b1d3980154e64776a6c93f4d530 100644 (file)
--- a/TPC/AliTPCClusterParam.cxx
+++ b/TPC/AliTPCClusterParam.cxx
@@ -98,6 +98,8 @@ AliTPCClusterParam::SetInstance(param);
 #include "AliTPCcalibDB.h"
 #include "AliTPCParam.h"
 
+#include "AliMathBase.h"
+
 ClassImp(AliTPCClusterParam)
 
 
@@ -1605,8 +1607,8 @@ Double_t  AliTPCClusterParam::GaussConvolution(Double_t x0, Double_t x1, Double_
   Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
   //
   Double_t sigmaErf =  2*s0*s1*TMath::Sqrt(2*sigma2);                       
-  Double_t erf0 = TMath::Erf( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
-  Double_t erf1 = TMath::Erf( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
+  Double_t erf0 = AliMathBase::ErfFast( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
+  Double_t erf1 = AliMathBase::ErfFast( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
   Double_t norm = 1./TMath::Sqrt(sigma2);
   norm/=2.*TMath::Sqrt(2.*TMath::Pi());
   Double_t val  = norm*exp0*(erf0+erf1);

diff --git a/TRD/AliTRDCalibraFit.cxx b/TRD/AliTRDCalibraFit.cxx
index 377b9fbd7759aede96a6de508f3ced43bbd774ae..e922d841b349994b000ddc205450c81fc55f977f 100644 (file)
--- a/TRD/AliTRDCalibraFit.cxx
+++ b/TRD/AliTRDCalibraFit.cxx
@@ -5971,9 +5971,9 @@ Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
   Double_t  sigma2  = par2save*par2save;
   Double_t  sqrt2   = TMath::Sqrt(2.0);
   Double_t  exp1    = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
-                               * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
+                               * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
   Double_t  exp2    = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
-                               * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
+                               * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
 
   //return par[0]*(exp1+par[4]*exp2);
   return par[0] * (exp1 + 1.00124 * exp2);