48 G4double alphaBeta=0.0, alphaReduced=0.0, betaReduced=0.0,
49 root1=0.0, root2=0.0, root3=0.0 ;
51 newton1=0.0, newton2=0.0, newton3=0.0, newton0=0.0,
52 temp=0.0, rootTemp=0.0 ;
58 for (i=1;i<=nJacobi;i++)
62 alphaReduced = alpha/nJacobi ;
63 betaReduced = beta/nJacobi ;
64 root1 = (1.0+
alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
65 0.767999*alphaReduced/nJacobi) ;
66 root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced
67 + 0.451998*alphaReduced*alphaReduced
68 + 0.83001*alphaReduced*betaReduced ;
69 root = 1.0-root1/root2 ;
73 root1=(4.1002+
alpha)/((1.0+alpha)*(1.0+0.155998*
alpha)) ;
74 root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*
alpha)/nJacobi ;
75 root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
76 root -= (1.0-root)*root1*root2*root3 ;
80 root1=(1.67001+0.27998*
alpha)/(1.0+0.37002*alpha) ;
81 root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
82 root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
83 root -= (
fAbscissa[0]-root)*root1*root2*root3 ;
85 else if (i == nJacobi-1)
87 root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
88 root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
89 root3=1.0/(1.0+20.0*alpha/((7.5+
alpha)*nJacobi*nJacobi)) ;
90 root += (root-
fAbscissa[nJacobi-4])*root1*root2*root3 ;
92 else if (i == nJacobi)
94 root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
95 root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
96 root3 = 1.0/(1.0+8.0*alpha/((6.28002+
alpha)*nJacobi*nJacobi)) ;
97 root += (root-
fAbscissa[nJacobi-3])*root1*root2*root3 ;
103 alphaBeta = alpha + beta ;
104 for (k=1;k<=maxNumber;k++)
106 temp = 2.0 + alphaBeta ;
107 newton1 = (alpha-beta+temp*root)/2.0 ;
109 for (
G4int j=2;j<=nJacobi;j++)
113 temp = 2*j+alphaBeta ;
114 a = 2*j*(j+alphaBeta)*(temp-2.0) ;
115 b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
116 c = 2.0*(j-1+
alpha)*(j-1+beta)*temp ;
117 newton1 = (b*newton2-c*newton3)/a ;
119 newton0 = (nJacobi*(alpha - beta - temp*root)*newton1 +
120 2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2)/
121 (temp*(1.0 - root*root)) ;
123 root = rootTemp - newton1/newton0 ;
124 if (std::fabs(root-rootTemp) <= tolerance)
131 G4Exception(
"G4GaussJacobiQ::G4GaussJacobiQ()",
"OutOfRange",
139 *temp*std::pow(2.0,alphaBeta)/(newton0*newton2) ;
std::vector< ExP01TrackerHit * > a
G4GaussJacobiQ(function pFunction, G4double alpha, G4double beta, G4int nJacobi)
G4double Integral() const
G4double GammaLogarithm(G4double xx)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static const G4double alpha