137 for ( std::vector<G4PhysicsTable*>::iterator
it =
fAngleBank.begin();
139 if ( (*
it) ) (*it)->clearAndDestroy();
163 for(jEl = 0 ; jEl < numOfEl; ++jEl)
173 G4cout<<
"G4NuclNuclDiffuseElastic::Initialise() the element: "
174 <<(*theElementTable)[jEl]->GetName()<<
G4endl;
223 if (iZ == 1 && iA == 1) theDef =
theProton;
226 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
227 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
241 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
243 if( cost >= 1.0 ) cost = 1.0;
244 else if( cost <= -1.0) cost = -1.0;
246 G4double thetaCMS = std::acos(cost);
278 if( z && (kRt > kRtC) )
310 if (iZ == 1 && iA == 1) theDef =
theProton;
313 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
314 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
328 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
330 if( cost >= 1.0 ) cost = 1.0;
331 else if( cost <= -1.0) cost = -1.0;
333 G4double thetaCMS = std::acos(cost);
360 if (iZ == 1 && iA == 1) theDef =
theProton;
363 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
364 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
378 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
380 if( cost >= 1.0 ) cost = 1.0;
381 else if( cost <= -1.0) cost = -1.0;
383 G4double thetaCMS = std::acos(cost);
403 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
416 bzero2 = bzero*bzero;
420 bonebyarg2 = bonebyarg*bonebyarg;
435 diffuse = 0.63*
fermi;
464 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
465 sigma += kr2*bonebyarg2;
483 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
496 bzero2 = bzero*bzero;
500 bonebyarg2 = bonebyarg*bonebyarg;
504 diffuse = 0.63*
fermi;
513 diffuse = 0.63*
fermi;
527 G4double sinHalfTheta = std::sin(0.5*theta);
528 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
554 sigma += mode2k2*bone2;
555 sigma += e2dk3t*bzero*bone;
558 sigma += kr2*bonebyarg2;
576 theta = std::sqrt(alpha);
580 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
593 bzero2 = bzero*bzero;
597 bonebyarg2 = bonebyarg*bonebyarg;
601 diffuse = 0.63*
fermi;
610 diffuse = 0.63*
fermi;
625 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
651 sigma += mode2k2*bone2;
652 sigma += e2dk3t*bzero*bone;
655 sigma += kr2*bonebyarg2;
713 G4double t = 2*p*p*( 1 - std::cos(theta) );
737 if (thetaMax >
pi) thetaMax =
pi;
746 for(i = 1; i <= iMax; i++)
748 theta1 = (i-1)*thetaMax/iMax;
749 theta2 = i*thetaMax/iMax;
754 result = 0.5*(theta1 + theta2);
758 if (i > iMax ) result = 0.5*(theta1 + theta2);
764 if(result < 0.) result = 0.;
765 if(result > thetaMax) result = thetaMax;
785 G4double totElab = std::sqrt(m1*m1+p*p);
859 G4int iMomentum, iAngle;
881 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
883 for( iMomentum = 0; iMomentum <
fEnergyBin; iMomentum++)
887 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) )
break;
892 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;
893 if ( iMomentum < 0 ) iMomentum = 0;
896 if (iMomentum == fEnergyBin -1 || iMomentum == 0 )
902 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
904 if( position < (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
906 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
921 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
924 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
926 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
944 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
947 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
949 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
963 randAngle = W1*theta1 + W2*theta2;
990 G4cout<<
"G4NuclNuclDiffuseElastic::Initialise() the element with Z = "
991 <<Z<<
"; and A = "<<A<<
G4endl;
1026 partMom = std::sqrt( kinE*(kinE + 2*m1) );
1032 if(alphaMax >
pi) alphaMax =
pi;
1061 alpha1 = alphaCoulomb + delth*(j-1);
1063 alpha2 = alpha1 + delth;
1092 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1099 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1101 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1102 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1104 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1105 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1107 if ( x1 == x2 ) randAngle =
x2;
1113 randAngle = x1 + ( position -
y1 )*( x2 - x1 )/( y2 -
y1 );
1152 t =
SampleT( theParticle, ptot, A);
1155 if(!(t < 0.0 || t >= 0.0))
1159 G4cout <<
"G4NuclNuclDiffuseElastic:WARNING: A = " << A
1160 <<
" mom(GeV)= " << plab/
GeV
1161 <<
" S-wave will be sampled"
1168 G4cout <<
" t= " << t <<
" tmax= " << tmax
1169 <<
" ptot= " << ptot <<
G4endl;
1182 else if( cost <= -1.0)
1189 sint = std::sqrt((1.0-cost)*(1.0+cost));
1193 G4cout <<
"cos(t)=" << cost <<
" std::sin(t)=" << sint <<
G4endl;
1236 G4double cost = std::cos(thetaCMS);
1244 else if( cost <= -1.0)
1251 sint = std::sqrt((1.0-cost)*(1.0+cost));
1255 G4cout <<
"cos(tcms)=" << cost <<
" std::sin(tcms)=" << sint <<
G4endl;
1297 G4double cost = std::cos(thetaLab);
1305 else if( cost <= -1.0)
1312 sint = std::sqrt((1.0-cost)*(1.0+cost));
1316 G4cout <<
"cos(tlab)=" << cost <<
" std::sin(tlab)=" << sint <<
G4endl;
1346 G4cout<<
"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = "
1347 <<Z<<
"; and A = "<<A<<
G4endl;
1356 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1357 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1358 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1372 alphaMax = kRmax*kRmax/kR2;
1374 if (alphaMax > 4.) alphaMax = 4.;
1376 alphaCoulomb = kRcoul*kRcoul/kR2;
1381 fBeta = a/std::sqrt(1+a*a);
1397 alpha1 = alphaMax*(j-1)/fAngleBin;
1398 alpha2 = alphaMax*( j )/fAngleBin;
1400 if( ( alpha2 > alphaCoulomb ) &&
z )
fAddCoulomb =
true;
1405 alpha1, alpha2,epsilon);
1415 <<sumL10<<
"\t"<<sumL96<<
"\t"<<sumAG<<
G4endl;
1417 angleVector->
PutValue( j-1 , alpha1, sumL10 );
1445 if ( n < 0 ) legPol = 0.;
1446 else if( n == 0 ) legPol = 1.;
1447 else if( n == 1 ) legPol =
x;
1448 else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
1449 else if( n == 3 ) legPol = (5.*x*x*x-3.*
x)/2.;
1450 else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
1451 else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*
x)/8.;
1452 else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
1457 legPol = std::sqrt( 2./(n*
CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*
CLHEP::pi );
1469 G4double n2, cofn, shny, chny, fn, gn;
1488 for( n = 1; n <= nMax; n++)
1492 cofn =
G4Exp(-0.5*n2)/(n2+twox2);
1494 chny = std::cosh(n*y);
1495 shny = std::sinh(n*y);
1497 fn = twox - twoxcos2xy*chny + n*sin2xy*shny;
1498 gn = twoxsin2xy*chny + n*cos2xy*shny;
1516 outRe +=
GetErf(x) + cof1*(1-cos2xy)/twox;
1517 outIm += cof1*sin2xy/twox;
1567 order /= std::sqrt(2.);
1570 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1571 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1572 G4complex out = gamma*(1. - a1*dTheta) - a0;
1595 order /= std::sqrt(2.);
1597 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1598 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1599 G4complex out = -gamma*(1. - a1*dTheta) - a0;
1637 G4double sindTheta = std::sin(dTheta);
1638 G4double persqrt2 = std::sqrt(0.5);
1671 for( n = 0; n <
fMaxL; n++)
1678 shiftN = std::exp( -0.5*(1.-im*
fEtaRatio)*T12b ) - 1.;
1694 G4double T12b,
a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
1695 G4double sinThetaH2 = sinThetaH*sinThetaH;
1704 for( n = 1; n <
fMaxL; n++)
1706 T12b = aTemp*
G4Exp(-b2/n)/
n;
1746 fBeta = a/std::sqrt(1+a*a);
1781 fBeta = a/std::sqrt(1+a*a);
1821 if( pN < 0. ) pN = 0.;
1824 if( tN < 0. ) tN = 0.;
1843 fBeta = a/std::sqrt(1+a*a);
1874 G4double proj_energy = proj_mass + pTkin;
1875 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
1890 if( proj_momentum >= 1.2 )
1894 else if( proj_momentum >= 0.6 )
1901 fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
1907 if( proj_momentum >= 10. )
1917 A0 = 100. - B0*
G4Log(3.0e7);
1919 xsection = A0 + B0*
G4Log(proj_energy) - 11
1921 0.93827*0.93827,-0.165);
1926 if(pParticle == tParticle)
1928 if( proj_momentum < 0.73 )
1932 else if( proj_momentum < 1.05 )
1934 hnXsc = 23 + 40*(
G4Log(proj_momentum/0.73))*
1935 (
G4Log(proj_momentum/0.73));
1946 if( proj_momentum < 0.8 )
1950 else if( proj_momentum < 1.4 )
1976 G4double sindTheta = std::sin(dTheta);
1993 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1994 out += ( cosFresnel + sinFresnel - 1. )*prof;
1998 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
2011 const G4double cof[6] = { 76.18009172947146, -86.50532032941677,
2012 24.01409824083091, -1.231739572450155,
2013 0.1208650973866179e-2, -0.5395239384953e-5 } ;
2017 tmp -= (z + 0.5) * std::log(tmp);
2020 for ( j = 0; j <= 5; j++ )
2025 return -tmp + std::log(2.5066282746310005*ser);
2035 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2037 modvalue = std::fabs(value);
2041 value2 = value*
value;
2043 fact1 = 57568490574.0 + value2*(-13362590354.0
2044 + value2*( 651619640.7
2045 + value2*(-11214424.18
2046 + value2*( 77392.33017
2047 + value2*(-184.9052456 ) ) ) ) );
2049 fact2 = 57568490411.0 + value2*( 1029532985.0
2050 + value2*( 9494680.718
2051 + value2*(59272.64853
2052 + value2*(267.8532712
2053 + value2*1.0 ) ) ) );
2055 bessel = fact1/fact2;
2063 shift = modvalue-0.785398164;
2065 fact1 = 1.0 + value2*(-0.1098628627e-2
2066 + value2*(0.2734510407e-4
2067 + value2*(-0.2073370639e-5
2068 + value2*0.2093887211e-6 ) ) );
2070 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
2071 + value2*(-0.6911147651e-5
2072 + value2*(0.7621095161e-6
2073 - value2*0.934945152e-7 ) ) );
2075 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
2087 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2089 modvalue = std::fabs(value);
2091 if ( modvalue < 8.0 )
2093 value2 = value*
value;
2095 fact1 = value*(72362614232.0 + value2*(-7895059235.0
2096 + value2*( 242396853.1
2097 + value2*(-2972611.439
2098 + value2*( 15704.48260
2099 + value2*(-30.16036606 ) ) ) ) ) );
2101 fact2 = 144725228442.0 + value2*(2300535178.0
2102 + value2*(18583304.74
2103 + value2*(99447.43394
2104 + value2*(376.9991397
2105 + value2*1.0 ) ) ) );
2106 bessel = fact1/fact2;
2114 shift = modvalue - 2.356194491;
2116 fact1 = 1.0 + value2*( 0.183105e-2
2117 + value2*(-0.3516396496e-4
2118 + value2*(0.2457520174e-5
2119 + value2*(-0.240337019e-6 ) ) ) );
2121 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
2122 + value2*( 0.8449199096e-5
2123 + value2*(-0.88228987e-6
2124 + value2*0.105787412e-6 ) ) );
2126 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
2128 if (value < 0.0) bessel = -bessel;