Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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