93 G4double G4GoudsmitSaundersonMscModel::ener[] = {-1.};
94 G4double G4GoudsmitSaundersonMscModel::TCSE[103][106] ;
95 G4double G4GoudsmitSaundersonMscModel::FTCSE[103][106] ;
96 G4double G4GoudsmitSaundersonMscModel::TCSP[103][106] ;
97 G4double G4GoudsmitSaundersonMscModel::FTCSP[103][106] ;
103 currentKinEnergy=currentRange=skindepth=par1=par2=par3
104 =zPathLength=truePathLength
105 =tausmall=taulim=tlimit=charge=lambdalimit=tPathLength=lambda0=lambda1
107 currentMaterialIndex = -1;
109 fr=0.02,rangeinit=0.,masslimite=0.6*
MeV,
110 particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*
mm;
111 tlimitmin=10.e-6*
mm;geombig=1.e50*
mm;geommin=1.e-3*
mm,tgeom=geombig;
112 tlimitminfix=1.e-6*
mm;stepmin=tlimitminfix;lambdalimit=1.*
mm;smallstep=1.e10;
114 inside=
false;insideskin=
false;
121 G4cout <<
"### G4GoudsmitSaundersonMscModel loading ELSEPA data" <<
G4endl;
122 LoadELSEPAXSections();
134 skindepth=
skin*stepmin;
146 if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
147 if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
150 CalculateIntegrals(p,Z,kinEnergy,cs0,cs);
160 G4double kineticEnergy = currentKinEnergy;
162 if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)||
168 if (tPathLength > currentRange*
dtrl) {
169 eloss = kineticEnergy -
170 GetEnergy(particle,currentRange-tPathLength,currentCouple);
172 eloss = tPathLength*
GetDEDX(particle,kineticEnergy,currentCouple);
181 kineticEnergy -= 0.5*eloss;
191 for(
G4int i=0;i<nelm;i++)
193 CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,
s0,s1);
194 lambda0 += (theAtomNumDensityVector[i]*
s0);
196 if(lambda0>0.0) { lambda0 =1./lambda0; }
201 if(lambda1>0.0) { g1 = lambda0/lambda1; }
206 for(
G4int i=0;i<1000;++i)
208 logx0 =
G4Log(1.+1./x0);
209 x1 = x0-(x0*((1.+x0)*logx0-1.0)-g1*0.5)/( (1.+2.*x0)*logx0-2.0);
212 if(x1 < 0.0) { x1 = 0.5*x0; }
213 else if(x1 > 2*x0) { x1 = 2*x0; }
214 else if(x1 < 0.5*x0) { x1 = 0.5*x0; }
215 delta = std::fabs( x1 - x0 );
217 if(delta < 1.0e-3*x1) {
break;}
223 if(lambda0>0.0) { lambdan=tPathLength/lambda0; }
232 G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
233 G4double cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
234 G4double us=0.0,vs=0.0,ws=1.0,wss=0.,x_coord=0.0,y_coord=0.0,z_coord=1.0;
241 else if((epsilon1<((1.+lambdan)*expn))||(lambdan<1.))
244 xi= 2.*scrA*xi/(1.-xi + scrA);
248 wss=std::sqrt(xi*(2.-xi));
257 SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
263 SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
269 us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
270 vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
271 ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2;
278 }
while((fabs(ws)>1.)&&(i<20));
279 if(i>=19)ws=cos(sqrtA);
280 wss=std::sqrt((1.-ws*ws));
281 us=wss*std::cos(phi1);
282 vs=wss*std::sin(phi1);
288 newDirection.
rotateUz(oldDirection);
292 if(Qn1<0.02) { z_coord = 1.0 - Qn1*(0.5 - Qn1/6.); }
293 else { z_coord = (1.-
G4Exp(-Qn1))/Qn1; }
294 G4double rr = zPathLength*std::sqrt((1.- z_coord*z_coord)/(1.-ws*ws));
300 z_coord *= zPathLength;
318 G4GoudsmitSaundersonMscModel::SampleCosineTheta(
G4double lambdan,
G4double scrA,
323 if(scrA < 10.) { Qn1 *= scrA*((1.+scrA)*
G4Log(1.+1./scrA)-1.); }
324 else { Qn1*= (1.0 - 0.5/scrA - 0.5/(scrA*scrA)) ; }
331 }
while(tet*r1*r1>sin(tet));
341 sint=sqrt(xi*(2.-xi));
356 if(kineticE<lowKEnergy)kineticE=lowKEnergy;
357 if(kineticE>highKEnergy)kineticE=highKEnergy;
362 if(iZ > 103) iZ = 103;
365 for(
G4int i=0;i<105;i++)
367 if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;
break;}
376 y1=TCSE[iZ-1][enerInd];
377 y2=TCSE[iZ-1][enerInd+1];
378 acoeff=(y2-y1)/(x2*x2-x1*x1);
379 bcoeff=y2-acoeff*x2*x2;
380 Sig0=acoeff*logE*logE+bcoeff;
382 y1=FTCSE[iZ-1][enerInd];
383 y2=FTCSE[iZ-1][enerInd+1];
384 acoeff=(y2-y1)/(x2*x2-x1*x1);
385 bcoeff=y2-acoeff*x2*x2;
386 Sig1=acoeff*logE*logE+bcoeff;
395 Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
399 Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
405 if(kinEnergy<=1.0e+9)
409 y1=TCSP[iZ-1][enerInd];
410 y2=TCSP[iZ-1][enerInd+1];
411 acoeff=(y2-y1)/(x2*x2-x1*x1);
412 bcoeff=y2-acoeff*x2*x2;
413 Sig0=acoeff*logE*logE+bcoeff;
415 y1=FTCSP[iZ-1][enerInd];
416 y2=FTCSP[iZ-1][enerInd+1];
417 acoeff=(y2-y1)/(x2*x2-x1*x1);
418 bcoeff=y2-acoeff*x2*x2;
419 Sig1=acoeff*logE*logE+bcoeff;
428 Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
432 Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
460 tPathLength = currentMinimalStep;
466 currentMaterialIndex = currentCouple->
GetIndex();
468 currentRange =
GetRange(particle,currentKinEnergy,currentCouple);
473 if(inside || tPathLength < tlimitminfix) {
476 if(tPathLength > currentRange) tPathLength = currentRange;
486 if(currentRange < presafety)
500 if(currentRange <= presafety)
511 rangeinit = currentRange;
512 if(firstStep) smallstep = 1.e10;
518 rat = 1.e-3/(rat*(10.+rat)) ;
520 stepmin = rat*lambda1;
521 skindepth =
skin*stepmin;
523 tlimitmin = 10.*stepmin;
524 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
529 if((geomlimit < geombig) && (geomlimit > geommin))
547 if(tlimit < tlimitmin) tlimit = tlimitmin;
549 if(tlimit > tgeom) tlimit = tgeom;
555 if((tPathLength < tlimit) && (tPathLength < presafety) &&
556 (smallstep >=
skin) && (tPathLength < geomlimit-0.999*skindepth))
565 else if(geomlimit < geombig)
567 if(geomlimit > skindepth)
569 if(tlimit > geomlimit-0.999*skindepth)
570 tlimit = geomlimit-0.999*skindepth;
575 if(tlimit > stepmin) tlimit = stepmin;
579 if(tlimit < stepmin) tlimit = stepmin;
581 if(tPathLength > tlimit) tPathLength = tlimit;
590 if((stepStatus !=
fGeomBoundary) && (presafety < tlimitminfix))
594 if(currentRange < presafety)
602 rangeinit = currentRange;
605 if(mass < masslimite)
607 if(lambda1 > currentRange)
609 if(lambda1 > lambdalimit)
610 fr *= 0.75+0.25*lambda1/lambdalimit;
615 rat = 1.e-3/(rat*(10.+rat)) ;
616 tlimitmin = 10.*lambda1*rat;
617 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
620 tlimit = fr*rangeinit;
626 if(tlimit < tlimitmin) tlimit = tlimitmin;
628 if(tPathLength > tlimit) tPathLength = tlimit;
636 if (currentRange > lambda1) tlimit =
facrange*currentRange;
639 if(tlimit < tlimitmin) tlimit = tlimitmin;
640 if(tPathLength > tlimit) tPathLength = tlimit;
657 zPathLength = tPathLength;
660 if(tPathLength < tlimitminfix) {
return zPathLength; }
664 if(tPathLength > currentRange)
665 { tPathLength = currentRange; }
667 G4double tau = tPathLength/lambda1 ;
669 if ((tau <= tausmall) || insideskin) {
670 zPathLength = tPathLength;
671 if(zPathLength > lambda1) { zPathLength = lambda1; }
676 if (tPathLength < currentRange*
dtrl) {
677 if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
678 else zmean = lambda1*(1.-
G4Exp(-tau));
679 }
else if(currentKinEnergy < mass || tPathLength == currentRange) {
680 par1 = 1./currentRange ;
681 par2 = 1./(par1*lambda1) ;
683 if(tPathLength < currentRange)
684 zmean = (1.-
G4Exp(par3*
G4Log(1.-tPathLength/currentRange)))/(par1*par3) ;
686 zmean = 1./(par1*par3) ;
692 par1 = (lambda1-lambda11)/(lambda1*tPathLength) ;
693 par2 = 1./(par1*lambda1) ;
695 zmean = (1.-
G4Exp(par3*
G4Log(lambda11/lambda1)))/(par1*par3) ;
698 zPathLength = zmean ;
705 if (tPathLength > stepmin && zt < ztmax) {
708 if(zt >= 0.333333333) {
710 G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
716 grej = exp(cz*
G4Log(u/u0))*(1.-u)/(1.-u0) ;
723 zPathLength = tPathLength*u ;
726 if(zPathLength > lambda1) zPathLength = lambda1;
738 if(geomStepLength == zPathLength && tPathLength <= currentRange)
742 zPathLength = geomStepLength;
743 tPathLength = geomStepLength;
744 if(geomStepLength < tlimitminfix)
return tPathLength;
747 if((geomStepLength > lambda1*tausmall) && !insideskin)
750 tPathLength = -lambda1*
G4Log(1.-geomStepLength/lambda1) ;
753 if(par1*par3*geomStepLength < 1.)
754 tPathLength = (1.-
G4Exp(
G4Log(1.-par1*par3*geomStepLength)/par3))/par1 ;
756 tPathLength = currentRange;
759 if(tPathLength < geomStepLength) tPathLength = geomStepLength;
768 void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()
770 G4String filename =
"XSECTIONS.dat";
772 char* path = getenv(
"G4LEDATA");
775 G4Exception(
"G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
"em0006",
777 "Environment variable G4LEDATA not defined");
782 G4String dirFile = pathString +
"/msc_GS/" + filename;
783 std::ifstream infile(dirFile, std::ios::in);
784 if( !infile.is_open()) {
786 ed <<
"Data file <" + dirFile +
"> is not opened!";
787 G4Exception(
"G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
794 for(
G4int i=0 ; i<106 ;i++){
798 ed <<
"Error reading <" + dirFile +
"> loop #1 i= " << i;
799 G4Exception(
"G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
803 if(aRead > 0.0) { aRead =
G4Log(aRead); }
804 else { aRead = 0.0; }
808 for(
G4int j=0;j<103;j++){
809 for(
G4int i=0;i<106;i++){
813 ed <<
"Error reading <" + dirFile +
"> loop #2 j= " << j
815 G4Exception(
"G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
819 if(aRead > 0.0) { aRead =
G4Log(aRead); }
820 else { aRead = 0.0; }
825 for(
G4int j=0;j<103;j++){
826 for(
G4int i=0;i<106;i++){
830 ed <<
"Error reading <" + dirFile +
"> loop #3 j= " << j
832 G4Exception(
"G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
836 if(aRead > 0.0) { aRead =
G4Log(aRead); }
837 else { aRead = 0.0; }
842 for(
G4int j=0;j<103;j++){
843 for(
G4int i=0;i<106;i++){
847 ed <<
"Error reading <" + dirFile +
"> loop #4 j= " << j
849 G4Exception(
"G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
853 if(aRead > 0.0) { aRead =
G4Log(aRead); }
854 else { aRead = 0.0; }
859 for(
G4int j=0;j<103;j++){
860 for(
G4int i=0;i<106;i++){
864 ed <<
"Error reading <" + dirFile +
"> loop #5 j= " << j
866 G4Exception(
"G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
870 if(aRead > 0.0) { aRead =
G4Log(aRead); }
871 else { aRead = 0.0; }
void set(double x, double y, double z)
G4GoudsmitSaundersonMscModel(const G4String &nam="GoudsmitSaunderson")
virtual G4double ComputeTrueStepLength(G4double geomStepLength)
static G4LossTableManager * Instance()
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
const G4DynamicParticle * GetDynamicParticle() const
G4StepStatus GetStepStatus() const
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4double SampleTheta(G4double, G4double, G4double)
G4ParticleDefinition * GetDefinition() const
const G4Step * GetStep() const
const G4ElementVector * GetElementVector() const
G4StepPoint * GetPreStepPoint() const
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
void StartTracking(G4Track *)
virtual G4double ComputeTruePathLengthLimit(const G4Track &track, G4double ¤tMinimalStep)
const G4double * GetVecNbOfAtomsPerVolume() const
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetPosition() const
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Hep3Vector & rotateUz(const Hep3Vector &)
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
G4ThreeVector fDisplacement
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double ComputeSafety(const G4ThreeVector &position, G4double limit)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
static G4Positron * Positron()
G4double GetDEDX(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
void SetCurrentCouple(const G4MaterialCutsCouple *)
virtual ~G4GoudsmitSaundersonMscModel()
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double, G4double, G4double)
G4double GetSafety() const
G4MscStepLimitType steppingAlgorithm
static G4Electron * Electron()
size_t GetNumberOfElements() const
virtual G4ThreeVector & SampleScattering(const G4ThreeVector &, G4double safety)
const G4Material * GetMaterial() const
virtual G4double ComputeGeomPathLength(G4double truePathLength)