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 8.3.p1)


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