00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071 #include "G4GoudsmitSaundersonMscModel.hh"
00072 #include "G4GoudsmitSaundersonTable.hh"
00073
00074 #include "G4PhysicalConstants.hh"
00075 #include "G4SystemOfUnits.hh"
00076
00077 #include "G4ParticleChangeForMSC.hh"
00078 #include "G4MaterialCutsCouple.hh"
00079 #include "G4DynamicParticle.hh"
00080 #include "G4Electron.hh"
00081 #include "G4Positron.hh"
00082
00083 #include "G4LossTableManager.hh"
00084 #include "G4Track.hh"
00085 #include "G4PhysicsTable.hh"
00086 #include "Randomize.hh"
00087
00088 using namespace std;
00089
00090 G4double G4GoudsmitSaundersonMscModel::ener[] = {-1.};
00091 G4double G4GoudsmitSaundersonMscModel::TCSE[103][106] ;
00092 G4double G4GoudsmitSaundersonMscModel::FTCSE[103][106] ;
00093 G4double G4GoudsmitSaundersonMscModel::TCSP[103][106] ;
00094 G4double G4GoudsmitSaundersonMscModel::FTCSP[103][106] ;
00095
00096
00097 G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam)
00098 : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(100.*TeV)
00099 {
00100 currentKinEnergy=currentRange=skindepth=par1=par2=par3
00101 =zPathLength=truePathLength
00102 =tausmall=taulim=tlimit=charge=lambdalimit=tPathLength=lambda0=lambda1
00103 =lambda11=mass=0.0;
00104 currentMaterialIndex = -1;
00105
00106 fr=0.02,rangeinit=0.,masslimite=0.6*MeV,
00107 particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm;
00108 tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig;
00109 tlimitminfix=1.e-6*mm;stepmin=tlimitminfix;lambdalimit=1.*mm;smallstep=1.e10;
00110 theManager=G4LossTableManager::Instance();
00111 inside=false;insideskin=false;
00112 samplez=false;
00113 firstStep = true;
00114
00115 GSTable = new G4GoudsmitSaundersonTable();
00116
00117 if(ener[0] < 0.0){
00118 G4cout << "### G4GoudsmitSaundersonMscModel loading ELSEPA data" << G4endl;
00119 LoadELSEPAXSections();
00120 }
00121 }
00122
00123 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel()
00124 {
00125 delete GSTable;
00126 }
00127
00128 void G4GoudsmitSaundersonMscModel::Initialise(const G4ParticleDefinition* p,
00129 const G4DataVector&)
00130 {
00131 skindepth=skin*stepmin;
00132 SetParticle(p);
00133 fParticleChange = GetParticleChangeForMSC(p);
00134 }
00135
00136
00137
00138 G4double
00139 G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
00140 G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
00141 {
00142 G4double kinEnergy = kineticEnergy;
00143 if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
00144 if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
00145
00146 G4double cs(0.0), cs0(0.0);
00147 CalculateIntegrals(p,Z,kinEnergy,cs0,cs);
00148
00149 return cs;
00150 }
00151
00152
00153 G4ThreeVector&
00154 G4GoudsmitSaundersonMscModel::SampleScattering(const G4ThreeVector& oldDirection, G4double)
00155 {
00156 fDisplacement.set(0.0,0.0,0.0);
00157 G4double kineticEnergy = currentKinEnergy;
00158
00159 if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)||
00160 (tPathLength/tausmall < lambda1)) { return fDisplacement; }
00161
00163
00164 G4double eloss = 0.0;
00165 if (tPathLength > currentRange*dtrl) {
00166 eloss = kineticEnergy -
00167 GetEnergy(particle,currentRange-tPathLength,currentCouple);
00168 } else {
00169 eloss = tPathLength*GetDEDX(particle,kineticEnergy,currentCouple);
00170 }
00171
00172
00173
00174
00175
00176
00177
00178 kineticEnergy -= 0.5*eloss;
00179
00181
00182 const G4Material* mat = currentCouple->GetMaterial();
00183 const G4ElementVector* theElementVector = mat->GetElementVector();
00184 const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
00185 G4int nelm = mat->GetNumberOfElements();
00186 G4double s0(0.0), s1(0.0);
00187 lambda0 = 0.0;
00188 for(G4int i=0;i<nelm;i++)
00189 {
00190 CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,s0,s1);
00191 lambda0 += (theAtomNumDensityVector[i]*s0);
00192 }
00193 if(lambda0>0.0) lambda0 =1./lambda0;
00194
00195
00196
00197 G4double g1=0.0;
00198 if(lambda1>0.0) { g1 = lambda0/lambda1; }
00199
00200 G4double logx0,x1,delta;
00201 G4double x0=g1*0.5;
00202
00203 for(G4int i=0;i<1000;++i)
00204 {
00205 logx0=std::log(1.+1./x0);
00206 x1 = x0-(x0*((1.+x0)*logx0-1.0)-g1*0.5)/( (1.+2.*x0)*logx0-2.0);
00207
00208
00209 if(x1 < 0.0) { x1 = 0.5*x0; }
00210 else if(x1 > 2*x0) { x1 = 2*x0; }
00211 else if(x1 < 0.5*x0) { x1 = 0.5*x0; }
00212 delta = std::fabs( x1 - x0 );
00213 x0 = x1;
00214 if(delta < 1.0e-3*x1) { break;}
00215 }
00216 G4double scrA = x1;
00217
00218 G4double lambdan=0.;
00219
00220 if(lambda0>0.0) { lambdan=tPathLength/lambda0; }
00221 if(lambdan<=1.0e-12) { return fDisplacement; }
00222
00223
00224
00225
00226 G4double Qn1 = lambdan *g1;
00227 G4double Qn12 = 0.5*Qn1;
00228
00229 G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
00230 G4double cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
00231 G4double us=0.0,vs=0.0,ws=1.0,wss=0.,x_coord=0.0,y_coord=0.0,z_coord=1.0;
00232
00233 G4double epsilon1=G4UniformRand();
00234 G4double expn = std::exp(-lambdan);
00235
00236 if(epsilon1<expn)
00237 { return fDisplacement; }
00238 else if((epsilon1<((1.+lambdan)*expn))||(lambdan<1.))
00239 {
00240 G4double xi=G4UniformRand();
00241 xi= 2.*scrA*xi/(1.-xi + scrA);
00242 if(xi<0.)xi=0.;
00243 else if(xi>2.)xi=2.;
00244 ws=(1. - xi);
00245 wss=std::sqrt(xi*(2.-xi));
00246 G4double phi0=CLHEP::twopi*G4UniformRand();
00247 us=wss*cos(phi0);
00248 vs=wss*sin(phi0);
00249 }
00250 else
00251 {
00252
00253
00254 SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
00255 G4double phi1 = CLHEP::twopi*G4UniformRand();
00256 cosPhi1 = cos(phi1);
00257 sinPhi1 = sin(phi1);
00258
00259
00260 SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
00261 G4double phi2 = CLHEP::twopi*G4UniformRand();
00262 cosPhi2 = cos(phi2);
00263 sinPhi2 = sin(phi2);
00264
00265
00266 us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
00267 vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
00268 ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2;
00269 G4double sqrtA=sqrt(scrA);
00270 if(acos(ws)<sqrtA)
00271 {
00272 G4int i=0;
00273 do{i++;
00274 ws=1.+Qn12*log(G4UniformRand());
00275 }while((fabs(ws)>1.)&&(i<20));
00276 if(i>=19)ws=cos(sqrtA);
00277 wss=std::sqrt((1.-ws*ws));
00278 us=wss*std::cos(phi1);
00279 vs=wss*std::sin(phi1);
00280 }
00281 }
00282
00283
00284 G4ThreeVector newDirection(us,vs,ws);
00285 newDirection.rotateUz(oldDirection);
00286 fParticleChange->ProposeMomentumDirection(newDirection);
00287
00288
00289 if(Qn1<0.02) { z_coord = 1.0 - Qn1*(0.5 - Qn1/6.); }
00290 else { z_coord = (1.-std::exp(-Qn1))/Qn1; }
00291 G4double rr = zPathLength*std::sqrt((1.- z_coord*z_coord)/(1.-ws*ws));
00292 x_coord = rr*us;
00293 y_coord = rr*vs;
00294
00295
00296 z_coord -= 1.0;
00297 z_coord *= zPathLength;
00298
00299
00300
00301
00302
00303
00304
00305
00306 fDisplacement.set(x_coord,y_coord,z_coord);
00307 fDisplacement.rotateUz(oldDirection);
00308
00309 return fDisplacement;
00310 }
00311
00312
00313
00314 void
00315 G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA,
00316 G4double &cost, G4double &sint)
00317 {
00318 G4double r1,tet,xi=0.;
00319 G4double Qn1 = 2.* lambdan;
00320 if(scrA < 10.) { Qn1 *= scrA*((1.+scrA)*log(1.+1./scrA)-1.); }
00321 else { Qn1*= (1.0 - 0.5/scrA - 0.5/(scrA*scrA)) ; }
00322 if (Qn1<0.001)
00323 {
00324 do{
00325 r1=G4UniformRand();
00326 xi=-0.5*Qn1*log(G4UniformRand());
00327 tet=acos(1.-xi);
00328 }while(tet*r1*r1>sin(tet));
00329 }
00330 else if(Qn1>0.5) { xi=2.*G4UniformRand(); }
00331 else{ xi=2.*(GSTable->SampleTheta(lambdan,scrA,G4UniformRand()));}
00332
00333
00334 if(xi<0.)xi=0.;
00335 else if(xi>2.)xi=2.;
00336
00337 cost=(1. - xi);
00338 sint=sqrt(xi*(2.-xi));
00339
00340 }
00341
00342
00343
00344
00345
00346 void
00347 G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z,
00348 G4double kinEnergy,G4double &Sig0,
00349 G4double &Sig1)
00350 {
00351 G4double x1,x2,y1,y2,acoeff,bcoeff;
00352 G4double kineticE = kinEnergy;
00353 if(kineticE<lowKEnergy)kineticE=lowKEnergy;
00354 if(kineticE>highKEnergy)kineticE=highKEnergy;
00355 kineticE /= eV;
00356 G4double logE=std::log(kineticE);
00357
00358 G4int iZ = G4int(Z);
00359 if(iZ > 103) iZ = 103;
00360
00361 G4int enerInd=0;
00362 for(G4int i=0;i<105;i++)
00363 {
00364 if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;break;}
00365 }
00366
00367 if(p==G4Electron::Electron())
00368 {
00369 if(kineticE<=1.0e+9)
00370 {
00371 x1=ener[enerInd];
00372 x2=ener[enerInd+1];
00373 y1=TCSE[iZ-1][enerInd];
00374 y2=TCSE[iZ-1][enerInd+1];
00375 acoeff=(y2-y1)/(x2*x2-x1*x1);
00376 bcoeff=y2-acoeff*x2*x2;
00377 Sig0=acoeff*logE*logE+bcoeff;
00378 Sig0 =std::exp(Sig0);
00379 y1=FTCSE[iZ-1][enerInd];
00380 y2=FTCSE[iZ-1][enerInd+1];
00381 acoeff=(y2-y1)/(x2*x2-x1*x1);
00382 bcoeff=y2-acoeff*x2*x2;
00383 Sig1=acoeff*logE*logE+bcoeff;
00384 Sig1=std::exp(Sig1);
00385 }
00386 else
00387 {
00388 x1=ener[104];
00389 x2=ener[105];
00390 y1=TCSE[iZ-1][104];
00391 y2=TCSE[iZ-1][105];
00392 Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
00393 Sig0=std::exp(Sig0);
00394 y1=FTCSE[iZ-1][104];
00395 y2=FTCSE[iZ-1][105];
00396 Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
00397 Sig1=std::exp(Sig1);
00398 }
00399 }
00400 if(p==G4Positron::Positron())
00401 {
00402 if(kinEnergy<=1.0e+9)
00403 {
00404 x1=ener[enerInd];
00405 x2=ener[enerInd+1];
00406 y1=TCSP[iZ-1][enerInd];
00407 y2=TCSP[iZ-1][enerInd+1];
00408 acoeff=(y2-y1)/(x2*x2-x1*x1);
00409 bcoeff=y2-acoeff*x2*x2;
00410 Sig0=acoeff*logE*logE+bcoeff;
00411 Sig0 =std::exp(Sig0);
00412 y1=FTCSP[iZ-1][enerInd];
00413 y2=FTCSP[iZ-1][enerInd+1];
00414 acoeff=(y2-y1)/(x2*x2-x1*x1);
00415 bcoeff=y2-acoeff*x2*x2;
00416 Sig1=acoeff*logE*logE+bcoeff;
00417 Sig1=std::exp(Sig1);
00418 }
00419 else
00420 {
00421 x1=ener[104];
00422 x2=ener[105];
00423 y1=TCSP[iZ-1][104];
00424 y2=TCSP[iZ-1][105];
00425 Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
00426 Sig0 =std::exp(Sig0);
00427 y1=FTCSP[iZ-1][104];
00428 y2=FTCSP[iZ-1][105];
00429 Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
00430 Sig1=std::exp(Sig1);
00431 }
00432 }
00433
00434 Sig0 *= barn;
00435 Sig1 *= barn;
00436
00437 }
00438
00439
00440
00441 void G4GoudsmitSaundersonMscModel::StartTracking(G4Track* track)
00442 {
00443 SetParticle(track->GetDynamicParticle()->GetDefinition());
00444 firstStep = true;
00445 inside = false;
00446 insideskin = false;
00447 tlimit = geombig;
00448 }
00449
00450
00451
00452
00453 G4double
00454 G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track,
00455 G4double& currentMinimalStep)
00456 {
00457 tPathLength = currentMinimalStep;
00458 const G4DynamicParticle* dp = track.GetDynamicParticle();
00459 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
00460 G4StepStatus stepStatus = sp->GetStepStatus();
00461 currentCouple = track.GetMaterialCutsCouple();
00462 SetCurrentCouple(currentCouple);
00463 currentMaterialIndex = currentCouple->GetIndex();
00464 currentKinEnergy = dp->GetKineticEnergy();
00465 currentRange = GetRange(particle,currentKinEnergy,currentCouple);
00466
00467 lambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
00468
00469
00470 if(inside || tPathLength < tlimitminfix) {
00471 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00472 }
00473 if(tPathLength > currentRange) tPathLength = currentRange;
00474
00475 G4double presafety = sp->GetSafety();
00476
00477
00478
00479
00480
00481
00482
00483 if(currentRange < presafety)
00484 {
00485 inside = true;
00486 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00487 }
00488
00489
00490
00491 if (steppingAlgorithm == fUseDistanceToBoundary)
00492 {
00493
00494 G4double geomlimit = ComputeGeomLimit(track, presafety, tPathLength);
00495
00496
00497 if(currentRange <= presafety)
00498 {
00499 inside = true;
00500 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00501 }
00502
00503 smallstep += 1.;
00504 insideskin = false;
00505
00506 if(firstStep || stepStatus == fGeomBoundary)
00507 {
00508 rangeinit = currentRange;
00509 if(firstStep) smallstep = 1.e10;
00510 else smallstep = 1.;
00511
00512
00513
00514 G4double rat = currentKinEnergy/MeV ;
00515 rat = 1.e-3/(rat*(10.+rat)) ;
00516
00517 stepmin = rat*lambda1;
00518 skindepth = skin*stepmin;
00519
00520 tlimitmin = 10.*stepmin;
00521 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
00522
00523
00524
00525
00526 if((geomlimit < geombig) && (geomlimit > geommin))
00527 {
00528 if(stepStatus == fGeomBoundary)
00529 tgeom = geomlimit/facgeom;
00530 else
00531 tgeom = 2.*geomlimit/facgeom;
00532 }
00533 else
00534 tgeom = geombig;
00535
00536 }
00537
00538
00539 tlimit = facrange*rangeinit;
00540 if(tlimit < facsafety*presafety)
00541 tlimit = facsafety*presafety;
00542
00543
00544 if(tlimit < tlimitmin) tlimit = tlimitmin;
00545
00546 if(tlimit > tgeom) tlimit = tgeom;
00547
00548
00549
00550
00551
00552 if((tPathLength < tlimit) && (tPathLength < presafety) &&
00553 (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
00554 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00555
00556
00557 if(smallstep < skin)
00558 {
00559 tlimit = stepmin;
00560 insideskin = true;
00561 }
00562 else if(geomlimit < geombig)
00563 {
00564 if(geomlimit > skindepth)
00565 {
00566 if(tlimit > geomlimit-0.999*skindepth)
00567 tlimit = geomlimit-0.999*skindepth;
00568 }
00569 else
00570 {
00571 insideskin = true;
00572 if(tlimit > stepmin) tlimit = stepmin;
00573 }
00574 }
00575
00576 if(tlimit < stepmin) tlimit = stepmin;
00577
00578 if(tPathLength > tlimit) tPathLength = tlimit;
00579
00580 }
00581
00582
00583 else if(steppingAlgorithm == fUseSafety)
00584 {
00585
00586
00587 if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
00588 presafety = ComputeSafety(sp->GetPosition(),tPathLength);
00589
00590
00591 if(currentRange < presafety)
00592 {
00593 inside = true;
00594 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00595 }
00596
00597 if(firstStep || stepStatus == fGeomBoundary)
00598 {
00599 rangeinit = currentRange;
00600 fr = facrange;
00601
00602 if(mass < masslimite)
00603 {
00604 if(lambda1 > currentRange)
00605 rangeinit = lambda1;
00606 if(lambda1 > lambdalimit)
00607 fr *= 0.75+0.25*lambda1/lambdalimit;
00608 }
00609
00610
00611 G4double rat = currentKinEnergy/MeV ;
00612 rat = 1.e-3/(rat*(10.+rat)) ;
00613 tlimitmin = 10.*lambda1*rat;
00614 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
00615 }
00616
00617 tlimit = fr*rangeinit;
00618
00619 if(tlimit < facsafety*presafety)
00620 tlimit = facsafety*presafety;
00621
00622
00623 if(tlimit < tlimitmin) tlimit = tlimitmin;
00624
00625 if(tPathLength > tlimit) tPathLength = tlimit;
00626 }
00627
00628
00629 else
00630 {
00631 if (stepStatus == fGeomBoundary)
00632 {
00633 if (currentRange > lambda1) tlimit = facrange*currentRange;
00634 else tlimit = facrange*lambda1;
00635
00636 if(tlimit < tlimitmin) tlimit = tlimitmin;
00637 if(tPathLength > tlimit) tPathLength = tlimit;
00638 }
00639 }
00640
00641
00642 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00643 }
00644
00645
00646
00647 G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double)
00648 {
00649 firstStep = false;
00650 par1 = -1. ;
00651 par2 = par3 = 0. ;
00652
00653
00654 zPathLength = tPathLength;
00655
00656
00657 if(tPathLength < tlimitminfix) { return zPathLength; }
00658
00659
00660
00661 if(tPathLength > currentRange)
00662 { tPathLength = currentRange; }
00663
00664 G4double tau = tPathLength/lambda1 ;
00665
00666 if ((tau <= tausmall) || insideskin) {
00667 zPathLength = tPathLength;
00668 if(zPathLength > lambda1) { zPathLength = lambda1; }
00669 return zPathLength;
00670 }
00671
00672 G4double zmean = tPathLength;
00673 if (tPathLength < currentRange*dtrl) {
00674 if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
00675 else zmean = lambda1*(1.-exp(-tau));
00676 } else if(currentKinEnergy < mass || tPathLength == currentRange) {
00677 par1 = 1./currentRange ;
00678 par2 = 1./(par1*lambda1) ;
00679 par3 = 1.+par2 ;
00680 if(tPathLength < currentRange)
00681 zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
00682 else
00683 zmean = 1./(par1*par3) ;
00684 } else {
00685 G4double T1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
00686
00687 lambda11 = GetTransportMeanFreePath(particle,T1);
00688
00689 par1 = (lambda1-lambda11)/(lambda1*tPathLength) ;
00690 par2 = 1./(par1*lambda1) ;
00691 par3 = 1.+par2 ;
00692 zmean = (1.-exp(par3*log(lambda11/lambda1)))/(par1*par3) ;
00693 }
00694
00695 zPathLength = zmean ;
00696
00697 if(samplez) {
00698
00699 const G4double ztmax = 0.99;
00700 G4double zt = zmean/tPathLength ;
00701
00702 if (tPathLength > stepmin && zt < ztmax) {
00703
00704 G4double u,cz1;
00705 if(zt >= 0.333333333) {
00706
00707 G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
00708 cz1 = 1.+cz ;
00709 G4double u0 = cz/cz1 ;
00710 G4double grej ;
00711 do {
00712 u = exp(log(G4UniformRand())/cz1) ;
00713 grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
00714 } while (grej < G4UniformRand()) ;
00715
00716 } else {
00717 cz1 = 1./zt-1.;
00718 u = 1.-exp(log(G4UniformRand())/cz1) ;
00719 }
00720 zPathLength = tPathLength*u ;
00721 }
00722 }
00723 if(zPathLength > lambda1) zPathLength = lambda1;
00724
00725
00726 return zPathLength;
00727 }
00728
00729
00730
00731 G4double
00732 G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength)
00733 {
00734
00735 if(geomStepLength == zPathLength && tPathLength <= currentRange)
00736 return tPathLength;
00737
00738
00739 zPathLength = geomStepLength;
00740 tPathLength = geomStepLength;
00741 if(geomStepLength < tlimitminfix) return tPathLength;
00742
00743
00744 if((geomStepLength > lambda1*tausmall) && !insideskin)
00745 {
00746 if(par1 < 0.)
00747 tPathLength = -lambda1*log(1.-geomStepLength/lambda1) ;
00748 else
00749 {
00750 if(par1*par3*geomStepLength < 1.)
00751 tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
00752 else
00753 tPathLength = currentRange;
00754 }
00755 }
00756 if(tPathLength < geomStepLength) tPathLength = geomStepLength;
00757
00758
00759 return tPathLength;
00760 }
00761
00762
00763
00764
00765 void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()
00766 {
00767 G4String filename = "XSECTIONS.dat";
00768
00769 char* path = getenv("G4LEDATA");
00770 if (!path)
00771 {
00772 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()","em0006",
00773 FatalException,
00774 "Environment variable G4LEDATA not defined");
00775 return;
00776 }
00777
00778 G4String pathString(path);
00779 G4String dirFile = pathString + "/msc_GS/" + filename;
00780 FILE *infile;
00781 infile = fopen(dirFile,"r");
00782 if (infile == 0)
00783 {
00784 G4ExceptionDescription ed;
00785 ed << "Data file <" + dirFile + "> is not opened!" << G4endl;
00786 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00787 "em0003",FatalException,ed);
00788 return;
00789 }
00790
00791
00792 G4float aRead;
00793 for(G4int i=0 ; i<106 ;i++){
00794 if(1 == fscanf(infile,"%f\t",&aRead)) {
00795 if(aRead > 0.0) { aRead = log(aRead); }
00796 else { aRead = 0.0; }
00797 } else {
00798 G4ExceptionDescription ed;
00799 ed << "Error reading <" + dirFile + "> loop #1 i= " << i << G4endl;
00800 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00801 "em0003",FatalException,ed);
00802 return;
00803 }
00804 ener[i]=aRead;
00805 }
00806 for(G4int j=0;j<103;j++){
00807 for(G4int i=0;i<106;i++){
00808 if(1 == fscanf(infile,"%f\t",&aRead)) {
00809 if(aRead > 0.0) { aRead = log(aRead); }
00810 else { aRead = 0.0; }
00811 } else {
00812 G4ExceptionDescription ed;
00813 ed << "Error reading <" + dirFile + "> loop #2 j= " << j
00814 << "; i= " << i << G4endl;
00815 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00816 "em0003",FatalException,ed);
00817 return;
00818 }
00819 TCSE[j][i]=aRead;
00820 }
00821 }
00822 for(G4int j=0;j<103;j++){
00823 for(G4int i=0;i<106;i++){
00824 if(1 == fscanf(infile,"%f\t",&aRead)) {
00825 if(aRead > 0.0) { aRead = log(aRead); }
00826 else { aRead = 0.0; }
00827 } else {
00828 G4ExceptionDescription ed;
00829 ed << "Error reading <" + dirFile + "> loop #3 j= " << j
00830 << "; i= " << i << G4endl;
00831 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00832 "em0003",FatalException,ed);
00833 return;
00834 }
00835 FTCSE[j][i]=aRead;
00836 }
00837 }
00838 for(G4int j=0;j<103;j++){
00839 for(G4int i=0;i<106;i++){
00840 if(1 == fscanf(infile,"%f\t",&aRead)) {
00841 if(aRead > 0.0) { aRead = log(aRead); }
00842 else { aRead = 0.0; }
00843 } else {
00844 G4ExceptionDescription ed;
00845 ed << "Error reading <" + dirFile + "> loop #4 j= " << j
00846 << "; i= " << i << G4endl;
00847 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00848 "em0003",FatalException,ed);
00849 return;
00850 }
00851 TCSP[j][i]=aRead;
00852 }
00853 }
00854 for(G4int j=0;j<103;j++){
00855 for(G4int i=0;i<106;i++){
00856 if(1 == fscanf(infile,"%f\t",&aRead)) {
00857 if(aRead > 0.0) { aRead = log(aRead); }
00858 else { aRead = 0.0; }
00859 } else {
00860 G4ExceptionDescription ed;
00861 ed << "Error reading <" + dirFile + "> loop #5 j= " << j
00862 << "; i= " << i << G4endl;
00863 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00864 "em0003",FatalException,ed);
00865 return;
00866 }
00867 FTCSP[j][i]=aRead;
00868 }
00869 }
00870
00871 fclose(infile);
00872
00873 }
00874
00875