Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/HEPNumerics/include/G4Integrator.icc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /global/HEPNumerics/include/G4Integrator.icc (Version 11.3.0) and /global/HEPNumerics/include/G4Integrator.icc (Version 7.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 // G4Integrator inline methods implementation  << 
 27 //                                                 23 //
 28 // Author: V.Grichine, 04.09.1999 - First impl <<  24 // $Id: G4Integrator.icc,v 1.11 2004/11/12 17:38:29 gcosmo Exp $
 29 //         G4SimpleIntegration class with H.P. <<  25 // GEANT4 tag $Name: geant4-07-00-patch-01 $
 30 //         E.TCherniaev advises                <<  26 //
 31 // ------------------------------------------- <<  27 // Implementation of G4Integrator methods. 
                                                   >>  28 //
                                                   >>  29 // 
 32                                                    30 
 33 //////////////////////////////////////////////     31 /////////////////////////////////////////////////////////////////////
 34 //                                                 32 //
 35 // Sympson integration method                      33 // Sympson integration method
 36 //                                                 34 //
 37 //////////////////////////////////////////////     35 /////////////////////////////////////////////////////////////////////
 38 //                                                 36 //
 39 // Integration of class member functions T::f  <<  37 // Integration of class member functions T::f by Simpson method. 
 40                                                << 
 41 template <class T, class F>                    << 
 42 G4double G4Integrator<T, F>::Simpson(T& typeT, << 
 43                                      G4double  << 
 44 {                                              << 
 45   G4int i;                                     << 
 46   G4double step  = (xFinal - xInitial) / itera << 
 47   G4double x     = xInitial;                   << 
 48   G4double xPlus = xInitial + 0.5 * step;      << 
 49   G4double mean  = ((typeT.*f)(xInitial) + (ty << 
 50   G4double sum   = (typeT.*f)(xPlus);          << 
 51                                                    38 
 52   for(i = 1; i < iterationNumber; ++i)         <<  39 template <class T, class F> 
 53   {                                            <<  40 G4double G4Integrator<T,F>::Simpson( T&       typeT, 
 54     x += step;                                 <<  41                                      F        f,
 55     xPlus += step;                             <<  42                                      G4double xInitial,
 56     mean += (typeT.*f)(x);                     <<  43                                      G4double xFinal,
 57     sum += (typeT.*f)(xPlus);                  <<  44                                      G4int    iterationNumber ) 
 58   }                                            <<  45 {
 59   mean += 2.0 * sum;                           <<  46    G4int    i ;
                                                   >>  47    G4double step = (xFinal - xInitial)/iterationNumber ;
                                                   >>  48    G4double x = xInitial ;
                                                   >>  49    G4double xPlus = xInitial + 0.5*step ;
                                                   >>  50    G4double mean = ( (typeT.*f)(xInitial) + (typeT.*f)(xFinal) )*0.5 ;
                                                   >>  51    G4double sum = (typeT.*f)(xPlus) ;
                                                   >>  52 
                                                   >>  53    for(i=1;i<iterationNumber;i++)
                                                   >>  54    {
                                                   >>  55       x     += step ;
                                                   >>  56       xPlus += step ;
                                                   >>  57       mean  += (typeT.*f)(x) ;
                                                   >>  58       sum   += (typeT.*f)(xPlus) ;
                                                   >>  59    }
                                                   >>  60    mean += 2.0*sum ;
 60                                                    61 
 61   return mean * step / 3.0;                    <<  62    return mean*step/3.0 ;   
 62 }                                                  63 }
 63                                                    64 
 64 //////////////////////////////////////////////     65 /////////////////////////////////////////////////////////////////////
 65 //                                                 66 //
 66 // Integration of class member functions T::f      67 // Integration of class member functions T::f by Simpson method.
 67 // Convenient to use with 'this' pointer           68 // Convenient to use with 'this' pointer
 68                                                    69 
 69 template <class T, class F>                    <<  70 template <class T, class F> 
 70 G4double G4Integrator<T, F>::Simpson(T* ptrT,  <<  71 G4double G4Integrator<T,F>::Simpson( T*       ptrT, 
 71                                      G4double  <<  72                                 F        f,
 72 {                                              <<  73                                 G4double xInitial,
 73   G4int i;                                     <<  74                                 G4double xFinal,
 74   G4double step  = (xFinal - xInitial) / itera <<  75                                 G4int    iterationNumber ) 
 75   G4double x     = xInitial;                   <<  76 {
 76   G4double xPlus = xInitial + 0.5 * step;      <<  77    G4int    i ;
 77   G4double mean  = ((ptrT->*f)(xInitial) + (pt <<  78    G4double step = (xFinal - xInitial)/iterationNumber ;
 78   G4double sum   = (ptrT->*f)(xPlus);          <<  79    G4double x = xInitial ;
                                                   >>  80    G4double xPlus = xInitial + 0.5*step ;
                                                   >>  81    G4double mean = ( (ptrT->*f)(xInitial) + (ptrT->*f)(xFinal) )*0.5 ;
                                                   >>  82    G4double sum = (ptrT->*f)(xPlus) ;
                                                   >>  83 
                                                   >>  84    for(i=1;i<iterationNumber;i++)
                                                   >>  85    {
                                                   >>  86       x     += step ;
                                                   >>  87       xPlus += step ;
                                                   >>  88       mean  += (ptrT->*f)(x) ;
                                                   >>  89       sum   += (ptrT->*f)(xPlus) ;
                                                   >>  90    }
                                                   >>  91    mean += 2.0*sum ;
 79                                                    92 
 80   for(i = 1; i < iterationNumber; ++i)         <<  93    return mean*step/3.0 ;   
 81   {                                            << 
 82     x += step;                                 << 
 83     xPlus += step;                             << 
 84     mean += (ptrT->*f)(x);                     << 
 85     sum += (ptrT->*f)(xPlus);                  << 
 86   }                                            << 
 87   mean += 2.0 * sum;                           << 
 88                                                << 
 89   return mean * step / 3.0;                    << 
 90 }                                                  94 }
 91                                                    95 
 92 //////////////////////////////////////////////     96 /////////////////////////////////////////////////////////////////////
 93 //                                                 97 //
 94 // Integration of class member functions T::f      98 // Integration of class member functions T::f by Simpson method.
 95 // Convenient to use, when function f is defin     99 // Convenient to use, when function f is defined in global scope, i.e. in main()
 96 // program                                        100 // program
 97                                                   101 
 98 template <class T, class F>                    << 102 template <class T, class F> 
 99 G4double G4Integrator<T, F>::Simpson(G4double  << 103 G4double G4Integrator<T,F>::Simpson( G4double (*f)(G4double),
100                                      G4double  << 104                                 G4double xInitial,
101 {                                              << 105                                 G4double xFinal,
102   G4int i;                                     << 106                                 G4int    iterationNumber ) 
103   G4double step  = (xFinal - xInitial) / itera << 107 {
104   G4double x     = xInitial;                   << 108    G4int    i ;
105   G4double xPlus = xInitial + 0.5 * step;      << 109    G4double step = (xFinal - xInitial)/iterationNumber ;
106   G4double mean  = ((*f)(xInitial) + (*f)(xFin << 110    G4double x = xInitial ;
107   G4double sum   = (*f)(xPlus);                << 111    G4double xPlus = xInitial + 0.5*step ;
                                                   >> 112    G4double mean = ( (*f)(xInitial) + (*f)(xFinal) )*0.5 ;
                                                   >> 113    G4double sum = (*f)(xPlus) ;
                                                   >> 114 
                                                   >> 115    for(i=1;i<iterationNumber;i++)
                                                   >> 116    {
                                                   >> 117       x     += step ;
                                                   >> 118       xPlus += step ;
                                                   >> 119       mean  += (*f)(x) ;
                                                   >> 120       sum   += (*f)(xPlus) ;
                                                   >> 121    }
                                                   >> 122    mean += 2.0*sum ;
108                                                   123 
109   for(i = 1; i < iterationNumber; ++i)         << 124    return mean*step/3.0 ;   
110   {                                            << 
111     x += step;                                 << 
112     xPlus += step;                             << 
113     mean += (*f)(x);                           << 
114     sum += (*f)(xPlus);                        << 
115   }                                            << 
116   mean += 2.0 * sum;                           << 
117                                                << 
118   return mean * step / 3.0;                    << 
119 }                                                 125 }
120                                                   126 
121 //////////////////////////////////////////////    127 //////////////////////////////////////////////////////////////////////////
122 //                                                128 //
123 // Adaptive Gauss method                          129 // Adaptive Gauss method
124 //                                                130 //
125 //////////////////////////////////////////////    131 //////////////////////////////////////////////////////////////////////////
126 //                                                132 //
127 //                                                133 //
128                                                   134 
129 template <class T, class F>                    << 135 template <class T, class F> 
130 G4double G4Integrator<T, F>::Gauss(T& typeT, F << 136 G4double G4Integrator<T,F>::Gauss( T& typeT, F f,
131                                    G4double xF << 137                                    G4double xInitial, G4double xFinal   ) 
132 {                                              << 138 {
133   static const G4double root = 1.0 / std::sqrt << 139    static G4double root = 1.0/std::sqrt(3.0) ;
134                                                << 140    
135   G4double xMean = (xInitial + xFinal) / 2.0;  << 141    G4double xMean = (xInitial + xFinal)/2.0 ;
136   G4double Step  = (xFinal - xInitial) / 2.0;  << 142    G4double Step = (xFinal - xInitial)/2.0 ;
137   G4double delta = Step * root;                << 143    G4double delta = Step*root ;
138   G4double sum   = ((typeT.*f)(xMean + delta)  << 144    G4double sum = ((typeT.*f)(xMean + delta) + 
139                                                << 145                    (typeT.*f)(xMean - delta)) ;
140   return sum * Step;                           << 146    
                                                   >> 147    return sum*Step ;   
141 }                                                 148 }
142                                                   149 
143 //////////////////////////////////////////////    150 //////////////////////////////////////////////////////////////////////
144 //                                                151 //
145 //                                                152 //
146                                                   153 
147 template <class T, class F>                    << 154 template <class T, class F> G4double 
148 G4double G4Integrator<T, F>::Gauss(T* ptrT, F  << 155 G4Integrator<T,F>::Gauss( T* ptrT, F f, G4double a, G4double b )
149 {                                                 156 {
150   return Gauss(*ptrT, f, a, b);                << 157   return Gauss(*ptrT,f,a,b) ;
151 }                                                 158 }
152                                                   159 
153 //////////////////////////////////////////////    160 ///////////////////////////////////////////////////////////////////////
154 //                                                161 //
155 //                                                162 //
156                                                   163 
157 template <class T, class F>                       164 template <class T, class F>
158 G4double G4Integrator<T, F>::Gauss(G4double (* << 165 G4double G4Integrator<T,F>::Gauss( G4double (*f)(G4double), 
159                                    G4double xF << 166                               G4double xInitial, G4double xFinal) 
160 {                                                 167 {
161   static const G4double root = 1.0 / std::sqrt << 168    static G4double root = 1.0/std::sqrt(3.0) ;
162                                                << 169    
163   G4double xMean = (xInitial + xFinal) / 2.0;  << 170    G4double xMean = (xInitial + xFinal)/2.0 ;
164   G4double Step  = (xFinal - xInitial) / 2.0;  << 171    G4double Step  = (xFinal - xInitial)/2.0 ;
165   G4double delta = Step * root;                << 172    G4double delta = Step*root ;
166   G4double sum   = ((*f)(xMean + delta) + (*f) << 173    G4double sum   = ( (*f)(xMean + delta) + (*f)(xMean - delta) ) ;
167                                                << 174    
168   return sum * Step;                           << 175    return sum*Step ;   
169 }                                                 176 }
170                                                   177 
171 //////////////////////////////////////////////    178 ///////////////////////////////////////////////////////////////////////////
172 //                                                179 //
173 //                                                180 //
174                                                   181 
175 template <class T, class F>                    << 182 template <class T, class F>  
176 void G4Integrator<T, F>::AdaptGauss(T& typeT,  << 183 void G4Integrator<T,F>::AdaptGauss( T& typeT, F f, G4double  xInitial,
177                                     G4double x << 184                                G4double  xFinal, G4double fTolerance,
178                                     G4double&  << 185                                G4double& sum,
179 {                                              << 186                                G4int&    depth      ) 
180   if(depth > 100)                              << 187 {
181   {                                            << 188    if(depth > 100)
182     G4cout << "G4Integrator<T,F>::AdaptGauss:  << 189    {
183     G4cout << "Function varies too rapidly to  << 190      G4cout<<"G4Integrator<T,F>::AdaptGauss: WARNING !!!"<<G4endl  ;
184            << G4endl;                          << 191 G4cout
185                                                << 192 <<"Function varies too rapidly to get stated accuracy in 100 steps "<<G4endl ;
186     return;                                    << 193 
187   }                                            << 194      return ;
188   G4double xMean     = (xInitial + xFinal) / 2 << 195    }
189   G4double leftHalf  = Gauss(typeT, f, xInitia << 196    G4double xMean = (xInitial + xFinal)/2.0 ;
190   G4double rightHalf = Gauss(typeT, f, xMean,  << 197    G4double leftHalf  = Gauss(typeT,f,xInitial,xMean) ;
191   G4double full      = Gauss(typeT, f, xInitia << 198    G4double rightHalf = Gauss(typeT,f,xMean,xFinal) ;
192   if(std::fabs(leftHalf + rightHalf - full) <  << 199    G4double full = Gauss(typeT,f,xInitial,xFinal) ;
193   {                                            << 200    if(std::fabs(leftHalf+rightHalf-full) < fTolerance)
194     sum += full;                               << 201    {
195   }                                            << 202       sum += full ;
196   else                                         << 203    }
197   {                                            << 204    else
198     ++depth;                                   << 205    {
199     AdaptGauss(typeT, f, xInitial, xMean, fTol << 206       depth++ ;
200     AdaptGauss(typeT, f, xMean, xFinal, fToler << 207       AdaptGauss(typeT,f,xInitial,xMean,fTolerance,sum,depth) ;
201   }                                            << 208       AdaptGauss(typeT,f,xMean,xFinal,fTolerance,sum,depth) ;
202 }                                              << 209    }
203                                                << 210 }
204 template <class T, class F>                    << 211 
205 void G4Integrator<T, F>::AdaptGauss(T* ptrT, F << 212 template <class T, class F>  
206                                     G4double x << 213 void G4Integrator<T,F>::AdaptGauss( T* ptrT, F f, G4double  xInitial,
207                                     G4double&  << 214                                G4double  xFinal, G4double fTolerance,
                                                   >> 215                                G4double& sum,
                                                   >> 216                                G4int&    depth      ) 
208 {                                                 217 {
209   AdaptGauss(*ptrT, f, xInitial, xFinal, fTole << 218   AdaptGauss(*ptrT,f,xInitial,xFinal,fTolerance,sum,depth) ;
210 }                                                 219 }
211                                                   220 
212 //////////////////////////////////////////////    221 /////////////////////////////////////////////////////////////////////////
213 //                                                222 //
214 //                                                223 //
215 template <class T, class F>                       224 template <class T, class F>
216 void G4Integrator<T, F>::AdaptGauss(G4double ( << 225 void G4Integrator<T,F>::AdaptGauss( G4double (*f)(G4double), 
217                                     G4double x << 226                                G4double xInitial, G4double xFinal, 
218                                     G4double&  << 227                                G4double fTolerance, G4double& sum, 
219 {                                              << 228                                G4int& depth ) 
220   if(depth > 100)                              << 229 {
221   {                                            << 230    if(depth > 100)
222     G4cout << "G4SimpleIntegration::AdaptGauss << 231    {
223     G4cout << "Function varies too rapidly to  << 232      G4cout<<"G4SimpleIntegration::AdaptGauss: WARNING !!!"<<G4endl  ;
224            << G4endl;                          << 233      G4cout<<"Function varies too rapidly to get stated accuracy in 100 steps "
225                                                << 234            <<G4endl ;
226     return;                                    << 235 
227   }                                            << 236      return ;
228   G4double xMean     = (xInitial + xFinal) / 2 << 237    }
229   G4double leftHalf  = Gauss(f, xInitial, xMea << 238    G4double xMean = (xInitial + xFinal)/2.0 ;
230   G4double rightHalf = Gauss(f, xMean, xFinal) << 239    G4double leftHalf  = Gauss(f,xInitial,xMean) ;
231   G4double full      = Gauss(f, xInitial, xFin << 240    G4double rightHalf = Gauss(f,xMean,xFinal) ;
232   if(std::fabs(leftHalf + rightHalf - full) <  << 241    G4double full = Gauss(f,xInitial,xFinal) ;
233   {                                            << 242    if(std::fabs(leftHalf+rightHalf-full) < fTolerance)
234     sum += full;                               << 243    {
235   }                                            << 244       sum += full ;
236   else                                         << 245    }
237   {                                            << 246    else
238     ++depth;                                   << 247    {
239     AdaptGauss(f, xInitial, xMean, fTolerance, << 248       depth++ ;
240     AdaptGauss(f, xMean, xFinal, fTolerance, s << 249       AdaptGauss(f,xInitial,xMean,fTolerance,sum,depth) ;
241   }                                            << 250       AdaptGauss(f,xMean,xFinal,fTolerance,sum,depth) ;
                                                   >> 251    }
242 }                                                 252 }
243                                                   253 
                                                   >> 254 
                                                   >> 255 
                                                   >> 256 
244 //////////////////////////////////////////////    257 ////////////////////////////////////////////////////////////////////////
245 //                                                258 //
246 // Adaptive Gauss integration with accuracy 'e    259 // Adaptive Gauss integration with accuracy 'e'
247 // Convenient for using with class object type    260 // Convenient for using with class object typeT
248                                                << 261        
249 template <class T, class F>                    << 262 template<class T, class F> G4double 
250 G4double G4Integrator<T, F>::AdaptiveGauss(T&  << 263 G4Integrator<T,F>::AdaptiveGauss(  T& typeT, F f, G4double xInitial,
251                                            G4d << 264                                              G4double xFinal, G4double e   ) 
252 {                                              << 265 {
253   G4int depth  = 0;                            << 266    G4int depth = 0 ;
254   G4double sum = 0.0;                          << 267    G4double sum = 0.0 ;
255   AdaptGauss(typeT, f, xInitial, xFinal, e, su << 268    AdaptGauss(typeT,f,xInitial,xFinal,e,sum,depth) ;
256   return sum;                                  << 269    return sum ;
257 }                                                 270 }
258                                                   271 
259 //////////////////////////////////////////////    272 ////////////////////////////////////////////////////////////////////////
260 //                                                273 //
261 // Adaptive Gauss integration with accuracy 'e    274 // Adaptive Gauss integration with accuracy 'e'
262 // Convenient for using with 'this' pointer       275 // Convenient for using with 'this' pointer
263                                                << 276        
264 template <class T, class F>                    << 277 template<class T, class F> G4double 
265 G4double G4Integrator<T, F>::AdaptiveGauss(T*  << 278 G4Integrator<T,F>::AdaptiveGauss(  T* ptrT, F f, G4double xInitial,
266                                            G4d << 279                                              G4double xFinal, G4double e   ) 
267 {                                                 280 {
268   return AdaptiveGauss(*ptrT, f, xInitial, xFi << 281   return AdaptiveGauss(*ptrT,f,xInitial,xFinal,e) ;
269 }                                                 282 }
270                                                   283 
271 //////////////////////////////////////////////    284 ////////////////////////////////////////////////////////////////////////
272 //                                                285 //
273 // Adaptive Gauss integration with accuracy 'e    286 // Adaptive Gauss integration with accuracy 'e'
274 // Convenient for using with global scope func    287 // Convenient for using with global scope function f
275                                                << 288        
276 template <class T, class F>                    << 289 template <class T, class F> G4double 
277 G4double G4Integrator<T, F>::AdaptiveGauss(G4d << 290 G4Integrator<T,F>::AdaptiveGauss( G4double (*f)(G4double), 
278                                            G4d << 291                              G4double xInitial, G4double xFinal, G4double e ) 
279                                            G4d << 292 {
280 {                                              << 293    G4int depth = 0 ;
281   G4int depth  = 0;                            << 294    G4double sum = 0.0 ;
282   G4double sum = 0.0;                          << 295    AdaptGauss(f,xInitial,xFinal,e,sum,depth) ;
283   AdaptGauss(f, xInitial, xFinal, e, sum, dept << 296    return sum ;
284   return sum;                                  << 
285 }                                                 297 }
286                                                   298 
287 //////////////////////////////////////////////    299 ////////////////////////////////////////////////////////////////////////////
288 // Gauss integration methods involving ortogon    300 // Gauss integration methods involving ortogonal polynomials
289 //////////////////////////////////////////////    301 ////////////////////////////////////////////////////////////////////////////
290 //                                                302 //
291 // Methods involving Legendre polynomials      << 303 // Methods involving Legendre polynomials  
292 //                                                304 //
293 //////////////////////////////////////////////    305 /////////////////////////////////////////////////////////////////////////
294 //                                                306 //
295 // The value nLegendre set the accuracy requir    307 // The value nLegendre set the accuracy required, i.e the number of points
296 // where the function pFunction will be evalua    308 // where the function pFunction will be evaluated during integration.
297 // The function creates the arrays for absciss << 309 // The function creates the arrays for abscissas and weights that used 
298 // in Gauss-Legendre quadrature method.        << 310 // in Gauss-Legendre quadrature method. 
299 // The values a and b are the limits of integr    311 // The values a and b are the limits of integration of the function  f .
300 // nLegendre MUST BE EVEN !!!                     312 // nLegendre MUST BE EVEN !!!
301 // Returns the integral of the function f betw << 313 // Returns the integral of the function f between a and b, by 2*fNumber point 
302 // Gauss-Legendre integration: the function is    314 // Gauss-Legendre integration: the function is evaluated exactly
303 // 2*fNumber times at interior points in the r << 315 // 2*fNumber times at interior points in the range of integration. 
304 // Since the weights and abscissas are, in thi << 316 // Since the weights and abscissas are, in this case, symmetric around 
305 // the midpoint of the range of integration, t << 317 // the midpoint of the range of integration, there are actually only 
306 // fNumber distinct values of each.               318 // fNumber distinct values of each.
307 // Convenient for using with some class object    319 // Convenient for using with some class object dataT
308                                                   320 
309 template <class T, class F>                    << 321 template <class T, class F> G4double 
310 G4double G4Integrator<T, F>::Legendre(T& typeT << 322 G4Integrator<T,F>::Legendre( T& typeT, F f, G4double a, G4double b, G4int nLegendre)
311                                       G4int nL << 
312 {                                                 323 {
313   G4double nwt, nwt1, temp1, temp2, temp3, tem << 324    G4double newton, newton1, temp1, temp2, temp3, temp ;
314   G4double xDiff, xMean, dx, integral;         << 325    G4double xDiff, xMean, dx, integral ;
315                                                << 
316   const G4double tolerance = 1.6e-10;          << 
317   G4int i, j, k = nLegendre;                   << 
318   G4int fNumber = (nLegendre + 1) / 2;         << 
319                                                << 
320   if(2 * fNumber != k)                         << 
321   {                                            << 
322     G4Exception("G4Integrator<T,F>::Legendre(T << 
323                 FatalException, "Invalid (odd) << 
324   }                                            << 
325                                                << 
326   G4double* fAbscissa = new G4double[fNumber]; << 
327   G4double* fWeight   = new G4double[fNumber]; << 
328                                                << 
329   for(i = 1; i <= fNumber; ++i)  // Loop over  << 
330   {                                            << 
331     nwt = std::cos(CLHEP::pi * (i - 0.25) /    << 
332                    (k + 0.5));  // Initial roo << 
333                                                << 
334     do  // loop of Newton's method             << 
335     {                                          << 
336       temp1 = 1.0;                             << 
337       temp2 = 0.0;                             << 
338       for(j = 1; j <= k; ++j)                  << 
339       {                                        << 
340         temp3 = temp2;                         << 
341         temp2 = temp1;                         << 
342         temp1 = ((2.0 * j - 1.0) * nwt * temp2 << 
343       }                                        << 
344       temp = k * (nwt * temp1 - temp2) / (nwt  << 
345       nwt1 = nwt;                              << 
346       nwt  = nwt1 - temp1 / temp;  // Newton's << 
347     } while(std::fabs(nwt - nwt1) > tolerance) << 
348                                                << 
349     fAbscissa[fNumber - i] = nwt;              << 
350     fWeight[fNumber - i]   = 2.0 / ((1.0 - nwt << 
351   }                                            << 
352                                                   326 
353   //                                           << 327    const G4double tolerance = 1.6e-10 ;
354   // Now we ready to get integral              << 328    G4int i, j,   k = nLegendre ;
355   //                                           << 329    G4int fNumber = (nLegendre + 1)/2 ;
356                                                << 330 
357   xMean    = 0.5 * (a + b);                    << 331    if(2*fNumber != k)
358   xDiff    = 0.5 * (b - a);                    << 332    {
359   integral = 0.0;                              << 333       G4Exception("Invalid (odd) n Legendre in G4Integrator<T,F>::Legendre") ;
360   for(i = 0; i < fNumber; ++i)                 << 334    }
361   {                                            << 335 
362     dx = xDiff * fAbscissa[i];                 << 336    G4double* fAbscissa = new G4double[fNumber] ;
363     integral += fWeight[i] * ((typeT.*f)(xMean << 337    G4double* fWeight   = new G4double[fNumber] ;
364   }                                            << 338       
365   delete[] fAbscissa;                          << 339    for(i=1;i<=fNumber;i++)      // Loop over the desired roots
366   delete[] fWeight;                            << 340    {
367   return integral *= xDiff;                    << 341       newton = std::cos(pi*(i - 0.25)/(k + 0.5)) ;  // Initial root approximation
368 }                                              << 342 
                                                   >> 343       do     // loop of Newton's method  
                                                   >> 344       {                           
                                                   >> 345          temp1 = 1.0 ;
                                                   >> 346          temp2 = 0.0 ;
                                                   >> 347          for(j=1;j<=k;j++)
                                                   >> 348          {
                                                   >> 349             temp3 = temp2 ;
                                                   >> 350             temp2 = temp1 ;
                                                   >> 351             temp1 = ((2.0*j - 1.0)*newton*temp2 - (j - 1.0)*temp3)/j ;
                                                   >> 352          }
                                                   >> 353          temp = k*(newton*temp1 - temp2)/(newton*newton - 1.0) ;
                                                   >> 354          newton1 = newton ;
                                                   >> 355          newton  = newton1 - temp1/temp ;       // Newton's method
                                                   >> 356       }
                                                   >> 357       while(std::fabs(newton - newton1) > tolerance) ;
                                                   >> 358          
                                                   >> 359       fAbscissa[fNumber-i] =  newton ;
                                                   >> 360       fWeight[fNumber-i] = 2.0/((1.0 - newton*newton)*temp*temp) ;
                                                   >> 361    }
                                                   >> 362 //
                                                   >> 363 // Now we ready to get integral 
                                                   >> 364 //
                                                   >> 365    
                                                   >> 366    xMean = 0.5*(a + b) ;
                                                   >> 367    xDiff = 0.5*(b - a) ;
                                                   >> 368    integral = 0.0 ;
                                                   >> 369    for(i=0;i<fNumber;i++)
                                                   >> 370    {
                                                   >> 371       dx = xDiff*fAbscissa[i] ;
                                                   >> 372       integral += fWeight[i]*( (typeT.*f)(xMean + dx) + 
                                                   >> 373                                (typeT.*f)(xMean - dx)    ) ;
                                                   >> 374    }
                                                   >> 375    delete[] fAbscissa;
                                                   >> 376    delete[] fWeight;
                                                   >> 377    return integral *= xDiff ;
                                                   >> 378 } 
369                                                   379 
370 //////////////////////////////////////////////    380 ///////////////////////////////////////////////////////////////////////
371 //                                                381 //
372 // Convenient for using with the pointer 'this    382 // Convenient for using with the pointer 'this'
373                                                   383 
374 template <class T, class F>                    << 384 template <class T, class F> G4double 
375 G4double G4Integrator<T, F>::Legendre(T* ptrT, << 385 G4Integrator<T,F>::Legendre( T* ptrT, F f, G4double a, G4double b, G4int nLegendre) 
376                                       G4int nL << 
377 {                                                 386 {
378   return Legendre(*ptrT, f, a, b, nLegendre);  << 387   return Legendre(*ptrT,f,a,b,nLegendre) ;
379 }                                                 388 }
380                                                   389 
381 //////////////////////////////////////////////    390 ///////////////////////////////////////////////////////////////////////
382 //                                                391 //
383 // Convenient for using with global scope func    392 // Convenient for using with global scope function f
384                                                   393 
385 template <class T, class F>                       394 template <class T, class F>
386 G4double G4Integrator<T, F>::Legendre(G4double << 395 G4double G4Integrator<T,F>::
387                                       G4double << 396 Legendre( G4double (*f)(G4double), G4double a, G4double b, G4int nLegendre) 
388 {                                                 397 {
389   G4double nwt, nwt1, temp1, temp2, temp3, tem << 398    G4double newton, newton1, temp1, temp2, temp3, temp ;
390   G4double xDiff, xMean, dx, integral;         << 399    G4double xDiff, xMean, dx, integral ;
391                                                << 
392   const G4double tolerance = 1.6e-10;          << 
393   G4int i, j, k = nLegendre;                   << 
394   G4int fNumber = (nLegendre + 1) / 2;         << 
395                                                   400 
396   if(2 * fNumber != k)                         << 401    const G4double tolerance = 1.6e-10 ;
397   {                                            << 402    G4int i, j,   k = nLegendre ;
398     G4Exception("G4Integrator<T,F>::Legendre(. << 403    G4int fNumber = (nLegendre + 1)/2 ;
399                 FatalException, "Invalid (odd) << 404 
400   }                                            << 405    if(2*fNumber != k)
401                                                << 406    {
402   G4double* fAbscissa = new G4double[fNumber]; << 407       G4Exception("Invalid (odd) n Legendre in G4Integrator<T,F>::Legendre") ;
403   G4double* fWeight   = new G4double[fNumber]; << 408    }
404                                                << 409 
405   for(i = 1; i <= fNumber; i++)  // Loop over  << 410    G4double* fAbscissa = new G4double[fNumber] ;
406   {                                            << 411    G4double* fWeight   = new G4double[fNumber] ;
407     nwt = std::cos(CLHEP::pi * (i - 0.25) /    << 412       
408                    (k + 0.5));  // Initial roo << 413    for(i=1;i<=fNumber;i++)      // Loop over the desired roots
409                                                << 414    {
410     do  // loop of Newton's method             << 415       newton = std::cos(pi*(i - 0.25)/(k + 0.5)) ;  // Initial root approximation
411     {                                          << 416 
412       temp1 = 1.0;                             << 417       do     // loop of Newton's method  
413       temp2 = 0.0;                             << 418       {                           
414       for(j = 1; j <= k; ++j)                  << 419          temp1 = 1.0 ;
415       {                                        << 420          temp2 = 0.0 ;
416         temp3 = temp2;                         << 421          for(j=1;j<=k;j++)
417         temp2 = temp1;                         << 422          {
418         temp1 = ((2.0 * j - 1.0) * nwt * temp2 << 423             temp3 = temp2 ;
419       }                                        << 424             temp2 = temp1 ;
420       temp = k * (nwt * temp1 - temp2) / (nwt  << 425             temp1 = ((2.0*j - 1.0)*newton*temp2 - (j - 1.0)*temp3)/j ;
421       nwt1 = nwt;                              << 426          }
422       nwt  = nwt1 - temp1 / temp;  // Newton's << 427          temp = k*(newton*temp1 - temp2)/(newton*newton - 1.0) ;
423     } while(std::fabs(nwt - nwt1) > tolerance) << 428          newton1 = newton ;
424                                                << 429          newton  = newton1 - temp1/temp ;       // Newton's method
425     fAbscissa[fNumber - i] = nwt;              << 430       }
426     fWeight[fNumber - i]   = 2.0 / ((1.0 - nwt << 431       while(std::fabs(newton - newton1) > tolerance) ;
427   }                                            << 432          
                                                   >> 433       fAbscissa[fNumber-i] =  newton ;
                                                   >> 434       fWeight[fNumber-i] = 2.0/((1.0 - newton*newton)*temp*temp) ;
                                                   >> 435    }
                                                   >> 436 //
                                                   >> 437 // Now we ready to get integral 
                                                   >> 438 //
                                                   >> 439    
                                                   >> 440    xMean = 0.5*(a + b) ;
                                                   >> 441    xDiff = 0.5*(b - a) ;
                                                   >> 442    integral = 0.0 ;
                                                   >> 443    for(i=0;i<fNumber;i++)
                                                   >> 444    {
                                                   >> 445       dx = xDiff*fAbscissa[i] ;
                                                   >> 446       integral += fWeight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx)    ) ;
                                                   >> 447    }
                                                   >> 448    delete[] fAbscissa;
                                                   >> 449    delete[] fWeight;
428                                                   450 
429   //                                           << 451    return integral *= xDiff ;
430   // Now we ready to get integral              << 452 } 
431   //                                           << 
432                                                << 
433   xMean    = 0.5 * (a + b);                    << 
434   xDiff    = 0.5 * (b - a);                    << 
435   integral = 0.0;                              << 
436   for(i = 0; i < fNumber; ++i)                 << 
437   {                                            << 
438     dx = xDiff * fAbscissa[i];                 << 
439     integral += fWeight[i] * ((*f)(xMean + dx) << 
440   }                                            << 
441   delete[] fAbscissa;                          << 
442   delete[] fWeight;                            << 
443                                                << 
444   return integral *= xDiff;                    << 
445 }                                              << 
446                                                   453 
447 //////////////////////////////////////////////    454 ////////////////////////////////////////////////////////////////////////////
448 //                                                455 //
449 // Returns the integral of the function to be     456 // Returns the integral of the function to be pointed by T::f between a and b,
450 // by ten point Gauss-Legendre integration: th    457 // by ten point Gauss-Legendre integration: the function is evaluated exactly
451 // ten times at interior points in the range o    458 // ten times at interior points in the range of integration. Since the weights
452 // and abscissas are, in this case, symmetric  << 459 // and abscissas are, in this case, symmetric around the midpoint of the 
453 // range of integration, there are actually on    460 // range of integration, there are actually only five distinct values of each
454 // Convenient for using with class object type    461 // Convenient for using with class object typeT
455                                                   462 
456 template <class T, class F>                    << 463  template <class T, class F>  
457 G4double G4Integrator<T, F>::Legendre10(T& typ << 464  G4double G4Integrator<T,F>::Legendre10( T& typeT, F f,G4double a, G4double b) 
458 {                                                 465 {
459   G4int i;                                     << 466    G4int i ;
460   G4double xDiff, xMean, dx, integral;         << 467    G4double xDiff, xMean, dx, integral ;
461                                                << 468    
462   // From Abramowitz M., Stegan I.A. 1964 , Ha << 469    // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
463                                                << 470    
464   static const G4double abscissa[] = { 0.14887 << 471    static G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
465                                        0.67940 << 472                                   0.679409568299024, 0.865063366688985,
466                                        0.97390 << 473                                   0.973906528517172                      } ;
467                                                << 474    
468   static const G4double weight[] = { 0.2955242 << 475    static G4double weight[] =   { 0.295524224714753, 0.269266719309996, 
469                                      0.2190863 << 476                                   0.219086362515982, 0.149451349150581,
470                                      0.0666713 << 477                                   0.066671344308688                      } ;
471   xMean                          = 0.5 * (a +  << 478    xMean = 0.5*(a + b) ;
472   xDiff                          = 0.5 * (b -  << 479    xDiff = 0.5*(b - a) ;
473   integral                       = 0.0;        << 480    integral = 0.0 ;
474   for(i = 0; i < 5; ++i)                       << 481    for(i=0;i<5;i++)
475   {                                            << 482    {
476     dx = xDiff * abscissa[i];                  << 483      dx = xDiff*abscissa[i] ;
477     integral += weight[i] * ((typeT.*f)(xMean  << 484      integral += weight[i]*( (typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ;
478   }                                            << 485    }
479   return integral *= xDiff;                    << 486    return integral *= xDiff ;
480 }                                                 487 }
481                                                   488 
482 //////////////////////////////////////////////    489 ///////////////////////////////////////////////////////////////////////////
483 //                                                490 //
484 // Convenient for using with the pointer 'this    491 // Convenient for using with the pointer 'this'
485                                                   492 
486 template <class T, class F>                    << 493 template <class T, class F>  
487 G4double G4Integrator<T, F>::Legendre10(T* ptr << 494 G4double G4Integrator<T,F>::Legendre10( T* ptrT, F f,G4double a, G4double b)
488 {                                                 495 {
489   return Legendre10(*ptrT, f, a, b);           << 496   return Legendre10(*ptrT,f,a,b) ;
490 }                                              << 497 } 
491                                                   498 
492 //////////////////////////////////////////////    499 //////////////////////////////////////////////////////////////////////////
493 //                                                500 //
494 // Convenient for using with global scope func    501 // Convenient for using with global scope functions
495                                                   502 
496 template <class T, class F>                    << 503 template <class T, class F> G4double 
497 G4double G4Integrator<T, F>::Legendre10(G4doub << 504 G4Integrator<T,F>::Legendre10( G4double (*f)(G4double), G4double a, G4double b) 
498                                         G4doub << 
499 {                                                 505 {
500   G4int i;                                     << 506    G4int i ;
501   G4double xDiff, xMean, dx, integral;         << 507    G4double xDiff, xMean, dx, integral ;
502                                                << 508    
503   // From Abramowitz M., Stegan I.A. 1964 , Ha << 509    // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
504                                                << 510    
505   static const G4double abscissa[] = { 0.14887 << 511    static G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
506                                        0.67940 << 512                                   0.679409568299024, 0.865063366688985,
507                                        0.97390 << 513                                   0.973906528517172                      } ;
508                                                << 514    
509   static const G4double weight[] = { 0.2955242 << 515    static G4double weight[] =   { 0.295524224714753, 0.269266719309996, 
510                                      0.2190863 << 516                                   0.219086362515982, 0.149451349150581,
511                                      0.0666713 << 517                                   0.066671344308688                      } ;
512   xMean                          = 0.5 * (a +  << 518    xMean = 0.5*(a + b) ;
513   xDiff                          = 0.5 * (b -  << 519    xDiff = 0.5*(b - a) ;
514   integral                       = 0.0;        << 520    integral = 0.0 ;
515   for(i = 0; i < 5; ++i)                       << 521    for(i=0;i<5;i++)
516   {                                            << 522    {
517     dx = xDiff * abscissa[i];                  << 523      dx = xDiff*abscissa[i] ;
518     integral += weight[i] * ((*f)(xMean + dx)  << 524      integral += weight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx)) ;
519   }                                            << 525    }
520   return integral *= xDiff;                    << 526    return integral *= xDiff ;
521 }                                                 527 }
522                                                   528 
523 //////////////////////////////////////////////    529 ///////////////////////////////////////////////////////////////////////
524 //                                                530 //
525 // Returns the integral of the function to be     531 // Returns the integral of the function to be pointed by T::f between a and b,
526 // by 96 point Gauss-Legendre integration: the    532 // by 96 point Gauss-Legendre integration: the function is evaluated exactly
527 // ten Times at interior points in the range o    533 // ten Times at interior points in the range of integration. Since the weights
528 // and abscissas are, in this case, symmetric  << 534 // and abscissas are, in this case, symmetric around the midpoint of the 
529 // range of integration, there are actually on    535 // range of integration, there are actually only five distinct values of each
530 // Convenient for using with some class object    536 // Convenient for using with some class object typeT
531                                                   537 
532 template <class T, class F>                    << 538 template <class T, class F>  
533 G4double G4Integrator<T, F>::Legendre96(T& typ << 539 G4double G4Integrator<T,F>::Legendre96( T& typeT, F f,G4double a, G4double b) 
534 {                                                 540 {
535   G4int i;                                     << 541    G4int i ;
536   G4double xDiff, xMean, dx, integral;         << 542    G4double xDiff, xMean, dx, integral ;
537                                                << 543    
538   // From Abramowitz M., Stegan I.A. 1964 , Ha << 544    // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
539                                                << 545    
540   static const G4double abscissa[] = {         << 546    static G4double 
541     0.016276744849602969579, 0.048812985136049 << 547    abscissa[] = { 
542     0.081297495464425558994, 0.113695850110665 << 548                   0.016276744849602969579, 0.048812985136049731112,
543     0.145973714654896941989, 0.178096882367618 << 549                   0.081297495464425558994, 0.113695850110665920911,
544                                                << 550                   0.145973714654896941989, 0.178096882367618602759,  // 6
545     0.210031310460567203603, 0.241743156163840 << 551                            
546     0.273198812591049141487, 0.304364944354496 << 552                   0.210031310460567203603, 0.241743156163840012328,
547     0.335208522892625422616, 0.365696861472313 << 553                   0.273198812591049141487, 0.304364944354496353024,
548                                                << 554                   0.335208522892625422616, 0.365696861472313635031,  // 12
549     0.395797649828908603285, 0.425478988407300 << 555                            
550     0.454709422167743008636, 0.483457973920596 << 556                   0.395797649828908603285, 0.425478988407300545365,
551     0.511694177154667673586, 0.539388108324357 << 557                   0.454709422167743008636, 0.483457973920596359768,
552                                                << 558                   0.511694177154667673586, 0.539388108324357436227,  // 18
553     0.566510418561397168404, 0.593032364777572 << 559                            
554     0.618925840125468570386, 0.644163403784967 << 560                   0.566510418561397168404, 0.593032364777572080684,
555     0.668718310043916153953, 0.692564536642171 << 561                   0.618925840125468570386, 0.644163403784967106798,
556                                                << 562                   0.668718310043916153953, 0.692564536642171561344,  // 24
557     0.715676812348967626225, 0.738030643744400 << 563                            
558     0.759602341176647498703, 0.780369043867433 << 564                   0.715676812348967626225, 0.738030643744400132851,
559     0.800308744139140817229, 0.819400310737931 << 565                   0.759602341176647498703, 0.780369043867433217604,
560                                                << 566                   0.800308744139140817229, 0.819400310737931675539,  // 30
561     0.837623511228187121494, 0.854959033434601 << 567                            
562     0.871388505909296502874, 0.886894517402420 << 568                   0.837623511228187121494, 0.854959033434601455463,
563     0.901460635315852341319, 0.915071423120898 << 569                   0.871388505909296502874, 0.886894517402420416057,
564                                                << 570                   0.901460635315852341319, 0.915071423120898074206,  // 36
565     0.927712456722308690965, 0.939370339752755 << 571                            
566     0.950032717784437635756, 0.959688291448742 << 572                   0.927712456722308690965, 0.939370339752755216932,
567     0.968326828463264212174, 0.975939174585136 << 573                   0.950032717784437635756, 0.959688291448742539300,
568                                                << 574                   0.968326828463264212174, 0.975939174585136466453,  // 42
569     0.982517263563014677447, 0.988054126329623 << 575                            
570     0.992543900323762624572, 0.995981842987209 << 576                   0.982517263563014677447, 0.988054126329623799481,
571     0.998364375863181677724, 0.999689503883230 << 577                   0.992543900323762624572, 0.995981842987209290650,
572   };                                           << 578                   0.998364375863181677724, 0.999689503883230766828   // 48
573                                                << 579                                                                             } ;
574   static const G4double weight[] = {           << 580    
575     0.032550614492363166242, 0.032516118713868 << 581    static G4double 
576     0.032447163714064269364, 0.032343822568575 << 582    weight[] = {  
577     0.032206204794030250669, 0.032034456231992 << 583                   0.032550614492363166242, 0.032516118713868835987,
578                                                << 584                   0.032447163714064269364, 0.032343822568575928429,
579     0.031828758894411006535, 0.031589330770727 << 585                   0.032206204794030250669, 0.032034456231992663218,  // 6
580     0.031316425596862355813, 0.031010332586313 << 586                            
581     0.030671376123669149014, 0.030299915420827 << 587                   0.031828758894411006535, 0.031589330770727168558,
582                                                << 588                   0.031316425596862355813, 0.031010332586313837423,
583     0.029896344136328385984, 0.029461089958167 << 589                   0.030671376123669149014, 0.030299915420827593794,  // 12
584     0.028994614150555236543, 0.028497411065085 << 590                            
585     0.027970007616848334440, 0.027412962726029 << 591                   0.029896344136328385984, 0.029461089958167905970,
586                                                << 592                   0.028994614150555236543, 0.028497411065085385646,
587     0.026826866725591762198, 0.026212340735672 << 593                   0.027970007616848334440, 0.027412962726029242823,  // 18
588     0.025570036005349361499, 0.024900633222483 << 594                            
589     0.024204841792364691282, 0.023483399085926 << 595                   0.026826866725591762198, 0.026212340735672413913,
590                                                << 596                   0.025570036005349361499, 0.024900633222483610288,
591     0.022737069658329374001, 0.021966644438744 << 597                   0.024204841792364691282, 0.023483399085926219842,  // 24
592     0.021172939892191298988, 0.020356797154333 << 598                            
593     0.019519081140145022410, 0.018660679627411 << 599                   0.022737069658329374001, 0.021966644438744349195,
594                                                << 600                   0.021172939892191298988, 0.020356797154333324595,
595     0.017782502316045260838, 0.016885479864245 << 601                   0.019519081140145022410, 0.018660679627411467385,  // 30
596     0.015970562902562291381, 0.015038721026994 << 602                            
597     0.014090941772314860916, 0.013128229566961 << 603                   0.017782502316045260838, 0.016885479864245172450,
598                                                << 604                   0.015970562902562291381, 0.015038721026994938006,
599     0.012151604671088319635, 0.011162102099838 << 605                   0.014090941772314860916, 0.013128229566961572637,  // 36
600     0.010160770535008415758, 0.009148671230783 << 606                            
601     0.008126876925698759217, 0.007096470791153 << 607                   0.012151604671088319635, 0.011162102099838498591,
602                                                << 608                   0.010160770535008415758, 0.009148671230783386633,
603     0.006058545504235961683, 0.005014202742927 << 609                   0.008126876925698759217, 0.007096470791153865269,  // 42
604     0.003964554338444686674, 0.002910731817934 << 610                            
605     0.001853960788946921732, 0.000796792065552 << 611                   0.006058545504235961683, 0.005014202742927517693,
606   };                                           << 612                   0.003964554338444686674, 0.002910731817934946408,
607   xMean    = 0.5 * (a + b);                    << 613                   0.001853960788946921732, 0.000796792065552012429   // 48
608   xDiff    = 0.5 * (b - a);                    << 614                                                                             } ;
609   integral = 0.0;                              << 615    xMean = 0.5*(a + b) ;
610   for(i = 0; i < 48; ++i)                      << 616    xDiff = 0.5*(b - a) ;
611   {                                            << 617    integral = 0.0 ;
612     dx = xDiff * abscissa[i];                  << 618    for(i=0;i<48;i++)
613     integral += weight[i] * ((typeT.*f)(xMean  << 619    {
614   }                                            << 620       dx = xDiff*abscissa[i] ;
615   return integral *= xDiff;                    << 621       integral += weight[i]*((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ;
                                                   >> 622    }
                                                   >> 623    return integral *= xDiff ;
616 }                                                 624 }
617                                                   625 
618 //////////////////////////////////////////////    626 ///////////////////////////////////////////////////////////////////////
619 //                                                627 //
620 // Convenient for using with the pointer 'this    628 // Convenient for using with the pointer 'this'
621                                                   629 
622 template <class T, class F>                    << 630 template <class T, class F>  
623 G4double G4Integrator<T, F>::Legendre96(T* ptr << 631 G4double G4Integrator<T,F>::Legendre96( T* ptrT, F f,G4double a, G4double b)
624 {                                                 632 {
625   return Legendre96(*ptrT, f, a, b);           << 633   return Legendre96(*ptrT,f,a,b) ;
626 }                                              << 634 } 
627                                                   635 
628 //////////////////////////////////////////////    636 ///////////////////////////////////////////////////////////////////////
629 //                                                637 //
630 // Convenient for using with global scope func << 638 // Convenient for using with global scope function f 
631                                                   639 
632 template <class T, class F>                    << 640 template <class T, class F> G4double 
633 G4double G4Integrator<T, F>::Legendre96(G4doub << 641 G4Integrator<T,F>::Legendre96( G4double (*f)(G4double), G4double a, G4double b) 
634                                         G4doub << 
635 {                                                 642 {
636   G4int i;                                     << 643    G4int i ;
637   G4double xDiff, xMean, dx, integral;         << 644    G4double xDiff, xMean, dx, integral ;
638                                                << 645    
639   // From Abramowitz M., Stegan I.A. 1964 , Ha << 646    // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
640                                                << 647    
641   static const G4double abscissa[] = {         << 648    static G4double 
642     0.016276744849602969579, 0.048812985136049 << 649    abscissa[] = { 
643     0.081297495464425558994, 0.113695850110665 << 650                   0.016276744849602969579, 0.048812985136049731112,
644     0.145973714654896941989, 0.178096882367618 << 651                   0.081297495464425558994, 0.113695850110665920911,
645                                                << 652                   0.145973714654896941989, 0.178096882367618602759,  // 6
646     0.210031310460567203603, 0.241743156163840 << 653                            
647     0.273198812591049141487, 0.304364944354496 << 654                   0.210031310460567203603, 0.241743156163840012328,
648     0.335208522892625422616, 0.365696861472313 << 655                   0.273198812591049141487, 0.304364944354496353024,
649                                                << 656                   0.335208522892625422616, 0.365696861472313635031,  // 12
650     0.395797649828908603285, 0.425478988407300 << 657                            
651     0.454709422167743008636, 0.483457973920596 << 658                   0.395797649828908603285, 0.425478988407300545365,
652     0.511694177154667673586, 0.539388108324357 << 659                   0.454709422167743008636, 0.483457973920596359768,
653                                                << 660                   0.511694177154667673586, 0.539388108324357436227,  // 18
654     0.566510418561397168404, 0.593032364777572 << 661                            
655     0.618925840125468570386, 0.644163403784967 << 662                   0.566510418561397168404, 0.593032364777572080684,
656     0.668718310043916153953, 0.692564536642171 << 663                   0.618925840125468570386, 0.644163403784967106798,
657                                                << 664                   0.668718310043916153953, 0.692564536642171561344,  // 24
658     0.715676812348967626225, 0.738030643744400 << 665                            
659     0.759602341176647498703, 0.780369043867433 << 666                   0.715676812348967626225, 0.738030643744400132851,
660     0.800308744139140817229, 0.819400310737931 << 667                   0.759602341176647498703, 0.780369043867433217604,
661                                                << 668                   0.800308744139140817229, 0.819400310737931675539,  // 30
662     0.837623511228187121494, 0.854959033434601 << 669                            
663     0.871388505909296502874, 0.886894517402420 << 670                   0.837623511228187121494, 0.854959033434601455463,
664     0.901460635315852341319, 0.915071423120898 << 671                   0.871388505909296502874, 0.886894517402420416057,
665                                                << 672                   0.901460635315852341319, 0.915071423120898074206,  // 36
666     0.927712456722308690965, 0.939370339752755 << 673                            
667     0.950032717784437635756, 0.959688291448742 << 674                   0.927712456722308690965, 0.939370339752755216932,
668     0.968326828463264212174, 0.975939174585136 << 675                   0.950032717784437635756, 0.959688291448742539300,
669                                                << 676                   0.968326828463264212174, 0.975939174585136466453,  // 42
670     0.982517263563014677447, 0.988054126329623 << 677                            
671     0.992543900323762624572, 0.995981842987209 << 678                   0.982517263563014677447, 0.988054126329623799481,
672     0.998364375863181677724, 0.999689503883230 << 679                   0.992543900323762624572, 0.995981842987209290650,
673   };                                           << 680                   0.998364375863181677724, 0.999689503883230766828   // 48
674                                                << 681                                                                             } ;
675   static const G4double weight[] = {           << 682    
676     0.032550614492363166242, 0.032516118713868 << 683    static G4double 
677     0.032447163714064269364, 0.032343822568575 << 684    weight[] = {  
678     0.032206204794030250669, 0.032034456231992 << 685                   0.032550614492363166242, 0.032516118713868835987,
679                                                << 686                   0.032447163714064269364, 0.032343822568575928429,
680     0.031828758894411006535, 0.031589330770727 << 687                   0.032206204794030250669, 0.032034456231992663218,  // 6
681     0.031316425596862355813, 0.031010332586313 << 688                            
682     0.030671376123669149014, 0.030299915420827 << 689                   0.031828758894411006535, 0.031589330770727168558,
683                                                << 690                   0.031316425596862355813, 0.031010332586313837423,
684     0.029896344136328385984, 0.029461089958167 << 691                   0.030671376123669149014, 0.030299915420827593794,  // 12
685     0.028994614150555236543, 0.028497411065085 << 692                            
686     0.027970007616848334440, 0.027412962726029 << 693                   0.029896344136328385984, 0.029461089958167905970,
687                                                << 694                   0.028994614150555236543, 0.028497411065085385646,
688     0.026826866725591762198, 0.026212340735672 << 695                   0.027970007616848334440, 0.027412962726029242823,  // 18
689     0.025570036005349361499, 0.024900633222483 << 696                            
690     0.024204841792364691282, 0.023483399085926 << 697                   0.026826866725591762198, 0.026212340735672413913,
691                                                << 698                   0.025570036005349361499, 0.024900633222483610288,
692     0.022737069658329374001, 0.021966644438744 << 699                   0.024204841792364691282, 0.023483399085926219842,  // 24
693     0.021172939892191298988, 0.020356797154333 << 700                            
694     0.019519081140145022410, 0.018660679627411 << 701                   0.022737069658329374001, 0.021966644438744349195,
695                                                << 702                   0.021172939892191298988, 0.020356797154333324595,
696     0.017782502316045260838, 0.016885479864245 << 703                   0.019519081140145022410, 0.018660679627411467385,  // 30
697     0.015970562902562291381, 0.015038721026994 << 704                            
698     0.014090941772314860916, 0.013128229566961 << 705                   0.017782502316045260838, 0.016885479864245172450,
699                                                << 706                   0.015970562902562291381, 0.015038721026994938006,
700     0.012151604671088319635, 0.011162102099838 << 707                   0.014090941772314860916, 0.013128229566961572637,  // 36
701     0.010160770535008415758, 0.009148671230783 << 708                            
702     0.008126876925698759217, 0.007096470791153 << 709                   0.012151604671088319635, 0.011162102099838498591,
703                                                << 710                   0.010160770535008415758, 0.009148671230783386633,
704     0.006058545504235961683, 0.005014202742927 << 711                   0.008126876925698759217, 0.007096470791153865269,  // 42
705     0.003964554338444686674, 0.002910731817934 << 712                            
706     0.001853960788946921732, 0.000796792065552 << 713                   0.006058545504235961683, 0.005014202742927517693,
707   };                                           << 714                   0.003964554338444686674, 0.002910731817934946408,
708   xMean    = 0.5 * (a + b);                    << 715                   0.001853960788946921732, 0.000796792065552012429   // 48
709   xDiff    = 0.5 * (b - a);                    << 716                                                                             } ;
710   integral = 0.0;                              << 717    xMean = 0.5*(a + b) ;
711   for(i = 0; i < 48; ++i)                      << 718    xDiff = 0.5*(b - a) ;
712   {                                            << 719    integral = 0.0 ;
713     dx = xDiff * abscissa[i];                  << 720    for(i=0;i<48;i++)
714     integral += weight[i] * ((*f)(xMean + dx)  << 721    {
715   }                                            << 722       dx = xDiff*abscissa[i] ;
716   return integral *= xDiff;                    << 723       integral += weight[i]*((*f)(xMean + dx) + (*f)(xMean - dx)) ;
                                                   >> 724    }
                                                   >> 725    return integral *= xDiff ;
717 }                                                 726 }
718                                                   727 
719 //////////////////////////////////////////////    728 //////////////////////////////////////////////////////////////////////////////
720 //                                                729 //
721 // Methods involving Chebyshev polynomials     << 730 // Methods involving Chebyshev polynomials 
722 //                                                731 //
723 //////////////////////////////////////////////    732 ///////////////////////////////////////////////////////////////////////////
724 //                                                733 //
725 // Integrates function pointed by T::f from a  << 734 // Integrates function pointed by T::f from a to b by Gauss-Chebyshev 
726 // quadrature method.                             735 // quadrature method.
727 // Convenient for using with class object type    736 // Convenient for using with class object typeT
728                                                   737 
729 template <class T, class F>                    << 738 template <class T, class F> G4double 
730 G4double G4Integrator<T, F>::Chebyshev(T& type << 739 G4Integrator<T,F>::Chebyshev( T& typeT, F f, G4double a, 
731                                        G4int n << 740                          G4double b, G4int nChebyshev ) 
732 {                                              << 741 {
733   G4int i;                                     << 742    G4int i ;
734   G4double xDiff, xMean, dx, integral = 0.0;   << 743    G4double xDiff, xMean, dx, integral = 0.0 ;
735                                                << 744    
736   G4int fNumber       = nChebyshev;  // Try to << 745    G4int fNumber = nChebyshev  ;   // Try to reduce fNumber twice ??
737   G4double cof        = CLHEP::pi / fNumber;   << 746    G4double cof = pi/fNumber ;
738   G4double* fAbscissa = new G4double[fNumber]; << 747    G4double* fAbscissa = new G4double[fNumber] ;
739   G4double* fWeight   = new G4double[fNumber]; << 748    G4double* fWeight   = new G4double[fNumber] ;
740   for(i = 0; i < fNumber; ++i)                 << 749    for(i=0;i<fNumber;i++)
741   {                                            << 750    {
742     fAbscissa[i] = std::cos(cof * (i + 0.5));  << 751       fAbscissa[i] = std::cos(cof*(i + 0.5)) ;
743     fWeight[i]   = cof * std::sqrt(1 - fAbscis << 752       fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ;
744   }                                            << 753    }
745                                                << 754 //
746   //                                           << 755 // Now we ready to estimate the integral
747   // Now we ready to estimate the integral     << 756 //
748   //                                           << 757    xMean = 0.5*(a + b) ;
749                                                << 758    xDiff = 0.5*(b - a) ;
750   xMean = 0.5 * (a + b);                       << 759    for(i=0;i<fNumber;i++)
751   xDiff = 0.5 * (b - a);                       << 760    {
752   for(i = 0; i < fNumber; ++i)                 << 761       dx = xDiff*fAbscissa[i] ;
753   {                                            << 762       integral += fWeight[i]*(typeT.*f)(xMean + dx)  ;
754     dx = xDiff * fAbscissa[i];                 << 763    }
755     integral += fWeight[i] * (typeT.*f)(xMean  << 764    delete[] fAbscissa;
756   }                                            << 765    delete[] fWeight;
757   delete[] fAbscissa;                          << 766    return integral *= xDiff ;
758   delete[] fWeight;                            << 
759   return integral *= xDiff;                    << 
760 }                                                 767 }
761                                                   768 
762 //////////////////////////////////////////////    769 ///////////////////////////////////////////////////////////////////////
763 //                                                770 //
764 // Convenient for using with 'this' pointer       771 // Convenient for using with 'this' pointer
765                                                   772 
766 template <class T, class F>                    << 773 template <class T, class F> G4double 
767 G4double G4Integrator<T, F>::Chebyshev(T* ptrT << 774 G4Integrator<T,F>::Chebyshev( T* ptrT, F f, G4double a, G4double b, G4int n)
768                                        G4int n << 
769 {                                                 775 {
770   return Chebyshev(*ptrT, f, a, b, n);         << 776   return Chebyshev(*ptrT,f,a,b,n) ;
771 }                                              << 777 } 
772                                                   778 
773 //////////////////////////////////////////////    779 ////////////////////////////////////////////////////////////////////////
774 //                                                780 //
775 // For use with global scope functions f       << 781 // For use with global scope functions f 
776                                                   782 
777 template <class T, class F>                    << 783 template <class T, class F> G4double 
778 G4double G4Integrator<T, F>::Chebyshev(G4doubl << 784 G4Integrator<T,F>::Chebyshev( G4double (*f)(G4double), 
779                                        G4doubl << 785                          G4double a, G4double b, G4int nChebyshev) 
780 {                                              << 786 {
781   G4int i;                                     << 787    G4int i ;
782   G4double xDiff, xMean, dx, integral = 0.0;   << 788    G4double xDiff, xMean, dx, integral = 0.0 ;
783                                                << 789    
784   G4int fNumber       = nChebyshev;  // Try to << 790    G4int fNumber = nChebyshev  ;   // Try to reduce fNumber twice ??
785   G4double cof        = CLHEP::pi / fNumber;   << 791    G4double cof = pi/fNumber ;
786   G4double* fAbscissa = new G4double[fNumber]; << 792    G4double* fAbscissa = new G4double[fNumber] ;
787   G4double* fWeight   = new G4double[fNumber]; << 793    G4double* fWeight   = new G4double[fNumber] ;
788   for(i = 0; i < fNumber; ++i)                 << 794    for(i=0;i<fNumber;i++)
789   {                                            << 795    {
790     fAbscissa[i] = std::cos(cof * (i + 0.5));  << 796       fAbscissa[i] = std::cos(cof*(i + 0.5)) ;
791     fWeight[i]   = cof * std::sqrt(1 - fAbscis << 797       fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ;
792   }                                            << 798    }
793                                                << 799 //
794   //                                           << 800 // Now we ready to estimate the integral
795   // Now we ready to estimate the integral     << 801 //
796   //                                           << 802    xMean = 0.5*(a + b) ;
797                                                << 803    xDiff = 0.5*(b - a) ;
798   xMean = 0.5 * (a + b);                       << 804    for(i=0;i<fNumber;i++)
799   xDiff = 0.5 * (b - a);                       << 805    {
800   for(i = 0; i < fNumber; ++i)                 << 806       dx = xDiff*fAbscissa[i] ;
801   {                                            << 807       integral += fWeight[i]*(*f)(xMean + dx)  ;
802     dx = xDiff * fAbscissa[i];                 << 808    }
803     integral += fWeight[i] * (*f)(xMean + dx); << 809    delete[] fAbscissa;
804   }                                            << 810    delete[] fWeight;
805   delete[] fAbscissa;                          << 811    return integral *= xDiff ;
806   delete[] fWeight;                            << 
807   return integral *= xDiff;                    << 
808 }                                                 812 }
809                                                   813 
810 //////////////////////////////////////////////    814 //////////////////////////////////////////////////////////////////////
811 //                                                815 //
812 // Method involving Laguerre polynomials          816 // Method involving Laguerre polynomials
813 //                                                817 //
814 //////////////////////////////////////////////    818 //////////////////////////////////////////////////////////////////////
815 //                                                819 //
816 // Integral from zero to infinity of std::pow( << 820 // Integral from zero to infinity of std::pow(x,alpha)*std::exp(-x)*f(x). 
817 // The value of nLaguerre sets the accuracy.      821 // The value of nLaguerre sets the accuracy.
818 // The function creates arrays fAbscissa[0,.., << 822 // The function creates arrays fAbscissa[0,..,nLaguerre-1] and 
819 // fWeight[0,..,nLaguerre-1] .                 << 823 // fWeight[0,..,nLaguerre-1] . 
820 // Convenient for using with class object 'typ    824 // Convenient for using with class object 'typeT' and (typeT.*f) function
821 // (T::f)                                         825 // (T::f)
822                                                   826 
823 template <class T, class F>                    << 827 template <class T, class F> G4double 
824 G4double G4Integrator<T, F>::Laguerre(T& typeT << 828 G4Integrator<T,F>::Laguerre( T& typeT, F f, G4double alpha, G4int nLaguerre ) 
825                                       G4int nL << 
826 {                                                 829 {
827   const G4double tolerance = 1.0e-10;          << 830    const G4double tolerance = 1.0e-10 ;
828   const G4int maxNumber    = 12;               << 831    const G4int maxNumber = 12 ;
829   G4int i, j, k;                               << 832    G4int i, j, k ;
830   G4double nwt      = 0., nwt1, temp1, temp2,  << 833    G4double newton=0., newton1, temp1, temp2, temp3, temp, cofi ;
831   G4double integral = 0.0;                     << 834    G4double integral = 0.0 ;
832                                                << 835 
833   G4int fNumber       = nLaguerre;             << 836    G4int fNumber = nLaguerre ;
834   G4double* fAbscissa = new G4double[fNumber]; << 837    G4double* fAbscissa = new G4double[fNumber] ;
835   G4double* fWeight   = new G4double[fNumber]; << 838    G4double* fWeight   = new G4double[fNumber] ;
                                                   >> 839       
                                                   >> 840    for(i=1;i<=fNumber;i++)      // Loop over the desired roots
                                                   >> 841    {
                                                   >> 842       if(i == 1)
                                                   >> 843       {
                                                   >> 844          newton = (1.0 + alpha)*(3.0 + 0.92*alpha)/(1.0 + 2.4*fNumber + 1.8*alpha) ;
                                                   >> 845       }
                                                   >> 846       else if(i == 2)
                                                   >> 847       {
                                                   >> 848          newton += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
                                                   >> 849       }
                                                   >> 850       else
                                                   >> 851       {
                                                   >> 852          cofi = i - 2 ;
                                                   >> 853          newton += ((1.0+2.55*cofi)/(1.9*cofi) + 1.26*cofi*alpha/(1.0+3.5*cofi))*
                                                   >> 854                    (newton - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
                                                   >> 855       }
                                                   >> 856       for(k=1;k<=maxNumber;k++)
                                                   >> 857       {
                                                   >> 858          temp1 = 1.0 ;
                                                   >> 859          temp2 = 0.0 ;
836                                                   860 
837   for(i = 1; i <= fNumber; ++i)  // Loop over  << 861          for(j=1;j<=fNumber;j++)
838   {                                            << 862          {
839     if(i == 1)                                 << 863             temp3 = temp2 ;
840     {                                          << 864             temp2 = temp1 ;
841       nwt = (1.0 + alpha) * (3.0 + 0.92 * alph << 865          temp1 = ((2*j - 1 + alpha - newton)*temp2 - (j - 1 + alpha)*temp3)/j ;
842             (1.0 + 2.4 * fNumber + 1.8 * alpha << 866          }
843     }                                          << 867          temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/newton ;
844     else if(i == 2)                            << 868          newton1 = newton ;
845     {                                          << 869          newton  = newton1 - temp1/temp ;
846       nwt += (15.0 + 6.25 * alpha) / (1.0 + 0. << 870 
847     }                                          << 871          if(std::fabs(newton - newton1) <= tolerance) 
848     else                                       << 872          {
849     {                                          << 873             break ;
850       cofi = i - 2;                            << 874          }
851       nwt += ((1.0 + 2.55 * cofi) / (1.9 * cof << 875       }
852               1.26 * cofi * alpha / (1.0 + 3.5 << 876       if(k > maxNumber)
853              (nwt - fAbscissa[i - 3]) / (1.0 + << 877       {
854     }                                          << 878          G4Exception("Too many (>12) iterations in G4Integration::Laguerre") ;
855     for(k = 1; k <= maxNumber; ++k)            << 879       }
856     {                                          << 880          
857       temp1 = 1.0;                             << 881       fAbscissa[i-1] =  newton ;
858       temp2 = 0.0;                             << 882       fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) - 
859                                                << 883                 GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
860       for(j = 1; j <= fNumber; ++j)            << 884    }
861       {                                        << 885 //
862         temp3 = temp2;                         << 886 // Integral evaluation
863         temp2 = temp1;                         << 887 //
864         temp1 =                                << 888    for(i=0;i<fNumber;i++)
865           ((2 * j - 1 + alpha - nwt) * temp2 - << 889    {
866       }                                        << 890       integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ;
867       temp = (fNumber * temp1 - (fNumber + alp << 891    }
868       nwt1 = nwt;                              << 892    delete[] fAbscissa;
869       nwt  = nwt1 - temp1 / temp;              << 893    delete[] fWeight;
870                                                << 894    return integral ;
871       if(std::fabs(nwt - nwt1) <= tolerance)   << 895 }
872       {                                        << 
873         break;                                 << 
874       }                                        << 
875     }                                          << 
876     if(k > maxNumber)                          << 
877     {                                          << 
878       G4Exception("G4Integrator<T,F>::Laguerre << 
879                   FatalException, "Too many (> << 
880     }                                          << 
881                                                << 
882     fAbscissa[i - 1] = nwt;                    << 
883     fWeight[i - 1]   = -std::exp(GammaLogarith << 
884                                GammaLogarithm( << 
885                      (temp * fNumber * temp2); << 
886   }                                            << 
887                                                   896 
888   //                                           << 
889   // Integral evaluation                       << 
890   //                                           << 
891                                                   897 
892   for(i = 0; i < fNumber; ++i)                 << 
893   {                                            << 
894     integral += fWeight[i] * (typeT.*f)(fAbsci << 
895   }                                            << 
896   delete[] fAbscissa;                          << 
897   delete[] fWeight;                            << 
898   return integral;                             << 
899 }                                              << 
900                                                   898 
901 //////////////////////////////////////////////    899 //////////////////////////////////////////////////////////////////////
902 //                                                900 //
903 //                                                901 //
904                                                   902 
905 template <class T, class F>                    << 903 template <class T, class F> G4double 
906 G4double G4Integrator<T, F>::Laguerre(T* ptrT, << 904 G4Integrator<T,F>::Laguerre( T* ptrT, F f, G4double alpha, G4int nLaguerre ) 
907                                       G4int nL << 
908 {                                                 905 {
909   return Laguerre(*ptrT, f, alpha, nLaguerre); << 906   return Laguerre(*ptrT,f,alpha,nLaguerre) ;
910 }                                                 907 }
911                                                   908 
912 //////////////////////////////////////////////    909 ////////////////////////////////////////////////////////////////////////
913 //                                                910 //
914 // For use with global scope functions f       << 911 // For use with global scope functions f 
915                                                << 
916 template <class T, class F>                    << 
917 G4double G4Integrator<T, F>::Laguerre(G4double << 
918                                       G4int nL << 
919 {                                              << 
920   const G4double tolerance = 1.0e-10;          << 
921   const G4int maxNumber    = 12;               << 
922   G4int i, j, k;                               << 
923   G4double nwt      = 0., nwt1, temp1, temp2,  << 
924   G4double integral = 0.0;                     << 
925                                                << 
926   G4int fNumber       = nLaguerre;             << 
927   G4double* fAbscissa = new G4double[fNumber]; << 
928   G4double* fWeight   = new G4double[fNumber]; << 
929                                                   912 
930   for(i = 1; i <= fNumber; ++i)  // Loop over  << 913 template <class T, class F> G4double 
931   {                                            << 914 G4Integrator<T,F>::Laguerre( G4double (*f)(G4double), 
932     if(i == 1)                                 << 915                              G4double alpha, G4int nLaguerre) 
933     {                                          << 916 {
934       nwt = (1.0 + alpha) * (3.0 + 0.92 * alph << 917    const G4double tolerance = 1.0e-10 ;
935             (1.0 + 2.4 * fNumber + 1.8 * alpha << 918    const G4int maxNumber = 12 ;
936     }                                          << 919    G4int i, j, k ;
937     else if(i == 2)                            << 920    G4double newton=0., newton1, temp1, temp2, temp3, temp, cofi ;
938     {                                          << 921    G4double integral = 0.0 ;
939       nwt += (15.0 + 6.25 * alpha) / (1.0 + 0. << 922 
940     }                                          << 923    G4int fNumber = nLaguerre ;
941     else                                       << 924    G4double* fAbscissa = new G4double[fNumber] ;
942     {                                          << 925    G4double* fWeight   = new G4double[fNumber] ;
943       cofi = i - 2;                            << 926       
944       nwt += ((1.0 + 2.55 * cofi) / (1.9 * cof << 927    for(i=1;i<=fNumber;i++)      // Loop over the desired roots
945               1.26 * cofi * alpha / (1.0 + 3.5 << 928    {
946              (nwt - fAbscissa[i - 3]) / (1.0 + << 929       if(i == 1)
947     }                                          << 930       {
948     for(k = 1; k <= maxNumber; ++k)            << 931          newton = (1.0 + alpha)*(3.0 + 0.92*alpha)
949     {                                          << 932                 / (1.0 + 2.4*fNumber + 1.8*alpha) ;
950       temp1 = 1.0;                             << 933       }
951       temp2 = 0.0;                             << 934       else if(i == 2)
952                                                << 935       {
953       for(j = 1; j <= fNumber; ++j)            << 936          newton += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
954       {                                        << 937       }
955         temp3 = temp2;                         << 938       else
956         temp2 = temp1;                         << 939       {
957         temp1 =                                << 940          cofi = i - 2 ;
958           ((2 * j - 1 + alpha - nwt) * temp2 - << 941          newton += ((1.0+2.55*cofi)/(1.9*cofi)
959       }                                        << 942                     + 1.26*cofi*alpha/(1.0+3.5*cofi))
960       temp = (fNumber * temp1 - (fNumber + alp << 943                  * (newton - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
961       nwt1 = nwt;                              << 944       }
962       nwt  = nwt1 - temp1 / temp;              << 945       for(k=1;k<=maxNumber;k++)
963                                                << 946       {
964       if(std::fabs(nwt - nwt1) <= tolerance)   << 947          temp1 = 1.0 ;
965       {                                        << 948          temp2 = 0.0 ;
966         break;                                 << 
967       }                                        << 
968     }                                          << 
969     if(k > maxNumber)                          << 
970     {                                          << 
971       G4Exception("G4Integrator<T,F>::Laguerre << 
972                   "Too many (>12) iterations." << 
973     }                                          << 
974                                                << 
975     fAbscissa[i - 1] = nwt;                    << 
976     fWeight[i - 1]   = -std::exp(GammaLogarith << 
977                                GammaLogarithm( << 
978                      (temp * fNumber * temp2); << 
979   }                                            << 
980                                                << 
981   //                                           << 
982   // Integral evaluation                       << 
983   //                                           << 
984                                                   949 
985   for(i = 0; i < fNumber; i++)                 << 950          for(j=1;j<=fNumber;j++)
986   {                                            << 951          {
987     integral += fWeight[i] * (*f)(fAbscissa[i] << 952             temp3 = temp2 ;
988   }                                            << 953             temp2 = temp1 ;
989   delete[] fAbscissa;                          << 954          temp1 = ((2*j - 1 + alpha - newton)*temp2 - (j - 1 + alpha)*temp3)/j ;
990   delete[] fWeight;                            << 955          }
991   return integral;                             << 956          temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/newton ;
                                                   >> 957          newton1 = newton ;
                                                   >> 958          newton  = newton1 - temp1/temp ;
                                                   >> 959 
                                                   >> 960          if(std::fabs(newton - newton1) <= tolerance) 
                                                   >> 961          {
                                                   >> 962             break ;
                                                   >> 963          }
                                                   >> 964       }
                                                   >> 965       if(k > maxNumber)
                                                   >> 966       {
                                                   >> 967          G4Exception("Too many (>12) iterations in G4Integration::Laguerre") ;
                                                   >> 968       }
                                                   >> 969          
                                                   >> 970       fAbscissa[i-1] =  newton ;
                                                   >> 971       fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) - 
                                                   >> 972                 GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
                                                   >> 973    }
                                                   >> 974 //
                                                   >> 975 // Integral evaluation
                                                   >> 976 //
                                                   >> 977    for(i=0;i<fNumber;i++)
                                                   >> 978    {
                                                   >> 979       integral += fWeight[i]*(*f)(fAbscissa[i]) ;
                                                   >> 980    }
                                                   >> 981    delete[] fAbscissa;
                                                   >> 982    delete[] fWeight;
                                                   >> 983    return integral ;
992 }                                                 984 }
993                                                   985 
994 //////////////////////////////////////////////    986 ///////////////////////////////////////////////////////////////////////
995 //                                                987 //
996 // Auxiliary function which returns the value     988 // Auxiliary function which returns the value of std::log(gamma-function(x))
997 // Returns the value ln(Gamma(xx) for xx > 0.  << 989 // Returns the value ln(Gamma(xx) for xx > 0.  Full accuracy is obtained for 
998 // xx > 1. For 0 < xx < 1. the reflection form    990 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
999 // (Adapted from Numerical Recipes in C)          991 // (Adapted from Numerical Recipes in C)
1000 //                                               992 //
1001                                                  993 
1002 template <class T, class F>                      994 template <class T, class F>
1003 G4double G4Integrator<T, F>::GammaLogarithm(G << 995 G4double G4Integrator<T,F>::GammaLogarithm(G4double xx)
1004 {                                                996 {
1005   static const G4double cof[6] = { 76.1800917 << 997   static G4double cof[6] = { 76.18009172947146,     -86.50532032941677,
1006                                    24.0140982 << 998                              24.01409824083091,      -1.231739572450155,
1007                                    0.12086509 << 999                               0.1208650973866179e-2, -0.5395239384953e-5  } ;
1008   G4int j;                                    << 1000   register G4int j;
1009   G4double x   = xx - 1.0;                    << 1001   G4double x = xx - 1.0 ;
1010   G4double tmp = x + 5.5;                     << 1002   G4double tmp = x + 5.5 ;
1011   tmp -= (x + 0.5) * std::log(tmp);           << 1003   tmp -= (x + 0.5) * std::log(tmp) ;
1012   G4double ser = 1.000000000190015;           << 1004   G4double ser = 1.000000000190015 ;
1013                                                  1005 
1014   for(j = 0; j <= 5; ++j)                     << 1006   for ( j = 0; j <= 5; j++ )
1015   {                                              1007   {
1016     x += 1.0;                                 << 1008     x += 1.0 ;
1017     ser += cof[j] / x;                        << 1009     ser += cof[j]/x ;
1018   }                                              1010   }
1019   return -tmp + std::log(2.5066282746310005 * << 1011   return -tmp + std::log(2.5066282746310005*ser) ;
1020 }                                                1012 }
1021                                                  1013 
1022 /////////////////////////////////////////////    1014 ///////////////////////////////////////////////////////////////////////
1023 //                                               1015 //
1024 // Method involving Hermite polynomials          1016 // Method involving Hermite polynomials
1025 //                                               1017 //
1026 /////////////////////////////////////////////    1018 ///////////////////////////////////////////////////////////////////////
1027 //                                               1019 //
1028 //                                               1020 //
1029 // Gauss-Hermite method for integration of st << 1021 // Gauss-Hermite method for integration of std::exp(-x*x)*f(x) 
1030 // from minus infinity to plus infinity .     << 1022 // from minus infinity to plus infinity . 
1031 //                                               1023 //
1032                                                  1024 
1033 template <class T, class F>                   << 1025 template <class T, class F>    
1034 G4double G4Integrator<T, F>::Hermite(T& typeT << 1026 G4double G4Integrator<T,F>::Hermite( T& typeT, F f, G4int nHermite) 
1035 {                                                1027 {
1036   const G4double tolerance = 1.0e-12;         << 1028    const G4double tolerance = 1.0e-12 ;
1037   const G4int maxNumber    = 12;              << 1029    const G4int maxNumber = 12 ;
1038                                               << 1030    
1039   G4int i, j, k;                              << 1031    G4int i, j, k ;
1040   G4double integral = 0.0;                    << 1032    G4double integral = 0.0 ;
1041   G4double nwt      = 0., nwt1, temp1, temp2, << 1033    G4double newton=0., newton1, temp1, temp2, temp3, temp ;
1042                                               << 
1043   G4double piInMinusQ =                       << 
1044     std::pow(CLHEP::pi, -0.25);  // 1.0/std:: << 
1045                                                  1034 
1046   G4int fNumber       = (nHermite + 1) / 2;   << 1035    G4double piInMinusQ = std::pow(pi,-0.25) ;    // 1.0/std::sqrt(std::sqrt(pi)) ??
1047   G4double* fAbscissa = new G4double[fNumber] << 
1048   G4double* fWeight   = new G4double[fNumber] << 
1049                                                  1036 
1050   for(i = 1; i <= fNumber; ++i)               << 1037    G4int fNumber = (nHermite +1)/2 ;
1051   {                                           << 1038    G4double* fAbscissa = new G4double[fNumber] ;
1052     if(i == 1)                                << 1039    G4double* fWeight   = new G4double[fNumber] ;
1053     {                                         << 
1054       nwt = std::sqrt((G4double)(2 * nHermite << 
1055             1.85575001 * std::pow((G4double)( << 
1056     }                                         << 
1057     else if(i == 2)                           << 
1058     {                                         << 
1059       nwt -= 1.14001 * std::pow((G4double) nH << 
1060     }                                         << 
1061     else if(i == 3)                           << 
1062     {                                         << 
1063       nwt = 1.86002 * nwt - 0.86002 * fAbscis << 
1064     }                                         << 
1065     else if(i == 4)                           << 
1066     {                                         << 
1067       nwt = 1.91001 * nwt - 0.91001 * fAbscis << 
1068     }                                         << 
1069     else                                      << 
1070     {                                         << 
1071       nwt = 2.0 * nwt - fAbscissa[i - 3];     << 
1072     }                                         << 
1073     for(k = 1; k <= maxNumber; ++k)           << 
1074     {                                         << 
1075       temp1 = piInMinusQ;                     << 
1076       temp2 = 0.0;                            << 
1077                                               << 
1078       for(j = 1; j <= nHermite; ++j)          << 
1079       {                                       << 
1080         temp3 = temp2;                        << 
1081         temp2 = temp1;                        << 
1082         temp1 = nwt * std::sqrt(2.0 / j) * te << 
1083                 std::sqrt(((G4double)(j - 1)) << 
1084       }                                       << 
1085       temp = std::sqrt((G4double) 2 * nHermit << 
1086       nwt1 = nwt;                             << 
1087       nwt  = nwt1 - temp1 / temp;             << 
1088                                               << 
1089       if(std::fabs(nwt - nwt1) <= tolerance)  << 
1090       {                                       << 
1091         break;                                << 
1092       }                                       << 
1093     }                                         << 
1094     if(k > maxNumber)                         << 
1095     {                                         << 
1096       G4Exception("G4Integrator<T,F>::Hermite << 
1097                   FatalException, "Too many ( << 
1098     }                                         << 
1099     fAbscissa[i - 1] = nwt;                   << 
1100     fWeight[i - 1]   = 2.0 / (temp * temp);   << 
1101   }                                           << 
1102                                                  1040 
1103   //                                          << 1041    for(i=1;i<=fNumber;i++)
1104   // Integral calculation                     << 1042    {
1105   //                                          << 1043       if(i == 1)
                                                   >> 1044       {
                                                   >> 1045          newton = std::sqrt((G4double)(2*nHermite + 1)) - 
                                                   >> 1046                   1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ;
                                                   >> 1047       }
                                                   >> 1048       else if(i == 2)
                                                   >> 1049       {
                                                   >> 1050          newton -= 1.14001*std::pow((G4double)nHermite,0.425999)/newton ;
                                                   >> 1051       }
                                                   >> 1052       else if(i == 3)
                                                   >> 1053       {
                                                   >> 1054          newton = 1.86002*newton - 0.86002*fAbscissa[0] ;
                                                   >> 1055       }
                                                   >> 1056       else if(i == 4)
                                                   >> 1057       {
                                                   >> 1058          newton = 1.91001*newton - 0.91001*fAbscissa[1] ;
                                                   >> 1059       }
                                                   >> 1060       else 
                                                   >> 1061       {
                                                   >> 1062          newton = 2.0*newton - fAbscissa[i - 3] ;
                                                   >> 1063       }
                                                   >> 1064       for(k=1;k<=maxNumber;k++)
                                                   >> 1065       {
                                                   >> 1066          temp1 = piInMinusQ ;
                                                   >> 1067          temp2 = 0.0 ;
1106                                                  1068 
1107   for(i = 0; i < fNumber; ++i)                << 1069          for(j=1;j<=nHermite;j++)
1108   {                                           << 1070          {
1109     integral +=                               << 1071             temp3 = temp2 ;
1110       fWeight[i] * ((typeT.*f)(fAbscissa[i])  << 1072             temp2 = temp1 ;
1111   }                                           << 1073             temp1 = newton*std::sqrt(2.0/j)*temp2 - 
1112   delete[] fAbscissa;                         << 1074                     std::sqrt(((G4double)(j - 1))/j)*temp3 ;
1113   delete[] fWeight;                           << 1075          }
1114   return integral;                            << 1076          temp = std::sqrt((G4double)2*nHermite)*temp2 ;
                                                   >> 1077          newton1 = newton ;
                                                   >> 1078          newton = newton1 - temp1/temp ;
                                                   >> 1079 
                                                   >> 1080          if(std::fabs(newton - newton1) <= tolerance) 
                                                   >> 1081          {
                                                   >> 1082             break ;
                                                   >> 1083          }
                                                   >> 1084       }
                                                   >> 1085       if(k > maxNumber)
                                                   >> 1086       {
                                                   >> 1087          G4Exception("Too many (>12) iterations in G4Integrator<T,F>::Hermite") ;
                                                   >> 1088       }
                                                   >> 1089       fAbscissa[i-1] =  newton ;
                                                   >> 1090       fWeight[i-1] = 2.0/(temp*temp) ;
                                                   >> 1091    }
                                                   >> 1092 //
                                                   >> 1093 // Integral calculation
                                                   >> 1094 //
                                                   >> 1095    for(i=0;i<fNumber;i++)
                                                   >> 1096    {
                                                   >> 1097      integral += fWeight[i]*( (typeT.*f)(fAbscissa[i]) + 
                                                   >> 1098                               (typeT.*f)(-fAbscissa[i])   ) ;
                                                   >> 1099    }
                                                   >> 1100    delete[] fAbscissa;
                                                   >> 1101    delete[] fWeight;
                                                   >> 1102    return integral ;
1115 }                                                1103 }
1116                                                  1104 
                                                   >> 1105 
1117 /////////////////////////////////////////////    1106 ////////////////////////////////////////////////////////////////////////
1118 //                                               1107 //
1119 // For use with 'this' pointer                   1108 // For use with 'this' pointer
1120                                                  1109 
1121 template <class T, class F>                   << 1110 template <class T, class F>    
1122 G4double G4Integrator<T, F>::Hermite(T* ptrT, << 1111 G4double G4Integrator<T,F>::Hermite( T* ptrT, F f, G4int n)
1123 {                                                1112 {
1124   return Hermite(*ptrT, f, n);                << 1113   return Hermite(*ptrT,f,n) ;
1125 }                                             << 1114 } 
1126                                                  1115 
1127 /////////////////////////////////////////////    1116 ////////////////////////////////////////////////////////////////////////
1128 //                                               1117 //
1129 // For use with global scope f                   1118 // For use with global scope f
1130                                                  1119 
1131 template <class T, class F>                      1120 template <class T, class F>
1132 G4double G4Integrator<T, F>::Hermite(G4double << 1121 G4double G4Integrator<T,F>::Hermite( G4double (*f)(G4double), G4int nHermite) 
1133 {                                                1122 {
1134   const G4double tolerance = 1.0e-12;         << 1123    const G4double tolerance = 1.0e-12 ;
1135   const G4int maxNumber    = 12;              << 1124    const G4int maxNumber = 12 ;
1136                                               << 1125    
1137   G4int i, j, k;                              << 1126    G4int i, j, k ;
1138   G4double integral = 0.0;                    << 1127    G4double integral = 0.0 ;
1139   G4double nwt      = 0., nwt1, temp1, temp2, << 1128    G4double newton=0., newton1, temp1, temp2, temp3, temp ;
1140                                               << 1129 
1141   G4double piInMinusQ =                       << 1130    G4double piInMinusQ = std::pow(pi,-0.25) ;    // 1.0/std::sqrt(std::sqrt(pi)) ??
1142     std::pow(CLHEP::pi, -0.25);  // 1.0/std:: << 1131 
1143                                               << 1132    G4int fNumber = (nHermite +1)/2 ;
1144   G4int fNumber       = (nHermite + 1) / 2;   << 1133    G4double* fAbscissa = new G4double[fNumber] ;
1145   G4double* fAbscissa = new G4double[fNumber] << 1134    G4double* fWeight   = new G4double[fNumber] ;
1146   G4double* fWeight   = new G4double[fNumber] << 1135 
1147                                               << 1136    for(i=1;i<=fNumber;i++)
1148   for(i = 1; i <= fNumber; ++i)               << 1137    {
1149   {                                           << 1138       if(i == 1)
1150     if(i == 1)                                << 1139       {
1151     {                                         << 1140          newton = std::sqrt((G4double)(2*nHermite + 1)) - 
1152       nwt = std::sqrt((G4double)(2 * nHermite << 1141                   1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ;
1153             1.85575001 * std::pow((G4double)( << 1142       }
1154     }                                         << 1143       else if(i == 2)
1155     else if(i == 2)                           << 1144       {
1156     {                                         << 1145          newton -= 1.14001*std::pow((G4double)nHermite,0.425999)/newton ;
1157       nwt -= 1.14001 * std::pow((G4double) nH << 1146       }
1158     }                                         << 1147       else if(i == 3)
1159     else if(i == 3)                           << 1148       {
1160     {                                         << 1149          newton = 1.86002*newton - 0.86002*fAbscissa[0] ;
1161       nwt = 1.86002 * nwt - 0.86002 * fAbscis << 1150       }
1162     }                                         << 1151       else if(i == 4)
1163     else if(i == 4)                           << 1152       {
1164     {                                         << 1153          newton = 1.91001*newton - 0.91001*fAbscissa[1] ;
1165       nwt = 1.91001 * nwt - 0.91001 * fAbscis << 1154       }
1166     }                                         << 1155       else 
1167     else                                      << 1156       {
1168     {                                         << 1157          newton = 2.0*newton - fAbscissa[i - 3] ;
1169       nwt = 2.0 * nwt - fAbscissa[i - 3];     << 1158       }
1170     }                                         << 1159       for(k=1;k<=maxNumber;k++)
1171     for(k = 1; k <= maxNumber; ++k)           << 1160       {
1172     {                                         << 1161          temp1 = piInMinusQ ;
1173       temp1 = piInMinusQ;                     << 1162          temp2 = 0.0 ;
1174       temp2 = 0.0;                            << 
1175                                               << 
1176       for(j = 1; j <= nHermite; ++j)          << 
1177       {                                       << 
1178         temp3 = temp2;                        << 
1179         temp2 = temp1;                        << 
1180         temp1 = nwt * std::sqrt(2.0 / j) * te << 
1181                 std::sqrt(((G4double)(j - 1)) << 
1182       }                                       << 
1183       temp = std::sqrt((G4double) 2 * nHermit << 
1184       nwt1 = nwt;                             << 
1185       nwt  = nwt1 - temp1 / temp;             << 
1186                                               << 
1187       if(std::fabs(nwt - nwt1) <= tolerance)  << 
1188       {                                       << 
1189         break;                                << 
1190       }                                       << 
1191     }                                         << 
1192     if(k > maxNumber)                         << 
1193     {                                         << 
1194       G4Exception("G4Integrator<T,F>::Hermite << 
1195                   "Too many (>12) iterations. << 
1196     }                                         << 
1197     fAbscissa[i - 1] = nwt;                   << 
1198     fWeight[i - 1]   = 2.0 / (temp * temp);   << 
1199   }                                           << 
1200                                               << 
1201   //                                          << 
1202   // Integral calculation                     << 
1203   //                                          << 
1204                                                  1163 
1205   for(i = 0; i < fNumber; ++i)                << 1164          for(j=1;j<=nHermite;j++)
1206   {                                           << 1165          {
1207     integral += fWeight[i] * ((*f)(fAbscissa[ << 1166             temp3 = temp2 ;
1208   }                                           << 1167             temp2 = temp1 ;
1209   delete[] fAbscissa;                         << 1168             temp1 = newton*std::sqrt(2.0/j)*temp2 - 
1210   delete[] fWeight;                           << 1169                     std::sqrt(((G4double)(j - 1))/j)*temp3 ;
1211   return integral;                            << 1170          }
                                                   >> 1171          temp = std::sqrt((G4double)2*nHermite)*temp2 ;
                                                   >> 1172          newton1 = newton ;
                                                   >> 1173          newton = newton1 - temp1/temp ;
                                                   >> 1174 
                                                   >> 1175          if(std::fabs(newton - newton1) <= tolerance) 
                                                   >> 1176          {
                                                   >> 1177             break ;
                                                   >> 1178          }
                                                   >> 1179       }
                                                   >> 1180       if(k > maxNumber)
                                                   >> 1181       {
                                                   >> 1182          G4Exception("Too many (>12) iterations in G4Integrator<T,F>::Hermite") ;
                                                   >> 1183       }
                                                   >> 1184       fAbscissa[i-1] =  newton ;
                                                   >> 1185       fWeight[i-1] = 2.0/(temp*temp) ;
                                                   >> 1186    }
                                                   >> 1187 //
                                                   >> 1188 // Integral calculation
                                                   >> 1189 //
                                                   >> 1190    for(i=0;i<fNumber;i++)
                                                   >> 1191    {
                                                   >> 1192      integral += fWeight[i]*( (*f)(fAbscissa[i]) + (*f)(-fAbscissa[i])   ) ;
                                                   >> 1193    }
                                                   >> 1194    delete[] fAbscissa;
                                                   >> 1195    delete[] fWeight;
                                                   >> 1196    return integral ;
1212 }                                                1197 }
1213                                                  1198 
1214 /////////////////////////////////////////////    1199 ////////////////////////////////////////////////////////////////////////////
1215 //                                               1200 //
1216 // Method involving Jacobi polynomials           1201 // Method involving Jacobi polynomials
1217 //                                               1202 //
1218 /////////////////////////////////////////////    1203 ////////////////////////////////////////////////////////////////////////////
1219 //                                               1204 //
1220 // Gauss-Jacobi method for integration of ((1    1205 // Gauss-Jacobi method for integration of ((1-x)^alpha)*((1+x)^beta)*f(x)
1221 // from minus unit to plus unit .                1206 // from minus unit to plus unit .
1222 //                                               1207 //
1223                                                  1208 
1224 template <class T, class F>                   << 1209 template <class T, class F> 
1225 G4double G4Integrator<T, F>::Jacobi(T& typeT, << 1210 G4double G4Integrator<T,F>::Jacobi( T& typeT, F f, G4double alpha, 
1226                                     G4double  << 1211                                               G4double beta, G4int nJacobi) 
1227 {                                             << 1212 {
1228   const G4double tolerance = 1.0e-12;         << 1213   const G4double tolerance = 1.0e-12 ;
1229   const G4double maxNumber = 12;              << 1214   const G4double maxNumber = 12 ;
1230   G4int i, k, j;                              << 1215   G4int i, k, j ;
1231   G4double alphaBeta, alphaReduced, betaReduc << 1216   G4double alphaBeta, alphaReduced, betaReduced, root1=0., root2=0., root3=0. ;
1232                                               << 1217   G4double a, b, c, newton1, newton2, newton3, newton, temp, root=0., rootTemp ;
1233   G4double a, b, c, nwt1, nwt2, nwt3, nwt, te << 1218 
1234                                               << 1219   G4int     fNumber   = nJacobi ;
1235   G4int fNumber       = nJacobi;              << 1220   G4double* fAbscissa = new G4double[fNumber] ;
1236   G4double* fAbscissa = new G4double[fNumber] << 1221   G4double* fWeight   = new G4double[fNumber] ;
1237   G4double* fWeight   = new G4double[fNumber] << 1222 
1238                                               << 1223   for (i=1;i<=nJacobi;i++)
1239   for(i = 1; i <= nJacobi; ++i)               << 1224   {
1240   {                                           << 1225      if (i == 1)
1241     if(i == 1)                                << 1226      {
1242     {                                         << 1227         alphaReduced = alpha/nJacobi ;
1243       alphaReduced = alpha / nJacobi;         << 1228         betaReduced = beta/nJacobi ;
1244       betaReduced  = beta / nJacobi;          << 1229         root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
1245       root1        = (1.0 + alpha) * (2.78002 << 1230               0.767999*alphaReduced/nJacobi) ;
1246                                0.767999 * alp << 1231         root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced +
1247       root2        = 1.0 + 1.48 * alphaReduce << 1232                 0.451998*alphaReduced*alphaReduced +
1248               0.451998 * alphaReduced * alpha << 1233                 0.83001*alphaReduced*betaReduced      ;
1249               0.83001 * alphaReduced * betaRe << 1234         root  = 1.0-root1/root2 ;
1250       root = 1.0 - root1 / root2;             << 1235      } 
1251     }                                         << 1236      else if (i == 2)
1252     else if(i == 2)                           << 1237      {
1253     {                                         << 1238         root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
1254       root1 = (4.1002 + alpha) / ((1.0 + alph << 1239         root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
1255       root2 = 1.0 + 0.06 * (nJacobi - 8.0) *  << 1240         root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
1256       root3 =                                 << 1241         root -= (1.0-root)*root1*root2*root3 ;
1257         1.0 + 0.012002 * beta * (1.0 + 0.2499 << 1242      } 
1258       root -= (1.0 - root) * root1 * root2 *  << 1243      else if (i == 3) 
1259     }                                         << 1244      {
1260     else if(i == 3)                           << 1245         root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
1261     {                                         << 1246         root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
1262       root1 = (1.67001 + 0.27998 * alpha) / ( << 1247         root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
1263       root2 = 1.0 + 0.22 * (nJacobi - 8.0) /  << 1248         root -= (fAbscissa[0]-root)*root1*root2*root3 ;
1264       root3 = 1.0 + 8.0 * beta / ((6.28001 +  << 1249      }
1265       root -= (fAbscissa[0] - root) * root1 * << 1250      else if (i == nJacobi-1)
1266     }                                         << 1251      {
1267     else if(i == nJacobi - 1)                 << 1252         root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
1268     {                                         << 1253         root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
1269       root1 = (1.0 + 0.235002 * beta) / (0.76 << 1254         root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
1270       root2 = 1.0 / (1.0 + 0.639002 * (nJacob << 1255         root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
1271                              (1.0 + 0.71001 * << 1256      } 
1272       root3 = 1.0 / (1.0 + 20.0 * alpha / ((7 << 1257      else if (i == nJacobi) 
1273       root += (root - fAbscissa[nJacobi - 4]) << 1258      {
1274     }                                         << 1259         root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
1275     else if(i == nJacobi)                     << 1260         root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
1276     {                                         << 1261         root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
1277       root1 = (1.0 + 0.37002 * beta) / (1.670 << 1262         root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
1278       root2 = 1.0 / (1.0 + 0.22 * (nJacobi -  << 1263      } 
1279       root3 =                                 << 1264      else
1280         1.0 / (1.0 + 8.0 * alpha / ((6.28002  << 1265      {
1281       root += (root - fAbscissa[nJacobi - 3]) << 1266         root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
1282     }                                         << 1267      }
1283     else                                      << 1268      alphaBeta = alpha + beta ;
1284     {                                         << 1269      for (k=1;k<=maxNumber;k++)
1285       root = 3.0 * fAbscissa[i - 2] - 3.0 * f << 1270      {
1286     }                                         << 1271         temp = 2.0 + alphaBeta ;
1287     alphaBeta = alpha + beta;                 << 1272         newton1 = (alpha-beta+temp*root)/2.0 ;
1288     for(k = 1; k <= maxNumber; ++k)           << 1273         newton2 = 1.0 ;
1289     {                                         << 1274         for (j=2;j<=nJacobi;j++)
1290       temp = 2.0 + alphaBeta;                 << 1275         {
1291       nwt1 = (alpha - beta + temp * root) / 2 << 1276            newton3 = newton2 ;
1292       nwt2 = 1.0;                             << 1277            newton2 = newton1 ;
1293       for(j = 2; j <= nJacobi; ++j)           << 1278            temp = 2*j+alphaBeta ;
1294       {                                       << 1279            a = 2*j*(j+alphaBeta)*(temp-2.0) ;
1295         nwt3 = nwt2;                          << 1280             b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
1296         nwt2 = nwt1;                          << 1281            c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
1297         temp = 2 * j + alphaBeta;             << 1282            newton1 = (b*newton2-c*newton3)/a ;
1298         a    = 2 * j * (j + alphaBeta) * (tem << 1283         }
1299         b    = (temp - 1.0) *                 << 1284         newton = (nJacobi*(alpha - beta - temp*root)*newton1 +
1300             (alpha * alpha - beta * beta + te << 1285               2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2)/
1301         c    = 2.0 * (j - 1 + alpha) * (j - 1 << 1286              (temp*(1.0 - root*root)) ;
1302         nwt1 = (b * nwt2 - c * nwt3) / a;     << 1287         rootTemp = root ;
1303       }                                       << 1288         root = rootTemp - newton1/newton ;
1304       nwt = (nJacobi * (alpha - beta - temp * << 1289         if (std::fabs(root-rootTemp) <= tolerance)
1305              2.0 * (nJacobi + alpha) * (nJaco << 1290         {
1306             (temp * (1.0 - root * root));     << 1291            break ;
1307       rootTemp = root;                        << 1292         }
1308       root     = rootTemp - nwt1 / nwt;       << 1293      }
1309       if(std::fabs(root - rootTemp) <= tolera << 1294      if (k > maxNumber) 
1310       {                                       << 1295      {
1311         break;                                << 1296         G4Exception("Too many iterations (>12) in G4Integrator<T,F>::Jacobi") ;
1312       }                                       << 1297      }
1313     }                                         << 1298      fAbscissa[i-1] = root ;
1314     if(k > maxNumber)                         << 1299      fWeight[i-1] = std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) + 
1315     {                                         << 1300                         GammaLogarithm((G4double)(beta+nJacobi)) - 
1316       G4Exception("G4Integrator<T,F>::Jacobi( << 1301                         GammaLogarithm((G4double)(nJacobi+1.0)) -
1317                   FatalException, "Too many ( << 1302                         GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
1318     }                                         << 1303                         *temp*std::pow(2.0,alphaBeta)/(newton*newton2)             ;
1319     fAbscissa[i - 1] = root;                  << 1304   }
1320     fWeight[i - 1] =                          << 1305 //
1321       std::exp(GammaLogarithm((G4double)(alph << 1306 // Calculation of the integral
1322                GammaLogarithm((G4double)(beta << 1307 //
1323                GammaLogarithm((G4double)(nJac << 1308    G4double integral = 0.0 ;
1324                GammaLogarithm((G4double)(nJac << 1309    for(i=0;i<fNumber;i++)
1325       temp * std::pow(2.0, alphaBeta) / (nwt  << 1310    {
1326   }                                           << 1311       integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ;
1327                                               << 1312    }
1328   //                                          << 1313    delete[] fAbscissa;
1329   // Calculation of the integral              << 1314    delete[] fWeight;
1330   //                                          << 1315    return integral ;
1331                                               << 
1332   G4double integral = 0.0;                    << 
1333   for(i = 0; i < fNumber; ++i)                << 
1334   {                                           << 
1335     integral += fWeight[i] * (typeT.*f)(fAbsc << 
1336   }                                           << 
1337   delete[] fAbscissa;                         << 
1338   delete[] fWeight;                           << 
1339   return integral;                            << 
1340 }                                                1316 }
1341                                                  1317 
                                                   >> 1318 
1342 /////////////////////////////////////////////    1319 /////////////////////////////////////////////////////////////////////////
1343 //                                               1320 //
1344 // For use with 'this' pointer                   1321 // For use with 'this' pointer
1345                                                  1322 
1346 template <class T, class F>                   << 1323 template <class T, class F>    
1347 G4double G4Integrator<T, F>::Jacobi(T* ptrT,  << 1324 G4double G4Integrator<T,F>::Jacobi( T* ptrT, F f, G4double alpha, 
1348                                     G4int n)  << 1325                                              G4double beta, G4int n)
1349 {                                                1326 {
1350   return Jacobi(*ptrT, f, alpha, beta, n);    << 1327   return Jacobi(*ptrT,f,alpha,beta,n) ;
1351 }                                             << 1328 } 
1352                                                  1329 
1353 /////////////////////////////////////////////    1330 /////////////////////////////////////////////////////////////////////////
1354 //                                               1331 //
1355 // For use with global scope f                << 1332 // For use with global scope f 
1356                                                  1333 
1357 template <class T, class F>                      1334 template <class T, class F>
1358 G4double G4Integrator<T, F>::Jacobi(G4double  << 1335 G4double G4Integrator<T,F>::Jacobi( G4double (*f)(G4double), G4double alpha, 
1359                                     G4double  << 1336                                            G4double beta, G4int nJacobi) 
1360 {                                                1337 {
1361   const G4double tolerance = 1.0e-12;         << 1338   const G4double tolerance = 1.0e-12 ;
1362   const G4double maxNumber = 12;              << 1339   const G4double maxNumber = 12 ;
1363   G4int i, k, j;                              << 1340   G4int i, k, j ;
1364   G4double alphaBeta, alphaReduced, betaReduc << 1341   G4double alphaBeta, alphaReduced, betaReduced, root1=0., root2=0., root3=0. ;
1365                                               << 1342   G4double a, b, c, newton1, newton2, newton3, newton, temp, root=0., rootTemp ;
1366   G4double a, b, c, nwt1, nwt2, nwt3, nwt, te << 1343 
1367                                               << 1344   G4int     fNumber   = nJacobi ;
1368   G4int fNumber       = nJacobi;              << 1345   G4double* fAbscissa = new G4double[fNumber] ;
1369   G4double* fAbscissa = new G4double[fNumber] << 1346   G4double* fWeight   = new G4double[fNumber] ;
1370   G4double* fWeight   = new G4double[fNumber] << 1347 
1371                                               << 1348   for (i=1;i<=nJacobi;i++)
1372   for(i = 1; i <= nJacobi; ++i)               << 1349   {
1373   {                                           << 1350      if (i == 1)
1374     if(i == 1)                                << 1351      {
1375     {                                         << 1352         alphaReduced = alpha/nJacobi ;
1376       alphaReduced = alpha / nJacobi;         << 1353         betaReduced = beta/nJacobi ;
1377       betaReduced  = beta / nJacobi;          << 1354         root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
1378       root1        = (1.0 + alpha) * (2.78002 << 1355               0.767999*alphaReduced/nJacobi) ;
1379                                0.767999 * alp << 1356         root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced +
1380       root2        = 1.0 + 1.48 * alphaReduce << 1357                 0.451998*alphaReduced*alphaReduced +
1381               0.451998 * alphaReduced * alpha << 1358                 0.83001*alphaReduced*betaReduced      ;
1382               0.83001 * alphaReduced * betaRe << 1359         root  = 1.0-root1/root2 ;
1383       root = 1.0 - root1 / root2;             << 1360      } 
1384     }                                         << 1361      else if (i == 2)
1385     else if(i == 2)                           << 1362      {
1386     {                                         << 1363         root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
1387       root1 = (4.1002 + alpha) / ((1.0 + alph << 1364         root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
1388       root2 = 1.0 + 0.06 * (nJacobi - 8.0) *  << 1365         root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
1389       root3 =                                 << 1366         root -= (1.0-root)*root1*root2*root3 ;
1390         1.0 + 0.012002 * beta * (1.0 + 0.2499 << 1367      } 
1391       root -= (1.0 - root) * root1 * root2 *  << 1368      else if (i == 3) 
1392     }                                         << 1369      {
1393     else if(i == 3)                           << 1370         root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
1394     {                                         << 1371         root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
1395       root1 = (1.67001 + 0.27998 * alpha) / ( << 1372         root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
1396       root2 = 1.0 + 0.22 * (nJacobi - 8.0) /  << 1373         root -= (fAbscissa[0]-root)*root1*root2*root3 ;
1397       root3 = 1.0 + 8.0 * beta / ((6.28001 +  << 1374      }
1398       root -= (fAbscissa[0] - root) * root1 * << 1375      else if (i == nJacobi-1)
1399     }                                         << 1376      {
1400     else if(i == nJacobi - 1)                 << 1377         root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
1401     {                                         << 1378         root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
1402       root1 = (1.0 + 0.235002 * beta) / (0.76 << 1379         root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
1403       root2 = 1.0 / (1.0 + 0.639002 * (nJacob << 1380         root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
1404                              (1.0 + 0.71001 * << 1381      } 
1405       root3 = 1.0 / (1.0 + 20.0 * alpha / ((7 << 1382      else if (i == nJacobi) 
1406       root += (root - fAbscissa[nJacobi - 4]) << 1383      {
1407     }                                         << 1384         root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
1408     else if(i == nJacobi)                     << 1385         root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
1409     {                                         << 1386         root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
1410       root1 = (1.0 + 0.37002 * beta) / (1.670 << 1387         root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
1411       root2 = 1.0 / (1.0 + 0.22 * (nJacobi -  << 1388      } 
1412       root3 =                                 << 1389      else
1413         1.0 / (1.0 + 8.0 * alpha / ((6.28002  << 1390      {
1414       root += (root - fAbscissa[nJacobi - 3]) << 1391         root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
1415     }                                         << 1392      }
1416     else                                      << 1393      alphaBeta = alpha + beta ;
1417     {                                         << 1394      for (k=1;k<=maxNumber;k++)
1418       root = 3.0 * fAbscissa[i - 2] - 3.0 * f << 1395      {
1419     }                                         << 1396         temp = 2.0 + alphaBeta ;
1420     alphaBeta = alpha + beta;                 << 1397         newton1 = (alpha-beta+temp*root)/2.0 ;
1421     for(k = 1; k <= maxNumber; ++k)           << 1398         newton2 = 1.0 ;
1422     {                                         << 1399         for (j=2;j<=nJacobi;j++)
1423       temp = 2.0 + alphaBeta;                 << 1400         {
1424       nwt1 = (alpha - beta + temp * root) / 2 << 1401            newton3 = newton2 ;
1425       nwt2 = 1.0;                             << 1402            newton2 = newton1 ;
1426       for(j = 2; j <= nJacobi; ++j)           << 1403            temp = 2*j+alphaBeta ;
1427       {                                       << 1404            a = 2*j*(j+alphaBeta)*(temp-2.0) ;
1428         nwt3 = nwt2;                          << 1405            b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
1429         nwt2 = nwt1;                          << 1406            c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
1430         temp = 2 * j + alphaBeta;             << 1407            newton1 = (b*newton2-c*newton3)/a ;
1431         a    = 2 * j * (j + alphaBeta) * (tem << 1408         }
1432         b    = (temp - 1.0) *                 << 1409         newton = (nJacobi*(alpha - beta - temp*root)*newton1 +
1433             (alpha * alpha - beta * beta + te << 1410              2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2) /
1434         c    = 2.0 * (j - 1 + alpha) * (j - 1 << 1411              (temp*(1.0 - root*root)) ;
1435         nwt1 = (b * nwt2 - c * nwt3) / a;     << 1412         rootTemp = root ;
1436       }                                       << 1413         root = rootTemp - newton1/newton ;
1437       nwt = (nJacobi * (alpha - beta - temp * << 1414         if (std::fabs(root-rootTemp) <= tolerance)
1438              2.0 * (nJacobi + alpha) * (nJaco << 1415         {
1439             (temp * (1.0 - root * root));     << 1416            break ;
1440       rootTemp = root;                        << 1417         }
1441       root     = rootTemp - nwt1 / nwt;       << 1418      }
1442       if(std::fabs(root - rootTemp) <= tolera << 1419      if (k > maxNumber) 
1443       {                                       << 1420      {
1444         break;                                << 1421         G4Exception("Too many iterations (>12) in G4Integrator<T,F>::Jacobi") ;
1445       }                                       << 1422      }
1446     }                                         << 1423      fAbscissa[i-1] = root ;
1447     if(k > maxNumber)                         << 1424      fWeight[i-1] =
1448     {                                         << 1425         std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) + 
1449       G4Exception("G4Integrator<T,F>::Jacobi( << 1426                  GammaLogarithm((G4double)(beta+nJacobi)) - 
1450                   "Too many (>12) iterations. << 1427                  GammaLogarithm((G4double)(nJacobi+1.0)) -
1451     }                                         << 1428                  GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
1452     fAbscissa[i - 1] = root;                  << 1429         *temp*std::pow(2.0,alphaBeta)/(newton*newton2);
1453     fWeight[i - 1] =                          << 1430   }
1454       std::exp(GammaLogarithm((G4double)(alph << 1431 //
1455                GammaLogarithm((G4double)(beta << 1432 // Calculation of the integral
1456                GammaLogarithm((G4double)(nJac << 1433 //
1457                GammaLogarithm((G4double)(nJac << 1434    G4double integral = 0.0 ;
1458       temp * std::pow(2.0, alphaBeta) / (nwt  << 1435    for(i=0;i<fNumber;i++)
1459   }                                           << 1436    {
1460                                               << 1437       integral += fWeight[i]*(*f)(fAbscissa[i]) ;
1461   //                                          << 1438    }
1462   // Calculation of the integral              << 1439    delete[] fAbscissa;
1463   //                                          << 1440    delete[] fWeight;
1464                                               << 1441    return integral ;
1465   G4double integral = 0.0;                    << 
1466   for(i = 0; i < fNumber; ++i)                << 
1467   {                                           << 
1468     integral += fWeight[i] * (*f)(fAbscissa[i << 
1469   }                                           << 
1470   delete[] fAbscissa;                         << 
1471   delete[] fWeight;                           << 
1472   return integral;                            << 
1473 }                                                1442 }
1474                                                  1443 
1475 //                                               1444 //
1476 //                                               1445 //
1477 /////////////////////////////////////////////    1446 ///////////////////////////////////////////////////////////////////
1478                                                  1447