TEvtGen/EvtGenBase/EvtDiLog.cxx

   1 // $Id: EvtDiLog.cc,v 1.1 2009/02/18 21:21:33 ryd Exp $
   2 // Include files
   3
   4 // local
   5 #include "EvtGenBase/EvtDiLog.hh"
   6
   7 #include <cmath>
   8
   9 //-----------------------------------------------------------------------------
  10 // Implementation file for class : EvtDiLog
  11 //
  12 // 2007-01-23 : Patrick Robbe
  13 //-----------------------------------------------------------------------------
  14
  15 double EvtDiLog::DiLog( double x ) {
  16
  17   double h , t , y , s , a , alfa , b0, b1, b2 ;
  18   if ( x == 1. ) h = PI6 ;
  19   else if ( x == -1. ) h = -PI12 ;
  20   else {
  21     t = -x ;
  22     if ( t <= -2. ) {
  23       y = -1./(1.+t) ;
  24       s = 1. ;
  25       a = -PI3 + HF * ( std::pow( log(-t) , 2 ) -
  26                         std::pow( log( 1. + 1./t ) , 2 ) ) ;
  27     } else if ( t < -1. ) {
  28       y = -1. - t ;
  29       s = -1. ;
  30       a = log( -t ) ;
  31       a = -PI6 + a * ( a + log( 1. + 1./t ) ) ;
  32     } else if ( t <= -HF ) {
  33       y = - (1. + t ) / t ;
  34       s = 1. ;
  35       a = log( -t ) ;
  36       a = -PI6 + a * ( -HF * a + log( 1. + t ) ) ;
  37     } else if ( t < 0 ) {
  38       y = -t / ( 1. + t ) ;
  39       s = -1. ;
  40       a = HF * std::pow( log( 1. + t ) , 2 ) ;
  41     } else if ( t <= 1. ) {
  42       y = t ;
  43       s = 1. ;
  44       a = 0. ;
  45     } else {
  46       y = 1. / t ;
  47       s = -1. ;
  48       a = PI6 + HF * std::pow( log( t ) , 2 ) ;
  49     }
  50
  51     h = y + y - 1. ;
  52     alfa = h + h ;
  53     b1 = 0. ;
  54     b2 = 0. ;
  55     for ( int i = 19 ; i >= 0 ; --i ) {
  56       b0 = C[ i ] + alfa * b1 - b2 ;
  57       b2 = b1 ;
  58       b1 = b0 ;
  59     }
  60
  61     h = -(s * ( b0 -h * b2 ) + a ) ;
  62   }
  63
  64   return h ;
  65 }