Line data    Source code

       1             : // $Id: EvtDiLog.cpp,v 1.3 2009-03-16 15:52:42 robbep 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           0 :   if ( x == 1. ) h = PI6 ;
      19           0 :   else if ( x == -1. ) h = -PI12 ;
      20             :   else {
      21           0 :     t = -x ;
      22           0 :     if ( t <= -2. ) {
      23           0 :       y = -1./(1.+t) ;
      24             :       s = 1. ;
      25           0 :       a = -PI3 + HF * ( std::pow( log(-t) , 2 ) - 
      26           0 :                         std::pow( log( 1. + 1./t ) , 2 ) ) ;
      27           0 :     } else if ( t < -1. ) {
      28           0 :       y = -1. - t ;
      29             :       s = -1. ;
      30           0 :       a = log( -t ) ;
      31           0 :       a = -PI6 + a * ( a + log( 1. + 1./t ) ) ;
      32           0 :     } else if ( t <= -HF ) {
      33           0 :       y = - (1. + t ) / t ;
      34             :       s = 1. ;
      35           0 :       a = log( -t ) ;
      36           0 :       a = -PI6 + a * ( -HF * a + log( 1. + t ) ) ;
      37           0 :     } else if ( t < 0 ) {
      38           0 :       y = -t / ( 1. + t ) ;
      39             :       s = -1. ;
      40           0 :       a = HF * std::pow( log( 1. + t ) , 2 ) ;
      41           0 :     } else if ( t <= 1. ) {
      42             :       y = t ;
      43             :       s = 1. ;
      44             :       a = 0. ;
      45           0 :     } else {
      46           0 :       y = 1. / t ;
      47             :       s = -1. ;
      48           0 :       a = PI6 + HF * std::pow( log( t ) , 2 ) ;
      49             :     }
      50             :     
      51           0 :     h = y + y - 1. ;
      52           0 :     alfa = h + h ;
      53             :     b1 = 0. ;
      54             :     b2 = 0. ;
      55           0 :     for ( int i = 19 ; i >= 0 ; --i ) {
      56           0 :       b0 = C[ i ] + alfa * b1 - b2 ;
      57             :       b2 = b1 ;
      58             :       b1 = b0 ;
      59             :     }
      60             :     
      61           0 :     h = -(s * ( b0 -h * b2 ) + a ) ;
      62             :   }
      63             :   
      64           0 :   return h ;
      65             : }