63 const G4double G4PAIySection::fDelta = 0.005;
64 const G4double G4PAIySection::fError = 0.005;
66 const G4int G4PAIySection::fMaxSplineSize = 500;
70 static const G4double cofBetaBohr = 4.0;
72 static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
82 fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
83 fIntervalNumber = fSplineNumber = 0;
87 fRePartDielectricConst =
G4DataVector(fMaxSplineSize,0.0);
88 fImPartDielectricConst =
G4DataVector(fMaxSplineSize,0.0);
98 for(
G4int i = 0; i < 500; ++i )
100 for(
G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
123 G4cout<<
"G4PAIySection::Initialize(...,G4SandiaTable* sandia)"<<
G4endl;
137 G4cout<<
"fDensity = "<<fDensity<<
"\t"<<fElectronDensity<<
"\t fIntervalNumber = "<<fIntervalNumber<<
G4endl;
145 for( i = 1; i <= fIntervalNumber; i++ )
154 fEnergyInterval[i] = maxEnergyTransfer;
166 G4cout<<i<<
"\t"<<fEnergyInterval[i]/
keV<<
"\t"<<fA1[i]<<
"\t"<<fA2[i]<<
"\t"
167 <<fA3[i]<<
"\t"<<fA4[i]<<
"\t"<<
G4endl;
170 if( fVerbose > 0 )
G4cout<<
"last i = "<<i<<
"; "<<
"fIntervalNumber = "<<fIntervalNumber<<
G4endl;
172 if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer )
175 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
179 for( i = 1; i <= fIntervalNumber; i++ )
181 G4cout<<i<<
"\t"<<fEnergyInterval[i]/
keV<<
"\t"<<fA1[i]<<
"\t"<<fA2[i]<<
"\t"
182 <<fA3[i]<<
"\t"<<fA4[i]<<
"\t"<<
G4endl;
185 if( fVerbose > 0 )
G4cout<<
"Now checking, if two borders are too close together"<<
G4endl;
187 for( i = 1; i < fIntervalNumber; i++ )
189 if( fEnergyInterval[i+1]-fEnergyInterval[i] >
190 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]) )
continue;
193 for( j = i; j < fIntervalNumber; j++ )
195 fEnergyInterval[j] = fEnergyInterval[j+1];
207 for( i = 1; i <= fIntervalNumber; i++ )
209 G4cout<<i<<
"\t"<<fEnergyInterval[i]/
keV<<
"\t"<<fA1[i]<<
"\t"<<fA2[i]<<
"\t"
210 <<fA3[i]<<
"\t"<<fA4[i]<<
"\t"<<
G4endl;
215 ComputeLowEnergyCof(material);
218 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
220 NormShift(betaGammaSqRef);
221 SplainPAI(betaGammaSqRef);
225 for( i = 1; i <= fSplineNumber; i++ )
227 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
229 if( fVerbose > 0 )
G4cout<<i<<
"; dNdxPAI = "<<fDifPAIySection[i]<<
G4endl;
231 IntegralPAIySection();
244 static const G4double p0 = 1.20923e+00;
245 static const G4double p1 = 3.53256e-01;
246 static const G4double p2 = -1.45052e-03;
251 for( i = 0; i < numberOfElements; i++ )
254 sumZ += thisMaterialZ[i];
255 thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
257 for( i = 0; i < numberOfElements; i++ )
259 sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
261 fLowEnergyCof = sumCof;
262 delete [] thisMaterialZ;
263 delete [] thisMaterialCof;
275 G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
276 fLorentzFactor[fRefGammaNumber] - 1;
280 NormShift(betaGammaSq);
281 SplainPAI(betaGammaSq);
283 IntegralPAIySection();
287 for( i = 0; i<= fSplineNumber; i++)
289 fPAItable[i][fRefGammaNumber] = fIntegralPAIySection[i];
291 if(i != 0) fPAItable[i][0] = fSplineEnergy[i];
293 fPAItable[0][0] = fSplineNumber;
295 for(
G4int j = 1; j < 112; j++)
297 if( j == fRefGammaNumber )
continue;
299 betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
301 for(i = 1; i <= fSplineNumber; i++)
303 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
304 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
305 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
307 IntegralPAIySection();
311 for(i = 0; i <= fSplineNumber; i++)
313 fPAItable[i][j] = fIntegralPAIySection[i];
327 for( i = 1; i <= fIntervalNumber-1; i++ )
329 for( j = 1; j <= 2; j++ )
331 fSplineNumber = (i-1)*2 + j;
333 if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
334 else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
339 fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
343 for(i=2;i<=fSplineNumber;i++)
345 if(fSplineEnergy[i]<fEnergyInterval[j+1])
347 fIntegralTerm[i] = fIntegralTerm[i-1] +
348 RutherfordIntegral(j,fSplineEnergy[i-1],
353 G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
354 fEnergyInterval[j+1] );
356 fIntegralTerm[i] = fIntegralTerm[i-1] + x +
357 RutherfordIntegral(j,fEnergyInterval[j],
363 fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
370 for(
G4int k=1;k<=fIntervalNumber-1;k++)
375 fImPartDielectricConst[i] = fNormalizationCof*
376 ImPartDielectricConst(k,fSplineEnergy[i]);
377 fRePartDielectricConst[i] = fNormalizationCof*
378 RePartDielectricConst(fSplineEnergy[i]);
379 fIntegralTerm[i] *= fNormalizationCof;
381 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
382 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
383 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
400 while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
402 if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
412 for(
G4int j = fSplineNumber; j >= i+2; j-- )
414 fSplineEnergy[j] = fSplineEnergy[j-1];
415 fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
416 fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
417 fIntegralTerm[j] = fIntegralTerm[j-1];
419 fDifPAIySection[j] = fDifPAIySection[j-1];
420 fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
421 fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
429 fSplineEnergy[i+1] = en1;
441 fImPartDielectricConst[i+1] = fNormalizationCof*
442 ImPartDielectricConst(k,fSplineEnergy[i+1]);
443 fRePartDielectricConst[i+1] = fNormalizationCof*
444 RePartDielectricConst(fSplineEnergy[i+1]);
445 fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
446 RutherfordIntegral(k,fSplineEnergy[i],
449 fDifPAIySection[i+1] = DifPAIySection(i+1,betaGammaSq);
450 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
451 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
456 G4double x = 2*(fDifPAIySection[i+1] - y)/(fDifPAIySection[i+1] + y);
458 G4double delta = 2.*(fSplineEnergy[i+1]-fSplineEnergy[i])/(fSplineEnergy[i+1]+fSplineEnergy[i]);
464 if( x > fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*fDelta )
486 c1 = (x2 - x1)/x1/x2;
487 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2;
488 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
491 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
504 G4double energy2,energy3,energy4,result;
506 energy2 = energy1*energy1;
507 energy3 = energy2*energy1;
508 energy4 = energy3*energy1;
510 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
511 result *=
hbarc/energy1;
526 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
527 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
532 for(
G4int i=1;i<=fIntervalNumber-1;i++)
534 x1 = fEnergyInterval[i];
535 x2 = fEnergyInterval[i+1];
546 xln3 = log((x2 + x0)/(x1 + x0));
551 c1 = (x2 - x1)/x1/x2;
552 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2;
553 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
555 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
556 result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
557 result -= fA3[i]*c2/2/x02;
558 result -= fA4[i]*c3/3/x02;
560 cof1 = fA1[i]/x02 + fA3[i]/x04;
561 cof2 = fA2[i]/x03 + fA4[i]/x05;
563 result += 0.5*(cof1 +cof2)*xln2;
564 result += 0.5*(cof1 - cof2)*xln3;
581 G4double beta, be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
586 be2 = betaGammaSq/(1 + betaGammaSq);
592 if( betaGammaSq < 0.01 ) x2 = log(be2);
595 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
596 (1/betaGammaSq - fRePartDielectricConst[i]) +
597 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
599 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
605 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
606 x5 = -1 - fRePartDielectricConst[i] +
607 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
608 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
610 x7 = atan2(fImPartDielectricConst[i],x3);
615 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/
hbarc;
617 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
618 fImPartDielectricConst[i]*fImPartDielectricConst[i];
620 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
621 if(result < 1.0e-8) result = 1.0e-8;
627 result *= (1 - exp(-beta/betaBohr/lowCof));
648 G4double logarithm, x3, x5, argument, modul2, dNdxC;
653 be2 = betaGammaSq/(1 + betaGammaSq);
656 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq);
659 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
660 (1/betaGammaSq - fRePartDielectricConst[i]) +
661 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
662 logarithm += log(1+1.0/betaGammaSq);
665 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
671 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
672 x5 = -1.0 - fRePartDielectricConst[i] +
673 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
674 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
675 if( x3 == 0.0 ) argument = 0.5*
pi;
676 else argument = atan2(fImPartDielectricConst[i],x3);
679 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/
hbarc;
681 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
685 dNdxC *= (1-exp(-be4/betaBohr4));
689 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
690 fImPartDielectricConst[i]*fImPartDielectricConst[i];
707 G4double cof, resonance, modul2, dNdxP;
712 be2 = betaGammaSq/(1 + betaGammaSq);
716 resonance *= fImPartDielectricConst[i]/
hbarc;
718 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
720 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
723 dNdxP *= (1-exp(-be4/betaBohr4));
727 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
728 fImPartDielectricConst[i]*fImPartDielectricConst[i];
745 fIntegralPAIySection[fSplineNumber] = 0;
746 fIntegralPAIdEdx[fSplineNumber] = 0;
747 fIntegralPAIySection[0] = 0;
748 G4int k = fIntervalNumber -1;
750 for(
G4int i = fSplineNumber-1; i >= 1; i--)
752 if(fSplineEnergy[i] >= fEnergyInterval[k])
754 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + SumOverInterval(i);
755 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
759 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] +
760 SumOverBorder(i+1,fEnergyInterval[k]);
761 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
762 SumOverBorderdEdx(i+1,fEnergyInterval[k]);
777 fIntegralCerenkov[fSplineNumber] = 0;
778 fIntegralCerenkov[0] = 0;
779 k = fIntervalNumber -1;
781 for( i = fSplineNumber-1; i >= 1; i-- )
783 if(fSplineEnergy[i] >= fEnergyInterval[k])
785 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
790 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
791 SumOverBordCerenkov(i+1,fEnergyInterval[k]);
807 fIntegralPlasmon[fSplineNumber] = 0;
808 fIntegralPlasmon[0] = 0;
809 G4int k = fIntervalNumber -1;
810 for(
G4int i=fSplineNumber-1;i>=1;i--)
812 if(fSplineEnergy[i] >= fEnergyInterval[k])
814 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
818 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
819 SumOverBordPlasmon(i+1,fEnergyInterval[k]);
836 x0 = fSplineEnergy[i];
837 x1 = fSplineEnergy[i+1];
839 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)
return 0.;
841 y0 = fDifPAIySection[i];
842 yy1 = fDifPAIySection[i+1];
846 a = log10(yy1/y0)/log10(c);
853 result = b*log(x1/x0);
857 result = y0*(x1*pow(c,a-1) - x0)/a;
862 fIntegralPAIySection[0] += b*log(x1/x0);
866 fIntegralPAIySection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
878 x0 = fSplineEnergy[i];
879 x1 = fSplineEnergy[i+1];
881 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)
return 0.;
883 y0 = fDifPAIySection[i];
884 yy1 = fDifPAIySection[i+1];
886 a = log10(yy1/y0)/log10(c);
892 result = b*log(x1/x0);
896 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
912 x0 = fSplineEnergy[i];
913 x1 = fSplineEnergy[i+1];
915 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)
return 0.;
917 y0 = fdNdxCerenkov[i];
918 yy1 = fdNdxCerenkov[i+1];
923 a = log10(yy1/y0)/log10(c);
925 if(a < 20.) b = y0/pow(x0,a);
928 if(a == 0) result = b*log(c);
929 else result = y0*(x1*pow(c,a-1) - x0)/a;
932 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
933 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
949 x0 = fSplineEnergy[i];
950 x1 = fSplineEnergy[i+1];
952 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)
return 0.;
954 y0 = fdNdxPlasmon[i];
955 yy1 = fdNdxPlasmon[i+1];
957 a = log10(yy1/y0)/log10(c);
960 if(a < 20.) b = y0/pow(x0,a);
963 if(a == 0) result = b*log(x1/x0);
964 else result = y0*(x1*pow(c,a-1) - x0)/a;
967 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
968 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
985 x0 = fSplineEnergy[i];
986 x1 = fSplineEnergy[i+1];
987 y0 = fDifPAIySection[i];
988 yy1 = fDifPAIySection[i+1];
992 a = log10(yy1/y0)/log10(x1/x0);
995 if(a < 20.) b = y0/pow(x0,a);
1000 result = b*log(x0/e0);
1004 result = y0*(x0 - e0*pow(d,a-1))/a;
1009 fIntegralPAIySection[0] += b*log(x0/e0);
1013 fIntegralPAIySection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1015 x0 = fSplineEnergy[i - 1];
1016 x1 = fSplineEnergy[i - 2];
1017 y0 = fDifPAIySection[i - 1];
1018 yy1 = fDifPAIySection[i - 2];
1022 a = log10(yy1/y0)/log10(x1/x0);
1028 result += b*log(e0/x0);
1032 result += y0*(e0*pow(d,a-1) - x0)/a;
1037 fIntegralPAIySection[0] += b*log(e0/x0);
1041 fIntegralPAIySection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1055 x0 = fSplineEnergy[i];
1056 x1 = fSplineEnergy[i+1];
1057 y0 = fDifPAIySection[i];
1058 yy1 = fDifPAIySection[i+1];
1062 a = log10(yy1/y0)/log10(x1/x0);
1065 if(a < 20.) b = y0/pow(x0,a);
1070 result = b*log(x0/e0);
1074 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1076 x0 = fSplineEnergy[i - 1];
1077 x1 = fSplineEnergy[i - 2];
1078 y0 = fDifPAIySection[i - 1];
1079 yy1 = fDifPAIySection[i - 2];
1083 a = log10(yy1/y0)/log10(x1/x0);
1085 if(a < 20.) b = y0/pow(x0,a);
1090 result += b*log(e0/x0);
1094 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1111 x0 = fSplineEnergy[i];
1112 x1 = fSplineEnergy[i+1];
1113 y0 = fdNdxCerenkov[i];
1114 yy1 = fdNdxCerenkov[i+1];
1121 a = log10(yy1/y0)/log10(c);
1124 if(a < 20.) b = y0/pow(x0,a);
1127 if( a == 0 ) result = b*log(x0/e0);
1128 else result = y0*(x0 - e0*pow(d,a-1))/a;
1131 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
1132 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1136 x0 = fSplineEnergy[i - 1];
1137 x1 = fSplineEnergy[i - 2];
1138 y0 = fdNdxCerenkov[i - 1];
1139 yy1 = fdNdxCerenkov[i - 2];
1146 a = log10(yy1/y0)/log10(x1/x0);
1149 if(a < 20.) b = y0/pow(x0,a);
1151 if(a > 20.0) b = 0.0;
1152 else b = y0/pow(x0,a);
1157 if( a == 0 ) result += b*log(e0/x0);
1158 else result += y0*(e0*pow(d,a-1) - x0 )/a;
1162 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
1163 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1182 x0 = fSplineEnergy[i];
1183 x1 = fSplineEnergy[i+1];
1184 y0 = fdNdxPlasmon[i];
1185 yy1 = fdNdxPlasmon[i+1];
1189 a = log10(yy1/y0)/log10(c);
1192 if(a < 20.) b = y0/pow(x0,a);
1195 if( a == 0 ) result = b*log(x0/e0);
1196 else result = y0*(x0 - e0*pow(d,a-1))/a;
1199 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
1200 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1202 x0 = fSplineEnergy[i - 1];
1203 x1 = fSplineEnergy[i - 2];
1204 y0 = fdNdxPlasmon[i - 1];
1205 yy1 = fdNdxPlasmon[i - 2];
1209 a = log10(yy1/y0)/log10(c);
1211 if(a < 20.) b = y0/pow(x0,a);
1214 if( a == 0 ) result += b*log(e0/x0);
1215 else result += y0*(e0*pow(d,a-1) - x0)/a;
1218 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
1219 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1240 meanNumber = fIntegralPAIySection[1]*step;
1241 numOfCollisions =
G4Poisson(meanNumber);
1245 while(numOfCollisions)
1249 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1251 if( position >= fIntegralPAIySection[iTransfer] )
break;
1253 loss += fSplineEnergy[iTransfer] ;
1276 meanNumber = fIntegralCerenkov[1]*step;
1277 numOfCollisions =
G4Poisson(meanNumber);
1281 while(numOfCollisions)
1285 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1287 if( position >= fIntegralCerenkov[iTransfer] )
break;
1289 loss += fSplineEnergy[iTransfer] ;
1312 meanNumber = fIntegralPlasmon[1]*step;
1313 numOfCollisions =
G4Poisson(meanNumber);
1317 while(numOfCollisions)
1321 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1323 if( position >= fIntegralPlasmon[iTransfer] )
break;
1325 loss += fSplineEnergy[iTransfer] ;
1336 void G4PAIySection::CallError(
G4int i,
const G4String& methodName)
const
1338 G4String head =
"G4PAIySection::" + methodName +
"()";
1340 ed <<
"Wrong index " << i <<
" fSplineNumber= " << fSplineNumber;
1349 G4int G4PAIySection::fNumberOfGammas = 111;
1351 const G4double G4PAIySection::fLorentzFactor[112] =
1354 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
1355 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00,
1356 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
1357 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00,
1358 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
1359 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00,
1360 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
1361 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01,
1362 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
1363 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01,
1364 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
1365 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02,
1366 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
1367 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02,
1368 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
1369 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03,
1370 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
1371 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03,
1372 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
1373 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04,
1374 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
1375 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04,
1385 G4int G4PAIySection::fRefGammaNumber = 29;
G4long G4Poisson(G4double mean)
G4double SumOverIntervaldEdx(G4int intervalNumber)
std::ostringstream G4ExceptionDescription
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
void NormShift(G4double betaGammaSq)
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
G4double SumOverInterCerenkov(G4int intervalNumber)
G4double GetDensity() const
G4int GetMaxInterval() const
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
G4double GetSandiaMatTablePAI(G4int, G4int)
void ComputeLowEnergyCof(const G4Material *material)
const G4Element * GetElement(G4int iel) const
G4double SumOverBorderdEdx(G4int intervalNumber, G4double energy)
G4double SumOverBordPlasmon(G4int intervalNumber, G4double energy)
G4GLOB_DLL std::ostream G4cout
G4double GetStepEnergyLoss(G4double step)
G4double GetElectronDensity() const
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
G4double GetStepCerenkovLoss(G4double step)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double SumOverBorder(G4int intervalNumber, G4double energy)
void SplainPAI(G4double betaGammaSq)
G4double SumOverBordCerenkov(G4int intervalNumber, G4double energy)
size_t GetNumberOfElements() const
G4double RePartDielectricConst(G4double energy)
G4double GetStepPlasmonLoss(G4double step)
G4double DifPAIySection(G4int intervalNumber, G4double betaGammaSq)
G4double SumOverInterval(G4int intervalNumber)
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
void IntegralPAIySection()
G4double SumOverInterPlasmon(G4int intervalNumber)