Geant4_10
G4GaussHermiteQ.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4GaussHermiteQ.cc 67970 2013-03-13 10:10:06Z gcosmo $
28 //
29 
30 #include "G4GaussHermiteQ.hh"
31 #include "G4PhysicalConstants.hh"
32 
33 // ----------------------------------------------------------
34 //
35 // Constructor for Gauss-Hermite
36 
37 G4GaussHermiteQ::G4GaussHermiteQ ( function pFunction,
38  G4int nHermite )
39  : G4VGaussianQuadrature(pFunction)
40 {
41  const G4double tolerance = 1.0e-12 ;
42  const G4int maxNumber = 12 ;
43 
44  G4int i=1, j=1, k=1 ;
45  G4double newton0=0.;
46  G4double newton1=0.0, temp1=0.0, temp2=0.0, temp3=0.0, temp=0.0 ;
47  G4double piInMinusQ = std::pow(pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ??
48 
49  fNumber = (nHermite +1)/2 ;
50  fAbscissa = new G4double[fNumber] ;
51  fWeight = new G4double[fNumber] ;
52 
53  for(i=1;i<=fNumber;i++)
54  {
55  if(i == 1)
56  {
57  newton0 = std::sqrt((G4double)(2*nHermite + 1)) -
58  1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ;
59  }
60  else if(i == 2)
61  {
62  newton0 -= 1.14001*std::pow((G4double)nHermite,0.425999)/newton0 ;
63  }
64  else if(i == 3)
65  {
66  newton0 = 1.86002*newton0 - 0.86002*fAbscissa[0] ;
67  }
68  else if(i == 4)
69  {
70  newton0 = 1.91001*newton0 - 0.91001*fAbscissa[1] ;
71  }
72  else
73  {
74  newton0 = 2.0*newton0 - fAbscissa[i - 3] ;
75  }
76  for(k=1;k<=maxNumber;k++)
77  {
78  temp1 = piInMinusQ ;
79  temp2 = 0.0 ;
80  for(j=1;j<=nHermite;j++)
81  {
82  temp3 = temp2 ;
83  temp2 = temp1 ;
84  temp1 = newton0*std::sqrt(2.0/j)*temp2
85  - std::sqrt(((G4double)(j - 1))/j)*temp3 ;
86  }
87  temp = std::sqrt((G4double)2*nHermite)*temp2 ;
88  newton1 = newton0 ;
89  newton0 = newton1 - temp1/temp ;
90  if(std::fabs(newton0 - newton1) <= tolerance)
91  {
92  break ;
93  }
94  }
95  if(k > maxNumber)
96  {
97  G4Exception("G4GaussHermiteQ::G4GaussHermiteQ()",
98  "OutOfRange", FatalException,
99  "Too many iterations in Gauss-Hermite constructor.") ;
100  }
101  fAbscissa[i-1] = newton0 ;
102  fWeight[i-1] = 2.0/(temp*temp) ;
103  }
104 }
105 
106 
107 // ----------------------------------------------------------
108 //
109 // Gauss-Hermite method for integration of std::exp(-x*x)*nFunction(x)
110 // from minus infinity to plus infinity .
111 
113 {
114  G4double integral = 0.0 ;
115  for(G4int i=0;i<fNumber;i++)
116  {
117  integral += fWeight[i]*(fFunction(fAbscissa[i])
118  + fFunction(-fAbscissa[i])) ;
119  }
120  return integral ;
121 }
int G4int
Definition: G4Types.hh:78
G4double Integral() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double G4double
Definition: G4Types.hh:76
G4GaussHermiteQ(function pFunction, G4int nHermite)