923 {
924 fIsNoScatteringInMSC = false;
925
926 G4double kineticEnergy = currentKinEnergy;
927
928
930
931 eloss = kineticEnergy -
GetEnergy(particle,currentRange-fTheTrueStepLenght,currentCouple);
932
933
934
935
940
941
943 G4double efStep = fTheTrueStepLenght;
944
945 G4double kineticEnergy0 = kineticEnergy;
946 if (gIsUseAccurate) {
947 kineticEnergy -= 0.5*eloss;
948
949 tau = kineticEnergy/electron_mass_c2;
950 tau2 = tau*tau;
951 eps0 = eloss/kineticEnergy0;
952 epsm = eloss/kineticEnergy;
953
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);
957 } else {
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);
965 }
966
967
968
970
972 if (fLambda0>0.0) {
973 lambdan=efStep/fLambda0;
974 }
975 if (lambdan<=1.0e-12) {
976 if (fIsEverythingWasDone) {
977 fTheZPathLenght = fTheTrueStepLenght;
978 }
979 fIsNoScatteringInMSC = true;
980 return;
981 }
982
984
985
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;
991
992
993 if (0.5*Qn1 > 7.0){
995 sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+cosTheta1));
997 sinTheta2 = std::sqrt((1.-cosTheta2)*(1.+cosTheta2));
998 } else {
999
1001 G4double pt2 = efEnergy*(efEnergy+2.0*CLHEP::electron_mass_c2);
1002 G4double beta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
1003
1004
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,
1011 true);
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;
1018 return;
1019 }
1020 }
1021
1023 sinPhi1 = std::sin(phi1);
1024 cosPhi1 = std::cos(phi1);
1026 sinPhi2 = std::sin(phi2);
1027 cosPhi2 = std::cos(phi2);
1028
1029
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;
1036
1037
1038 fTheNewDirection.
set(uss,vss,wss);
1039
1040
1041
1042 if(fIsNoDisplace && fIsEverythingWasDone)
1043 fTheZPathLenght = fTheTrueStepLenght;
1044
1045
1046 if(fIsNoDisplace)
1047 return;
1048
1049
1050
1051 Qn1 *= fMCtoQ1;
1052 if (gIsUseAccurate) {
1053
1055 if(Qn1<0.7) par = 1.;
1056 else if (Qn1<7.0) par = -0.031376*Qn1+1.01356;
1057 else par = 0.79;
1058
1059
1060
1061
1063 G4double gamma = 6.0*fScrA*(1.0 + fScrA)*(loga*(1.0 + 2.0*fScrA) - 2.0)/fG1;
1064 gamma *= fMCtoG2PerG1;
1065
1068
1069
1070
1071
1072 G4double delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1;
1073
1074
1076 G4double temp = (2.0+tau*temp1)/((tau+1.0)*temp1);
1077
1078 temp = temp - (tau+1.0)/((tau+2.0)*(loga*(1.0+fScrA)-1.0));
1079 temp = temp * epsm;
1080 temp1 = 1.0 - temp;
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);
1085
1086
1087
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;
1092
1093
1094 ut *=par;
1095 vt *=par;
1096 wt *=par;
1097
1098
1099
1100 x_coord = ut*fTheTrueStepLenght;
1101 y_coord = vt*fTheTrueStepLenght;
1102 z_coord = wt*fTheTrueStepLenght;
1103
1104 if(fIsEverythingWasDone){
1105
1106
1107
1108 G4double transportDistance = std::sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
1109
1110 if(transportDistance>fTheTrueStepLenght)
1111 transportDistance = fTheTrueStepLenght;
1112 fTheZPathLenght = transportDistance;
1113 }
1114
1115
1116 fTheDisplacementVector.
set(x_coord,y_coord,z_coord-fTheZPathLenght);
1117 } else {
1118
1119
1120
1121
1123 if(fIsEverythingWasDone){
1124
1125
1126 if(Qn1<0.1) {
1127 zz = 1.0 - Qn1*(0.5 - Qn1*(0.166666667 - 0.041666667*Qn1));
1128 } else {
1129 zz = (1.-
G4Exp(-Qn1))/Qn1;
1130 }
1131 } else {
1132
1133
1134 zz = fTheZPathLenght/fTheTrueStepLenght;
1135 }
1136
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;
1143
1144 if(fIsEverythingWasDone){
1145 G4double transportDistance = std::sqrt(x_coord*x_coord + y_coord*y_coord + z_coord*z_coord);
1146 fTheZPathLenght = transportDistance;
1147 }
1148
1149 fTheDisplacementVector.
set(x_coord,y_coord,z_coord- fTheZPathLenght);
1150 }
1151}
G4bool Sampling(G4double lambdaval, G4double qval, G4double scra, G4double &cost, G4double &sint, G4double lekin, G4double beta2, G4int matindx, GSMSCAngularDtr **gsDtr, G4int &mcekini, G4int &mcdelti, G4double &transfPar, G4bool isfirst)