480 currentRange =
GetRange(particle,currentKinEnergy,currentCouple,
488 fTheTrueStepLenght = currentMinimalStep;
489 fTheTransportDistance = currentMinimalStep;
490 fTheZPathLenght = currentMinimalStep;
491 fTheDisplacementVector.
set(0.,0.,0.);
492 fTheNewDirection.
set(0.,0.,1.);
495 fIsEverythingWasDone =
false;
497 fIsMultipleSacettring =
false;
499 fIsSingleScattering =
false;
501 fIsNoScatteringInMSC =
false;
504 fIsNoDisplace =
false;
506 presafety = sp->GetSafety();
511 distance *= (1.20-fZeff*(1.62e-2-9.22e-5*fZeff));
519 if (gIsOptimizationOn && (distance<presafety)) {
521 fIsMultipleSacettring =
true;
522 fIsNoDisplace =
true;
538 fIsWasOnBoundary =
true;
541 skindepth =
skin*fLambda0;
544 fIsInsideSkin =
false;
550 if ((stepStatus==
fGeomBoundary) || (presafety<skindepth) || (fTheTrueStepLenght<skindepth)) {
552 if ((stepStatus ==
fGeomBoundary) || (presafety < skindepth)) {
553 fIsInsideSkin =
true;
554 fIsWasOnBoundary =
true;
567 if (sslimit<fTheTrueStepLenght) {
568 fTheTrueStepLenght = sslimit;
569 fIsSingleScattering =
true;
572 fTheZPathLenght = fTheTrueStepLenght;
576 fIsEverythingWasDone =
true;
584 fIsMultipleSacettring =
true;
587 fIsFirstRealStep =
false;
591 if (fIsWasOnBoundary && !fIsInsideSkin) {
593 fIsWasOnBoundary =
false;
594 fIsFirstRealStep =
true;
603 if (firstStep || fIsFirstRealStep || rangeinit>1.e+20) {
604 rangeinit = currentRange;
615 if (geomlimit<geombig) {
619 if ((1.-geomlimit/fLambda1)> 0.) {
620 geomlimit = -fLambda1*
G4Log(1.-geomlimit/fLambda1);
638 tlimit = std::min(tlimit,tgeom);
648 if (geomlimit<geombig) {
649 tlimit = std::min(tlimit, geomlimit-0.999*skindepth);
652 if (firstStep || fIsFirstRealStep) {
653 fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit());
655 fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit);
659 presafety =
ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
660 geomlimit = presafety;
662 skindepth =
skin*fLambda0;
668 if ((stepStatus==
fGeomBoundary) || (presafety<skindepth) || (fTheTrueStepLenght<skindepth)) {
680 if (sslimit<fTheTrueStepLenght) {
681 fTheTrueStepLenght = sslimit;
682 fIsSingleScattering =
true;
685 fTheZPathLenght = fTheTrueStepLenght;
689 fIsEverythingWasDone =
true;
694 fIsMultipleSacettring =
true;
695 fIsEverythingWasDone =
true;
697 fTheTrueStepLenght = std::min(fTheTrueStepLenght,
facrange*currentRange);
700 if (fTheTrueStepLenght>presafety) {
701 fTheTrueStepLenght = std::min(fTheTrueStepLenght, presafety);
707 fTheTrueStepLenght = std::min(fTheTrueStepLenght, fLambda1*0.5);
719 fIsMultipleSacettring =
true;
721 presafety =
ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
724 if ((distance<presafety) && (gIsOptimizationOn)) {
725 fIsNoDisplace =
true;
728 if (firstStep || (stepStatus==
fGeomBoundary) || rangeinit>1.e+20) {
729 rangeinit = currentRange;
740 tlimit = std::max(fr*rangeinit,
facsafety*presafety);
743 fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit());
745 fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit);
753 if (fIsEverythingWasDone) {
754 if (fIsSingleScattering) {
758 G4double pt2 = currentKinEnergy*(currentKinEnergy+2.0*CLHEP::electron_mass_c2);
759 G4double beta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
762 if (cost<-1.) cost = -1.;
763 if (cost> 1.) cost = 1.;
766 G4double sint = std::sqrt(dum*(2.-dum));
770 fTheNewDirection.
set(sint*cosPhi,sint*sinPhi,cost);
771 }
else if (fIsMultipleSacettring) {
790 if (!fIsEverythingWasDone) {
793 fTheTrueStepLenght = std::min(fTheTrueStepLenght,currentRange);
795 fTheZPathLenght = fTheTrueStepLenght;
797 if (fTheTrueStepLenght<tlimitminfix2) {
798 return fTheZPathLenght;
800 G4double tau = fTheTrueStepLenght/fLambda1;
802 fTheZPathLenght = std::min(fTheTrueStepLenght, fLambda1);
803 }
else if (fTheTrueStepLenght<currentRange*
dtrl) {
804 if (tau<taulim) fTheZPathLenght = fTheTrueStepLenght*(1.-0.5*tau) ;
805 else fTheZPathLenght = fLambda1*(1.-
G4Exp(-tau));
806 }
else if (currentKinEnergy<mass || fTheTrueStepLenght==currentRange) {
807 par1 = 1./currentRange ;
808 par2 = 1./(par1*fLambda1) ;
810 if (fTheTrueStepLenght<currentRange) {
811 fTheZPathLenght = 1./(par1*par3) * (1.-std::pow(1.-par1*fTheTrueStepLenght,par3));
813 fTheZPathLenght = 1./(par1*par3);
816 G4double rfin = std::max(currentRange-fTheTrueStepLenght, 0.01*currentRange);
818 G4double lambda1 = GetTransportMeanFreePathOnly(particle,T1);
820 par1 = (fLambda1-lambda1)/(fLambda1*fTheTrueStepLenght);
821 par2 = 1./(par1*fLambda1);
824 fTheZPathLenght = 1./(par1*par3) * (1.-g4calc->
powA(1.-par1*fTheTrueStepLenght,par3));
827 fTheZPathLenght = std::min(fTheZPathLenght, fLambda1);
829 return fTheZPathLenght;
924 fIsNoScatteringInMSC =
false;
926 G4double kineticEnergy = currentKinEnergy;
931 eloss = kineticEnergy -
GetEnergy(particle,currentRange-fTheTrueStepLenght,currentCouple);
943 G4double efStep = fTheTrueStepLenght;
945 G4double kineticEnergy0 = kineticEnergy;
946 if (gIsUseAccurate) {
947 kineticEnergy -= 0.5*eloss;
949 tau = kineticEnergy/electron_mass_c2;
951 eps0 = eloss/kineticEnergy0;
952 epsm = eloss/kineticEnergy;
954 efEnergy = kineticEnergy * (1.-epsm*epsm*(6.+10.*tau+5.*tau2)/(24.*tau2+48.*tau+72.));
955 G4double dum = 0.166666*(4.+tau*(6.+tau*(7.+tau*(4.+tau))))*(epsm/((tau+1.)*(tau+2.)))*(epsm/((tau+1.)*(tau+2.)));
956 efStep = fTheTrueStepLenght*(1.-dum);
958 kineticEnergy -= 0.5*eloss;
959 efEnergy = kineticEnergy;
960 G4double factor = 1./(1.+0.9784671*kineticEnergy);
961 eps0 = eloss/kineticEnergy0;
962 epsm = eps0/(1.-0.5*eps0);
963 G4double temp = 0.3*(1 -factor*(1.-0.333333*factor))*eps0*eps0;
964 efStep = fTheTrueStepLenght*(1.+temp);
973 lambdan=efStep/fLambda0;
975 if (lambdan<=1.0e-12) {
976 if (fIsEverythingWasDone) {
977 fTheZPathLenght = fTheTrueStepLenght;
979 fIsNoScatteringInMSC =
true;
986 G4double cosTheta1 = 1.0, sinTheta1 = 0.0, cosTheta2 = 1.0, sinTheta2 = 0.0;
987 G4double cosPhi1 = 1.0, sinPhi1 = 0.0, cosPhi2 = 1.0, sinPhi2 = 0.0;
988 G4double uss = 0.0, vss = 0.0, wss = 1.0;
989 G4double x_coord = 0.0, y_coord = 0.0, z_coord = 1.0;
995 sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+cosTheta1));
997 sinTheta2 = std::sqrt((1.-cosTheta2)*(1.+cosTheta2));
1001 G4double pt2 = efEnergy*(efEnergy+2.0*CLHEP::electron_mass_c2);
1002 G4double beta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
1006 G4int mcEkinIdx = -1;
1007 G4int mcDeltIdx = -1;
1009 G4bool isMsc = fGSTable->
Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta1, sinTheta1, lekin, beta2,
1010 currentMaterialIndex, &gsDtr, mcEkinIdx, mcDeltIdx, transfPar,
1012 fGSTable->
Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta2, sinTheta2, lekin, beta2,
1013 currentMaterialIndex, &gsDtr, mcEkinIdx, mcDeltIdx, transfPar, !isMsc);
1014 if (cosTheta1+cosTheta2==2.) {
1015 if (fIsEverythingWasDone)
1016 fTheZPathLenght = fTheTrueStepLenght;
1017 fIsNoScatteringInMSC =
true;
1023 sinPhi1 = std::sin(phi1);
1024 cosPhi1 = std::cos(phi1);
1026 sinPhi2 = std::sin(phi2);
1027 cosPhi2 = std::cos(phi2);
1030 u2 = sinTheta2*cosPhi2;
1031 v2 = sinTheta2*sinPhi2;
1032 G4double u2p = cosTheta1*u2 + sinTheta1*cosTheta2;
1033 uss = u2p*cosPhi1 - v2*sinPhi1;
1034 vss = u2p*sinPhi1 + v2*cosPhi1;
1035 wss = cosTheta1*cosTheta2 - sinTheta1*u2;
1038 fTheNewDirection.
set(uss,vss,wss);
1042 if(fIsNoDisplace && fIsEverythingWasDone)
1043 fTheZPathLenght = fTheTrueStepLenght;
1052 if (gIsUseAccurate) {
1055 if(Qn1<0.7) par = 1.;
1056 else if (Qn1<7.0) par = -0.031376*Qn1+1.01356;
1063 G4double gamma = 6.0*fScrA*(1.0 + fScrA)*(loga*(1.0 + 2.0*fScrA) - 2.0)/fG1;
1064 gamma *= fMCtoG2PerG1;
1072 G4double delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1;
1076 G4double temp = (2.0+tau*temp1)/((tau+1.0)*temp1);
1078 temp = temp - (tau+1.0)/((tau+2.0)*(loga*(1.0+fScrA)-1.0));
1081 delta = delta + 0.40824829*(eps0*(tau+1.0)/((tau+2.0)*
1082 (loga*(1.0+fScrA)-1.0)*(loga*(1.0+2.0*fScrA)-2.0)) - 0.25*temp*temp);
1089 G4double ut = b*sinTheta1*cosPhi1 + c*(cosPhi1*u2 - sinPhi1*w1v2) + eta1*uss*temp1;
1090 G4double vt = b*sinTheta1*sinPhi1 + c*(sinPhi1*u2 + cosPhi1*w1v2) + eta1*vss*temp1;
1091 G4double wt = eta1*(1+temp) + b*cosTheta1 + c*cosTheta2 + eta1*wss*temp1;
1100 x_coord = ut*fTheTrueStepLenght;
1101 y_coord = vt*fTheTrueStepLenght;
1102 z_coord = wt*fTheTrueStepLenght;
1104 if(fIsEverythingWasDone){
1108 G4double transportDistance = std::sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
1110 if(transportDistance>fTheTrueStepLenght)
1111 transportDistance = fTheTrueStepLenght;
1112 fTheZPathLenght = transportDistance;
1116 fTheDisplacementVector.
set(x_coord,y_coord,z_coord-fTheZPathLenght);
1123 if(fIsEverythingWasDone){
1127 zz = 1.0 - Qn1*(0.5 - Qn1*(0.166666667 - 0.041666667*Qn1));
1129 zz = (1.-
G4Exp(-Qn1))/Qn1;
1134 zz = fTheZPathLenght/fTheTrueStepLenght;
1137 G4double rr = (1.-zz*zz)/(1.-wss*wss);
1138 if(rr >= 0.25) rr = 0.25;
1139 G4double rperp = fTheTrueStepLenght*std::sqrt(rr);
1140 x_coord = rperp*uss;
1141 y_coord = rperp*vss;
1142 z_coord = zz*fTheTrueStepLenght;
1144 if(fIsEverythingWasDone){
1145 G4double transportDistance = std::sqrt(x_coord*x_coord + y_coord*y_coord + z_coord*z_coord);
1146 fTheZPathLenght = transportDistance;
1149 fTheDisplacementVector.
set(x_coord,y_coord,z_coord- fTheZPathLenght);