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 5.1)


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