Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4GoudsmitSaundersonMscModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id$
27//
28// -------------------------------------------------------------------
29//
30// GEANT4 Class file
31//
32// File name: G4GoudsmitSaundersonMscModel
33//
34// Author: Omrane Kadri
35//
36// Creation date: 20.02.2009
37//
38// Modifications:
39// 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
40//
41// 15.04.2009 O.Kadri: cleanup: discard no scattering and single scattering theta
42// sampling from SampleCosineTheta() which means the splitting
43// step into two sub-steps occur only for msc regime
44//
45// 12.06.2009 O.Kadri: linear log-log extrapolation of lambda0 & lambda1 between 1 GeV - 100 TeV
46// adding a theta min limit due to screening effect of the atomic nucleus
47// 26.08.2009 O.Kadri: Cubic Spline interpolation was replaced with polynomial method
48// within CalculateIntegrals method
49// 05.10.2009 O.Kadri: tuning small angle theta distributions
50// assuming the case of lambdan<1 as single scattering regime
51// tuning theta sampling for theta below the screening angle
52// 08.02.2010 O.Kadri: bugfix in compound xsection calculation and small angle computation
53// adding a rejection condition to hard collision angular sampling
54// ComputeTruePathLengthLimit was taken from G4WentzelVIModel
55// 26.03.2010 O.Kadri: direct xsection calculation not inverse of the inverse
56// angular sampling without large angle rejection method
57// longitudinal displacement is computed exactly from <z>
58// 12.05.2010 O.Kadri: exchange between target and projectile has as a condition the particle type (e-/e-)
59// some cleanup to minimize time consuming (adding lamdan12 & Qn12, changing the error to 1.0e-12 for scrA)
60//
61//
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63//REFERENCES:
64//Ref.1:E. Benedito et al.,"Mixed simulation ... cross-sections", NIMB 174 (2001) pp 91-110;
65//Ref.2:I. Kawrakow et al.,"On the condensed ... transport",NIMB 142 (1998) pp 253-280;
66//Ref.3:I. Kawrakow et al.,"On the representation ... calculations",NIMB 134 (1998) pp 325-336;
67//Ref.4:Bielajew et al.,".....", NIMB 173 (2001) 332-343;
68//Ref.5:F. Salvat et al.,"ELSEPA--Dirac partial ...molecules", Comp.Phys.Comm.165 (2005) pp 157-190;
69//Ref.6:G4UrbanMscModel G4 9.2;
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71#include "G4GoudsmitSaundersonMscModel.hh"
72#include "G4GoudsmitSaundersonTable.hh"
73
74#include "G4PhysicalConstants.hh"
75#include "G4SystemOfUnits.hh"
76
77#include "G4ParticleChangeForMSC.hh"
78#include "G4MaterialCutsCouple.hh"
79#include "G4DynamicParticle.hh"
80#include "G4Electron.hh"
81#include "G4Positron.hh"
82
83#include "G4LossTableManager.hh"
84#include "G4Track.hh"
85#include "G4PhysicsTable.hh"
86#include "Randomize.hh"
87
88using namespace std;
89
90G4double G4GoudsmitSaundersonMscModel::ener[] = {-1.};
91G4double G4GoudsmitSaundersonMscModel::TCSE[103][106] ;
92G4double G4GoudsmitSaundersonMscModel::FTCSE[103][106] ;
93G4double G4GoudsmitSaundersonMscModel::TCSP[103][106] ;
94G4double G4GoudsmitSaundersonMscModel::FTCSP[103][106] ;
95
96//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam)
98 : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(100.*TeV)
99{
100 currentKinEnergy=currentRange=skindepth=par1=par2=par3
101 =zPathLength=truePathLength
102 =tausmall=taulim=tlimit=charge=lambdalimit=tPathLength=lambda0=lambda1
103 =lambda11=mass=0.0;
104 currentMaterialIndex = -1;
105
106 fr=0.02,rangeinit=0.,masslimite=0.6*MeV,
107 particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm;
108 tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig;
109 tlimitminfix=1.e-6*mm;stepmin=tlimitminfix;lambdalimit=1.*mm;smallstep=1.e10;
110 theManager=G4LossTableManager::Instance();
111 inside=false;insideskin=false;
112 samplez=false;
113 firstStep = true;
114
115 GSTable = new G4GoudsmitSaundersonTable();
116
117 if(ener[0] < 0.0){
118 G4cout << "### G4GoudsmitSaundersonMscModel loading ELSEPA data" << G4endl;
119 LoadELSEPAXSections();
120 }
121}
122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel()
124{
125 delete GSTable;
126}
127//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128void G4GoudsmitSaundersonMscModel::Initialise(const G4ParticleDefinition* p,
129 const G4DataVector&)
130{
131 skindepth=skin*stepmin;
132 SetParticle(p);
133 fParticleChange = GetParticleChangeForMSC(p);
134}
135
136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137
138G4double
139G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
140 G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
141{
142 G4double kinEnergy = kineticEnergy;
143 if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
144 if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
145
146 G4double cs(0.0), cs0(0.0);
147 CalculateIntegrals(p,Z,kinEnergy,cs0,cs);
148
149 return cs;
150}
151//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152
153G4ThreeVector&
154G4GoudsmitSaundersonMscModel::SampleScattering(const G4DynamicParticle* dynParticle, G4double)
155{
156 fDisplacement.set(0.0,0.0,0.0);
157 G4double kineticEnergy = dynParticle->GetKineticEnergy();
158 if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)||
159 (tPathLength/tausmall < lambda1)) { return fDisplacement; }
160
161 ///////////////////////////////////////////
162 // Effective energy
163 G4double eloss = 0.0;
164 if (tPathLength > currentRange*dtrl) {
165 eloss = kineticEnergy -
166 GetEnergy(particle,currentRange-tPathLength,currentCouple);
167 } else {
168 eloss = tPathLength*GetDEDX(particle,kineticEnergy,currentCouple);
169 }
170 /*
171 G4double ttau = kineticEnergy/electron_mass_c2;
172 G4double ttau2 = ttau*ttau;
173 G4double epsilonpp = eloss/kineticEnergy;
174 G4double cst1 = epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72);
175 kineticEnergy *= (1 - cst1);
176 */
177 kineticEnergy -= 0.5*eloss;
178
179 ///////////////////////////////////////////
180 // additivity rule for mixture and compound xsection's
181 const G4Material* mat = currentCouple->GetMaterial();
182 const G4ElementVector* theElementVector = mat->GetElementVector();
183 const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
184 G4int nelm = mat->GetNumberOfElements();
185 G4double s0(0.0), s1(0.0);
186 lambda0 = 0.0;
187 for(G4int i=0;i<nelm;i++)
188 {
189 CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,s0,s1);
190 lambda0 += (theAtomNumDensityVector[i]*s0);
191 }
192 if(lambda0>0.0) lambda0 =1./lambda0;
193
194 // Newton-Raphson root's finding method of scrA from:
195 // Sig1(PWA)/Sig0(PWA)=g1=2*scrA*((1+scrA)*log(1+1/scrA)-1)
196 G4double g1=0.0;
197 if(lambda1>0.0) { g1 = lambda0/lambda1; }
198
199 G4double logx0,x1,delta;
200 G4double x0=g1*0.5;
201 // V.Ivanchenko added limit of the loop
202 for(G4int i=0;i<1000;++i)
203 {
204 logx0=std::log(1.+1./x0);
205 x1 = x0-(x0*((1.+x0)*logx0-1.0)-g1*0.5)/( (1.+2.*x0)*logx0-2.0);
206
207 // V.Ivanchenko cut step size of iterative procedure
208 if(x1 < 0.0) { x1 = 0.5*x0; }
209 else if(x1 > 2*x0) { x1 = 2*x0; }
210 else if(x1 < 0.5*x0) { x1 = 0.5*x0; }
211 delta = std::fabs( x1 - x0 );
212 x0 = x1;
213 if(delta < 1.0e-3*x1) { break;}
214 }
215 G4double scrA = x1;
216
217 G4double lambdan=0.;
218
219 if(lambda0>0.0) { lambdan=tPathLength/lambda0; }
220 if(lambdan<=1.0e-12) { return fDisplacement; }
221
222 //G4cout << "E(eV)= " << kineticEnergy/eV << " L0= " << lambda0
223 // << " L1= " << lambda1 << G4endl;
224
225 G4double Qn1 = lambdan *g1;//2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
226 G4double Qn12 = 0.5*Qn1;
227
228 G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
229 G4double cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
230 G4double us=0.0,vs=0.0,ws=1.0,wss=0.,x_coord=0.0,y_coord=0.0,z_coord=1.0;
231
232 G4double epsilon1=G4UniformRand();
233 G4double expn = std::exp(-lambdan);
234
235 if(epsilon1<expn)// no scattering
236 { return fDisplacement; }
237 else if((epsilon1<((1.+lambdan)*expn))||(lambdan<1.))//single or plural scattering (Rutherford DCS's)
238 {
239 G4double xi=G4UniformRand();
240 xi= 2.*scrA*xi/(1.-xi + scrA);
241 if(xi<0.)xi=0.;
242 else if(xi>2.)xi=2.;
243 ws=(1. - xi);
244 wss=std::sqrt(xi*(2.-xi));
245 G4double phi0=CLHEP::twopi*G4UniformRand();
246 us=wss*cos(phi0);
247 vs=wss*sin(phi0);
248 }
249 else // multiple scattering
250 {
251 // Ref.2 subsection 4.4 "The best solution found"
252 // Sample first substep scattering angle
253 SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
254 G4double phi1 = CLHEP::twopi*G4UniformRand();
255 cosPhi1 = cos(phi1);
256 sinPhi1 = sin(phi1);
257
258 // Sample second substep scattering angle
259 SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
260 G4double phi2 = CLHEP::twopi*G4UniformRand();
261 cosPhi2 = cos(phi2);
262 sinPhi2 = sin(phi2);
263
264 // Overall scattering direction
265 us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
266 vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
267 ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2;
268 G4double sqrtA=sqrt(scrA);
269 if(acos(ws)<sqrtA)//small angle approximation for theta less than screening angle
270 {
271 G4int i=0;
272 do{i++;
273 ws=1.+Qn12*log(G4UniformRand());
274 }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run
275 if(i>=19)ws=cos(sqrtA);
276 wss=std::sqrt((1.-ws*ws));
277 us=wss*std::cos(phi1);
278 vs=wss*std::sin(phi1);
279 }
280 }
281
282 G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
283 G4ThreeVector newDirection(us,vs,ws);
284 newDirection.rotateUz(oldDirection);
285 fParticleChange->ProposeMomentumDirection(newDirection);
286
287 // corresponding to error less than 1% in the exact formula of <z>
288 if(Qn1<0.02) { z_coord = 1.0 - Qn1*(0.5 - Qn1/6.); }
289 else { z_coord = (1.-std::exp(-Qn1))/Qn1; }
290 G4double rr = zPathLength*std::sqrt((1.- z_coord*z_coord)/(1.-ws*ws));
291 x_coord = rr*us;
292 y_coord = rr*vs;
293
294 // displacement is computed relatively to the end point
295 z_coord -= 1.0;
296 z_coord *= zPathLength;
297 /*
298 G4cout << "G4GS::SampleSecondaries: e(MeV)= " << kineticEnergy
299 << " sinTheta= " << sqrt(1.0 - ws*ws)
300 << " trueStep(mm)= " << tPathLength
301 << " geomStep(mm)= " << zPathLength
302 << G4endl;
303 */
304
305 fDisplacement.set(x_coord,y_coord,z_coord);
306 fDisplacement.rotateUz(oldDirection);
307
308 return fDisplacement;
309}
310
311//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312
313void
314G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA,
315 G4double &cost, G4double &sint)
316{
317 G4double r1,tet,xi=0.;
318 G4double Qn1 = 2.* lambdan;
319 if(scrA < 10.) { Qn1 *= scrA*((1.+scrA)*log(1.+1./scrA)-1.); }
320 else { Qn1*= (1.0 - 0.5/scrA - 0.5/(scrA*scrA)) ; }
321 if (Qn1<0.001)
322 {
323 do{
324 r1=G4UniformRand();
325 xi=-0.5*Qn1*log(G4UniformRand());
326 tet=acos(1.-xi);
327 }while(tet*r1*r1>sin(tet));
328 }
329 else if(Qn1>0.5) { xi=2.*G4UniformRand(); }//isotropic distribution
330 else{ xi=2.*(GSTable->SampleTheta(lambdan,scrA,G4UniformRand()));}
331
332
333 if(xi<0.)xi=0.;
334 else if(xi>2.)xi=2.;
335
336 cost=(1. - xi);
337 sint=sqrt(xi*(2.-xi));
338
339}
340
341//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
342// Polynomial log-log interpolation of Lambda0 and Lambda1 between 100 eV - 1 GeV
343// linear log-log extrapolation between 1 GeV - 100 TeV
344
345void
346G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z,
347 G4double kinEnergy,G4double &Sig0,
348 G4double &Sig1)
349{
350 G4double x1,x2,y1,y2,acoeff,bcoeff;
351 G4double kineticE = kinEnergy;
352 if(kineticE<lowKEnergy)kineticE=lowKEnergy;
353 if(kineticE>highKEnergy)kineticE=highKEnergy;
354 kineticE /= eV;
355 G4double logE=std::log(kineticE);
356
357 G4int iZ = G4int(Z);
358 if(iZ > 103) iZ = 103;
359
360 G4int enerInd=0;
361 for(G4int i=0;i<105;i++)
362 {
363 if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;break;}
364 }
365
366 if(p==G4Electron::Electron())
367 {
368 if(kineticE<=1.0e+9)//Interpolation of the form y=ax�+b
369 {
370 x1=ener[enerInd];
371 x2=ener[enerInd+1];
372 y1=TCSE[iZ-1][enerInd];
373 y2=TCSE[iZ-1][enerInd+1];
374 acoeff=(y2-y1)/(x2*x2-x1*x1);
375 bcoeff=y2-acoeff*x2*x2;
376 Sig0=acoeff*logE*logE+bcoeff;
377 Sig0 =std::exp(Sig0);
378 y1=FTCSE[iZ-1][enerInd];
379 y2=FTCSE[iZ-1][enerInd+1];
380 acoeff=(y2-y1)/(x2*x2-x1*x1);
381 bcoeff=y2-acoeff*x2*x2;
382 Sig1=acoeff*logE*logE+bcoeff;
383 Sig1=std::exp(Sig1);
384 }
385 else //Interpolation of the form y=ax+b
386 {
387 x1=ener[104];
388 x2=ener[105];
389 y1=TCSE[iZ-1][104];
390 y2=TCSE[iZ-1][105];
391 Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
392 Sig0=std::exp(Sig0);
393 y1=FTCSE[iZ-1][104];
394 y2=FTCSE[iZ-1][105];
395 Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
396 Sig1=std::exp(Sig1);
397 }
398 }
399 if(p==G4Positron::Positron())
400 {
401 if(kinEnergy<=1.0e+9)
402 {
403 x1=ener[enerInd];
404 x2=ener[enerInd+1];
405 y1=TCSP[iZ-1][enerInd];
406 y2=TCSP[iZ-1][enerInd+1];
407 acoeff=(y2-y1)/(x2*x2-x1*x1);
408 bcoeff=y2-acoeff*x2*x2;
409 Sig0=acoeff*logE*logE+bcoeff;
410 Sig0 =std::exp(Sig0);
411 y1=FTCSP[iZ-1][enerInd];
412 y2=FTCSP[iZ-1][enerInd+1];
413 acoeff=(y2-y1)/(x2*x2-x1*x1);
414 bcoeff=y2-acoeff*x2*x2;
415 Sig1=acoeff*logE*logE+bcoeff;
416 Sig1=std::exp(Sig1);
417 }
418 else
419 {
420 x1=ener[104];
421 x2=ener[105];
422 y1=TCSP[iZ-1][104];
423 y2=TCSP[iZ-1][105];
424 Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
425 Sig0 =std::exp(Sig0);
426 y1=FTCSP[iZ-1][104];
427 y2=FTCSP[iZ-1][105];
428 Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
429 Sig1=std::exp(Sig1);
430 }
431 }
432
433 Sig0 *= barn;
434 Sig1 *= barn;
435
436}
437
438//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
439
440void G4GoudsmitSaundersonMscModel::StartTracking(G4Track* track)
441{
442 SetParticle(track->GetDynamicParticle()->GetDefinition());
443 firstStep = true;
444 inside = false;
445 insideskin = false;
446 tlimit = geombig;
447}
448
449//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
450//t->g->t step transformations taken from Ref.6
451
452G4double
453G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track,
454 G4double& currentMinimalStep)
455{
456 tPathLength = currentMinimalStep;
457 const G4DynamicParticle* dp = track.GetDynamicParticle();
458 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
459 G4StepStatus stepStatus = sp->GetStepStatus();
460 currentCouple = track.GetMaterialCutsCouple();
461 SetCurrentCouple(currentCouple);
462 currentMaterialIndex = currentCouple->GetIndex();
463 currentKinEnergy = dp->GetKineticEnergy();
464 currentRange = GetRange(particle,currentKinEnergy,currentCouple);
465
466 lambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
467
468 // stop here if small range particle
469 if(inside || tPathLength < tlimitminfix) {
470 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
471 }
472 if(tPathLength > currentRange) tPathLength = currentRange;
473
474 G4double presafety = sp->GetSafety();
475
476 //G4cout << "G4GS::StepLimit tPathLength= "
477 // <<tPathLength<<" safety= " << presafety
478 // << " range= " <<currentRange<< " lambda= "<<lambda1
479 // << " Alg: " << steppingAlgorithm <<G4endl;
480
481 // far from geometry boundary
482 if(currentRange < presafety)
483 {
484 inside = true;
485 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
486 }
487
488 // standard version
489 //
490 if (steppingAlgorithm == fUseDistanceToBoundary)
491 {
492 //compute geomlimit and presafety
493 G4double geomlimit = ComputeGeomLimit(track, presafety, tPathLength);
494
495 // is far from boundary
496 if(currentRange <= presafety)
497 {
498 inside = true;
499 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
500 }
501
502 smallstep += 1.;
503 insideskin = false;
504
505 if(firstStep || stepStatus == fGeomBoundary)
506 {
507 rangeinit = currentRange;
508 if(firstStep) smallstep = 1.e10;
509 else smallstep = 1.;
510
511 //define stepmin here (it depends on lambda!)
512 //rough estimation of lambda_elastic/lambda_transport
513 G4double rat = currentKinEnergy/MeV ;
514 rat = 1.e-3/(rat*(10.+rat)) ;
515 //stepmin ~ lambda_elastic
516 stepmin = rat*lambda1;
517 skindepth = skin*stepmin;
518 //define tlimitmin
519 tlimitmin = 10.*stepmin;
520 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
521
522 //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
523 // << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
524 // constraint from the geometry
525 if((geomlimit < geombig) && (geomlimit > geommin))
526 {
527 if(stepStatus == fGeomBoundary)
528 tgeom = geomlimit/facgeom;
529 else
530 tgeom = 2.*geomlimit/facgeom;
531 }
532 else
533 tgeom = geombig;
534
535 }
536
537 //step limit
538 tlimit = facrange*rangeinit;
539 if(tlimit < facsafety*presafety)
540 tlimit = facsafety*presafety;
541
542 //lower limit for tlimit
543 if(tlimit < tlimitmin) tlimit = tlimitmin;
544
545 if(tlimit > tgeom) tlimit = tgeom;
546
547 //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
548 // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
549
550 // shortcut
551 if((tPathLength < tlimit) && (tPathLength < presafety) &&
552 (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
553 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
554
555 // step reduction near to boundary
556 if(smallstep < skin)
557 {
558 tlimit = stepmin;
559 insideskin = true;
560 }
561 else if(geomlimit < geombig)
562 {
563 if(geomlimit > skindepth)
564 {
565 if(tlimit > geomlimit-0.999*skindepth)
566 tlimit = geomlimit-0.999*skindepth;
567 }
568 else
569 {
570 insideskin = true;
571 if(tlimit > stepmin) tlimit = stepmin;
572 }
573 }
574
575 if(tlimit < stepmin) tlimit = stepmin;
576
577 if(tPathLength > tlimit) tPathLength = tlimit;
578
579 }
580 // for 'normal' simulation with or without magnetic field
581 // there no small step/single scattering at boundaries
582 else if(steppingAlgorithm == fUseSafety)
583 {
584 // compute presafety again if presafety <= 0 and no boundary
585 // i.e. when it is needed for optimization purposes
586 if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
587 presafety = ComputeSafety(sp->GetPosition(),tPathLength);
588
589 // is far from boundary
590 if(currentRange < presafety)
591 {
592 inside = true;
593 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
594 }
595
596 if(firstStep || stepStatus == fGeomBoundary)
597 {
598 rangeinit = currentRange;
599 fr = facrange;
600 // 9.1 like stepping for e+/e- only (not for muons,hadrons)
601 if(mass < masslimite)
602 {
603 if(lambda1 > currentRange)
604 rangeinit = lambda1;
605 if(lambda1 > lambdalimit)
606 fr *= 0.75+0.25*lambda1/lambdalimit;
607 }
608
609 //lower limit for tlimit
610 G4double rat = currentKinEnergy/MeV ;
611 rat = 1.e-3/(rat*(10.+rat)) ;
612 tlimitmin = 10.*lambda1*rat;
613 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
614 }
615 //step limit
616 tlimit = fr*rangeinit;
617
618 if(tlimit < facsafety*presafety)
619 tlimit = facsafety*presafety;
620
621 //lower limit for tlimit
622 if(tlimit < tlimitmin) tlimit = tlimitmin;
623
624 if(tPathLength > tlimit) tPathLength = tlimit;
625 }
626
627 // version similar to 7.1 (needed for some experiments)
628 else
629 {
630 if (stepStatus == fGeomBoundary)
631 {
632 if (currentRange > lambda1) tlimit = facrange*currentRange;
633 else tlimit = facrange*lambda1;
634
635 if(tlimit < tlimitmin) tlimit = tlimitmin;
636 if(tPathLength > tlimit) tPathLength = tlimit;
637 }
638 }
639 //G4cout << "tPathLength= " << tPathLength
640 // << " currentMinimalStep= " << currentMinimalStep << G4endl;
641 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
642}
643
644//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
645// taken from Ref.6
646G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double)
647{
648 firstStep = false;
649 par1 = -1. ;
650 par2 = par3 = 0. ;
651
652 // do the true -> geom transformation
653 zPathLength = tPathLength;
654
655 // z = t for very small tPathLength
656 if(tPathLength < tlimitminfix) { return zPathLength; }
657
658 // this correction needed to run MSC with eIoni and eBrem inactivated
659 // and makes no harm for a normal run
660 if(tPathLength > currentRange)
661 { tPathLength = currentRange; }
662
663 G4double tau = tPathLength/lambda1 ;
664
665 if ((tau <= tausmall) || insideskin) {
666 zPathLength = tPathLength;
667 if(zPathLength > lambda1) { zPathLength = lambda1; }
668 return zPathLength;
669 }
670
671 G4double zmean = tPathLength;
672 if (tPathLength < currentRange*dtrl) {
673 if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
674 else zmean = lambda1*(1.-exp(-tau));
675 } else if(currentKinEnergy < mass || tPathLength == currentRange) {
676 par1 = 1./currentRange ;
677 par2 = 1./(par1*lambda1) ;
678 par3 = 1.+par2 ;
679 if(tPathLength < currentRange)
680 zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
681 else
682 zmean = 1./(par1*par3) ;
683 } else {
684 G4double T1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
685
686 lambda11 = GetTransportMeanFreePath(particle,T1);
687
688 par1 = (lambda1-lambda11)/(lambda1*tPathLength) ;
689 par2 = 1./(par1*lambda1) ;
690 par3 = 1.+par2 ;
691 zmean = (1.-exp(par3*log(lambda11/lambda1)))/(par1*par3) ;
692 }
693
694 zPathLength = zmean ;
695 // sample z
696 if(samplez) {
697
698 const G4double ztmax = 0.99;
699 G4double zt = zmean/tPathLength ;
700
701 if (tPathLength > stepmin && zt < ztmax) {
702
703 G4double u,cz1;
704 if(zt >= 0.333333333) {
705
706 G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
707 cz1 = 1.+cz ;
708 G4double u0 = cz/cz1 ;
709 G4double grej ;
710 do {
711 u = exp(log(G4UniformRand())/cz1) ;
712 grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
713 } while (grej < G4UniformRand()) ;
714
715 } else {
716 cz1 = 1./zt-1.;
717 u = 1.-exp(log(G4UniformRand())/cz1) ;
718 }
719 zPathLength = tPathLength*u ;
720 }
721 }
722 if(zPathLength > lambda1) zPathLength = lambda1;
723 //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda1 << G4endl;
724
725 return zPathLength;
726}
727
728//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
729// taken from Ref.6
730G4double
731G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength)
732{
733 // step defined other than transportation
734 if(geomStepLength == zPathLength && tPathLength <= currentRange)
735 return tPathLength;
736
737 // t = z for very small step
738 zPathLength = geomStepLength;
739 tPathLength = geomStepLength;
740 if(geomStepLength < tlimitminfix) return tPathLength;
741
742 // recalculation
743 if((geomStepLength > lambda1*tausmall) && !insideskin)
744 {
745 if(par1 < 0.)
746 tPathLength = -lambda1*log(1.-geomStepLength/lambda1) ;
747 else
748 {
749 if(par1*par3*geomStepLength < 1.)
750 tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
751 else
752 tPathLength = currentRange;
753 }
754 }
755 if(tPathLength < geomStepLength) tPathLength = geomStepLength;
756 //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
757
758 return tPathLength;
759}
760
761//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
762//Total & first transport x sections for e-/e+ generated from ELSEPA code
763
764void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()
765{
766 G4String filename = "XSECTIONS.dat";
767
768 char* path = getenv("G4LEDATA");
769 if (!path)
770 {
771 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()","em0006",
772 FatalException,
773 "Environment variable G4LEDATA not defined");
774 return;
775 }
776
777 G4String pathString(path);
778 G4String dirFile = pathString + "/msc_GS/" + filename;
779 FILE *infile;
780 infile = fopen(dirFile,"r");
781 if (infile == 0)
782 {
783 G4ExceptionDescription ed;
784 ed << "Data file <" + dirFile + "> is not opened!" << G4endl;
785 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
786 "em0003",FatalException,ed);
787 return;
788 }
789
790 // Read parameters from tables and take logarithms
791 G4float aRead;
792 for(G4int i=0 ; i<106 ;i++){
793 if(1 == fscanf(infile,"%f\t",&aRead)) {
794 if(aRead > 0.0) { aRead = log(aRead); }
795 else { aRead = 0.0; }
796 } else {
797 G4ExceptionDescription ed;
798 ed << "Error reading <" + dirFile + "> loop #1 i= " << i << G4endl;
799 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
800 "em0003",FatalException,ed);
801 return;
802 }
803 ener[i]=aRead;
804 }
805 for(G4int j=0;j<103;j++){
806 for(G4int i=0;i<106;i++){
807 if(1 == fscanf(infile,"%f\t",&aRead)) {
808 if(aRead > 0.0) { aRead = log(aRead); }
809 else { aRead = 0.0; }
810 } else {
811 G4ExceptionDescription ed;
812 ed << "Error reading <" + dirFile + "> loop #2 j= " << j
813 << "; i= " << i << G4endl;
814 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
815 "em0003",FatalException,ed);
816 return;
817 }
818 TCSE[j][i]=aRead;
819 }
820 }
821 for(G4int j=0;j<103;j++){
822 for(G4int i=0;i<106;i++){
823 if(1 == fscanf(infile,"%f\t",&aRead)) {
824 if(aRead > 0.0) { aRead = log(aRead); }
825 else { aRead = 0.0; }
826 } else {
827 G4ExceptionDescription ed;
828 ed << "Error reading <" + dirFile + "> loop #3 j= " << j
829 << "; i= " << i << G4endl;
830 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
831 "em0003",FatalException,ed);
832 return;
833 }
834 FTCSE[j][i]=aRead;
835 }
836 }
837 for(G4int j=0;j<103;j++){
838 for(G4int i=0;i<106;i++){
839 if(1 == fscanf(infile,"%f\t",&aRead)) {
840 if(aRead > 0.0) { aRead = log(aRead); }
841 else { aRead = 0.0; }
842 } else {
843 G4ExceptionDescription ed;
844 ed << "Error reading <" + dirFile + "> loop #4 j= " << j
845 << "; i= " << i << G4endl;
846 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
847 "em0003",FatalException,ed);
848 return;
849 }
850 TCSP[j][i]=aRead;
851 }
852 }
853 for(G4int j=0;j<103;j++){
854 for(G4int i=0;i<106;i++){
855 if(1 == fscanf(infile,"%f\t",&aRead)) {
856 if(aRead > 0.0) { aRead = log(aRead); }
857 else { aRead = 0.0; }
858 } else {
859 G4ExceptionDescription ed;
860 ed << "Error reading <" + dirFile + "> loop #5 j= " << j
861 << "; i= " << i << G4endl;
862 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
863 "em0003",FatalException,ed);
864 return;
865 }
866 FTCSP[j][i]=aRead;
867 }
868 }
869
870 fclose(infile);
871
872}
873
874//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:44
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
fUseSafety
@ fUseSafety
Definition: G4MscStepLimitType.hh:50
fUseDistanceToBoundary
@ fUseDistanceToBoundary
Definition: G4MscStepLimitType.hh:51
G4StepStatus
G4StepStatus
Definition: G4StepStatus.hh:51
fGeomBoundary
@ fGeomBoundary
Definition: G4StepStatus.hh:54
G4double
double G4double
Definition: G4Types.hh:64
G4float
float G4float
Definition: G4Types.hh:65
G4int
int G4int
Definition: G4Types.hh:66
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit
virtual G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep)
Definition: G4GoudsmitSaundersonMscModel.cc:453
G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel
G4GoudsmitSaundersonMscModel(const G4String &nam="GoudsmitSaunderson")
Definition: G4GoudsmitSaundersonMscModel.cc:97
G4GoudsmitSaundersonMscModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4GoudsmitSaundersonMscModel.cc:128
G4GoudsmitSaundersonMscModel::ComputeTrueStepLength
virtual G4double ComputeTrueStepLength(G4double geomStepLength)
Definition: G4GoudsmitSaundersonMscModel.cc:731
G4GoudsmitSaundersonMscModel::SampleScattering
virtual G4ThreeVector & SampleScattering(const G4DynamicParticle *, G4double safety)
Definition: G4GoudsmitSaundersonMscModel.cc:154
G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double, G4double, G4double)
Definition: G4GoudsmitSaundersonMscModel.cc:139
G4GoudsmitSaundersonMscModel::ComputeGeomPathLength
virtual G4double ComputeGeomPathLength(G4double truePathLength)
Definition: G4GoudsmitSaundersonMscModel.cc:646
G4GoudsmitSaundersonTable::SampleTheta
G4double SampleTheta(G4double, G4double, G4double)
Definition: G4GoudsmitSaundersonTable.cc:203
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:107
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:189
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
G4Material::GetVecNbOfAtomsPerVolume
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:205
G4ParticleChangeForMSC::ProposeMomentumDirection
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4Step::GetPreStepPoint
G4StepPoint * GetPreStepPoint() const
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4Track::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4Track::GetStep
const G4Step * GetStep() const
G4VEmModel::SetCurrentCouple
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:370
G4VMscModel::dtrl
G4double dtrl
Definition: G4VMscModel.hh:180
G4VMscModel::GetDEDX
G4double GetDEDX(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:273
G4VMscModel::facrange
G4double facrange
Definition: G4VMscModel.hh:176
G4VMscModel::ComputeGeomLimit
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:256
G4VMscModel::samplez
G4bool samplez
Definition: G4VMscModel.hh:188
G4VMscModel::skin
G4double skin
Definition: G4VMscModel.hh:179
G4VMscModel::GetTransportMeanFreePath
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:332
G4VMscModel::GetEnergy
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:304
G4VMscModel::GetRange
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:288
G4VMscModel::ComputeSafety
G4double ComputeSafety(const G4ThreeVector &position, G4double limit)
Definition: G4VMscModel.hh:238
G4VMscModel::steppingAlgorithm
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:186
G4VMscModel::GetParticleChangeForMSC
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
Definition: G4VMscModel.cc:89
G4VMscModel::ConvertTrueToGeom
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:246
G4VMscModel::facsafety
G4double facsafety
Definition: G4VMscModel.hh:178
G4VMscModel::fDisplacement
G4ThreeVector fDisplacement
Definition: G4VMscModel.hh:185
G4VMscModel::facgeom
G4double facgeom
Definition: G4VMscModel.hh:177
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76