85 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
99 for(
G4int i = 0; i < 500; ++i )
126 G4cout<<
"G4PAIySection::Initialize(...,G4SandiaTable* sandia)"<<
G4endl;
175 G4cout<<
"last i = "<<i<<
"; "<<
"fIntervalNumber = "
192 G4cout<<
"Now checking, if two borders are too close together"<<
G4endl;
250 static const G4double p0 = 1.20923e+00;
251 static const G4double p1 = 3.53256e-01;
252 static const G4double p2 = -1.45052e-03;
257 for( i = 0; i < numberOfElements; ++i )
259 thisMaterialZ[i] =
material->GetElement(i)->GetZ();
260 sumZ += thisMaterialZ[i];
261 thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
263 for( i = 0; i < numberOfElements; ++i )
265 sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
268 delete [] thisMaterialZ;
269 delete [] thisMaterialCof;
301 for(
G4int j = 1; j < 112; ++j)
335 for( j = 1; j <= 2; ++j)
441 G4double a = log10(y2/yy1)/log10(x2/x1);
442 G4double b = log10(yy1) - a*log10(x1);
472 if( x >
fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*
fDelta )
497 c2 = (x2 - x1)*(x2 + x1)/(x12*x12);
498 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/(x12*x12*x12);
501 return fA1[k]*log(x2/x1) +
fA2[k]*c1 +
fA3[k]*c2/2 +
fA4[k]*c3/3;
514 G4double energy2,energy3,energy4,result;
516 energy2 = energy1*energy1;
517 energy3 = energy2*energy1;
518 energy4 = energy3*energy1;
520 result =
fA1[k]/energy1+
fA2[k]/energy2+
fA3[k]/energy3+
fA4[k]/energy4;
521 result *=
hbarc/energy1;
536 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
537 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
556 xln3 = log((x2 + x0)/(x1 + x0));
563 c2 = (x2 - x1)*(x2 +x1)/(x12*x12);
564 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/(x12*x12*x12);
566 result -= (
fA1[i]/x02 +
fA3[i]/x04)*xln1;
567 result -= (
fA2[i]/x02 +
fA4[i]/x04)*c1;
568 result -=
fA3[i]*c2/2/x02;
569 result -=
fA4[i]*c3/3/x02;
571 cof1 =
fA1[i]/x02 +
fA3[i]/x04;
572 cof2 =
fA2[i]/x03 +
fA4[i]/x05;
574 result += 0.5*(cof1 +cof2)*xln2;
575 result += 0.5*(cof1 - cof2)*xln3;
592 G4double beta, be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
597 be2 = betaGammaSq/(1 + betaGammaSq);
603 if( betaGammaSq < 0.01 ) x2 = log(be2);
659 G4double logarithm, x3, x5, argument, modul2, dNdxC;
664 be2 = betaGammaSq/(1 + betaGammaSq);
667 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq);
673 logarithm += log(1+1.0/betaGammaSq);
686 if( x3 == 0.0 ) argument = 0.5*
pi;
692 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
718 G4double cof, resonance, modul2, dNdxP;
723 be2 = betaGammaSq/(1 + betaGammaSq);
850 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)
return 0.;
857 a = log10(yy1/y0)/log10(c);
864 result = b*log(x1/x0);
868 result = y0*(x1*pow(c,a-1) - x0)/a;
892 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)
return 0.;
897 a = log10(yy1/y0)/log10(c);
903 result = b*log(x1/x0);
907 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
926 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)
return 0.;
934 a = log10(yy1/y0)/log10(c);
936 if(a < 20.) b = y0/pow(x0,a);
939 if(a == 0) result = b*log(c);
940 else result = y0*(x1*pow(c,a-1) - x0)/a;
963 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)
return 0.;
968 a = log10(yy1/y0)/log10(c);
971 if(a < 20.) b = y0/pow(x0,a);
974 if(a == 0) result = b*log(x1/x0);
975 else result = y0*(x1*pow(c,a-1) - x0)/a;
993 G4double x0,x1,y0,yy1,a,d,e0,result;
1003 a = log10(yy1/y0)/log10(x1/x0);
1006 if(a < 20.) b = y0/pow(x0,a);
1011 result = b*log(x0/e0);
1015 result = y0*(x0 - e0*pow(d,a-1))/a;
1033 a = log10(yy1/y0)/log10(x1/x0);
1039 result += b*log(e0/x0);
1043 result += y0*(e0*pow(d,a-1) - x0)/a;
1063 G4double x0,x1,y0,yy1,a,d,e0,result;
1073 a = log10(yy1/y0)/log10(x1/x0);
1076 if(a < 20.) b = y0/pow(x0,a);
1081 result = b*log(x0/e0);
1085 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1094 a = log10(yy1/y0)/log10(x1/x0);
1096 if(a < 20.) b = y0/pow(x0,a);
1101 result += b*log(e0/x0);
1105 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1119 G4double x0,x1,y0,yy1,a,e0,c,d,result;
1132 a = log10(yy1/y0)/log10(c);
1135 if(a < 20.) b = y0/pow(x0,a);
1138 if( a == 0 ) result = b*log(x0/e0);
1139 else result = y0*(x0 - e0*pow(d,a-1))/a;
1157 a = log10(yy1/y0)/log10(x1/x0);
1160 if(a > 20.0) b = 0.0;
1161 else b = y0/pow(x0,a);
1166 if( a == 0 ) result += b*log(e0/x0);
1167 else result += y0*(e0*pow(d,a-1) - x0 )/a;
1188 G4double x0,x1,y0,yy1,a,c,d,e0,result;
1198 a = log10(yy1/y0)/log10(c);
1201 if(a < 20.) b = y0/pow(x0,a);
1204 if( a == 0 ) result = b*log(x0/e0);
1205 else result = y0*(x0 - e0*pow(d,a-1))/a;
1218 a = log10(yy1/y0)/log10(c);
1220 if(a < 20.) b = y0/pow(x0,a);
1223 if( a == 0 ) result += b*log(e0/x0);
1224 else result += y0*(e0*pow(d,a-1) - x0)/a;
1250 numOfCollisions =
G4Poisson(meanNumber);
1254 while(numOfCollisions)
1287 numOfCollisions =
G4Poisson(meanNumber);
1291 while(numOfCollisions)
1324 numOfCollisions =
G4Poisson(meanNumber);
1328 while(numOfCollisions)
1350 G4String head =
"G4PAIySection::" + methodName +
"()";
1352 ed <<
"Wrong index " << i <<
" fSplineNumber= " <<
fSplineNumber;
13661.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
13671.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00,
13681.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
13691.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00,
13702.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
13713.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00,
13725.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
13738.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01,
13741.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
13752.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01,
13765.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
13771.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02,
13781.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
13793.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02,
13806.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
13811.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03,
13822.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
13834.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03,
13848.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
13851.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04,
13863.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
13875.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04,
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
G4long G4Poisson(G4double mean)
static constexpr double keV
static constexpr double eV
static constexpr double pi
G4GLOB_DLL std::ostream G4cout
static const G4int fMaxSplineSize
G4DataVector fIntegralPAIdEdx
G4double GetStepCerenkovLoss(G4double step)
G4double RePartDielectricConst(G4double energy)
G4DataVector fImPartDielectricConst
static G4int fNumberOfGammas
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
G4double GetStepPlasmonLoss(G4double step)
G4double SumOverBordPlasmon(G4int intervalNumber, G4double energy)
static const G4double fDelta
void NormShift(G4double betaGammaSq)
void ComputeLowEnergyCof(const G4Material *material)
G4double SumOverBorderdEdx(G4int intervalNumber, G4double energy)
G4DataVector fdNdxCerenkov
G4double SumOverInterPlasmon(G4int intervalNumber)
G4DataVector fDifPAIySection
G4DataVector fdNdxPlasmon
G4DataVector fSplineEnergy
void IntegralPAIySection()
static const G4int fRefGammaNumber
G4DataVector fIntegralPAIySection
G4double GetStepEnergyLoss(G4double step)
static const G4double fError
G4double SumOverInterCerenkov(G4int intervalNumber)
G4DataVector fIntegralCerenkov
G4DataVector fIntegralPlasmon
G4DataVector fIntegralTerm
G4double SumOverInterval(G4int intervalNumber)
G4double SumOverBorder(G4int intervalNumber, G4double energy)
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
G4double SumOverBordCerenkov(G4int intervalNumber, G4double energy)
static const G4double fLorentzFactor[112]
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
G4double DifPAIySection(G4int intervalNumber, G4double betaGammaSq)
G4DataVector fRePartDielectricConst
void CallError(G4int i, const G4String &methodName) const
G4double SumOverIntervaldEdx(G4int intervalNumber)
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
G4double GetLorentzFactor(G4int j) const
G4double fNormalizationCof
G4DataVector fEnergyInterval
G4double fElectronDensity
G4double fPAItable[500][112]
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
void SplainPAI(G4double betaGammaSq)
G4int GetMaxInterval() const
G4double GetSandiaMatTablePAI(G4int, G4int) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments