source: trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc@ 1172

Last change on this file since 1172 was 1058, checked in by garnier, 17 years ago

file release beta
  • Property svn:executable set to *
File size: 25.1 KB
Line 
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: G4GoudsmitSaundersonMscModel.cc,v 1.7 2009/04/20 19:22:29 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33// File name: G4GoudsmitSaundersonMscModel
34//
35// Author: Omrane Kadri
36//
37// Creation date: 20.02.2009
38//
39// Modifications:
40// 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
41//
42// 15.04.2009 O.Kadri: cleanup: discard no scattering and single scattering theta
43// sampling from SampleCosineTheta() which means the splitting
44// step into two sub-steps occur only for msc regime
45//
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48//REFERENCES:
49//Ref.1:E. Benedito et al.,"Mixed simulation ... cross-sections", NIMB 174 (2001) pp 91-110;
50//Ref.2:I. Kawrakow et al.,"On the condensed ... transport",NIMB 142 (1998) pp 253-280;
51//Ref.3:I. Kawrakow et al.,"On the representation ... calculations",NIMB 134 (1998) pp 325-336;
52//Ref.4:Bielajew et al.,".....", NIMB 173 (2001) 332-343;
53//Ref.5:F. Salvat et al.,"ELSEPA--Dirac partial ...molecules", Comp.Phys.Comm.165 (2005) pp 157-190;
54//Ref.6:G4UrbanMscModel G4_v9.1Ref09;
55//Ref.7:G4eCoulombScatteringModel G4_v9.1Ref09.
56//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57
58#include "G4GoudsmitSaundersonMscModel.hh"
59#include "G4GoudsmitSaundersonTable.hh"
60
61#include "G4ParticleChangeForMSC.hh"
62#include "G4MaterialCutsCouple.hh"
63#include "G4DynamicParticle.hh"
64#include "G4DataInterpolation.hh"
65#include "G4Electron.hh"
66#include "G4Positron.hh"
67
68#include "G4LossTableManager.hh"
69#include "G4Track.hh"
70#include "G4PhysicsTable.hh"
71#include "Randomize.hh"
72#include "G4Poisson.hh"
73
74using namespace std;
75
76G4double G4GoudsmitSaundersonMscModel::ener[] = {-1.};
77G4double G4GoudsmitSaundersonMscModel::TCSE[103][106] ;
78G4double G4GoudsmitSaundersonMscModel::FTCSE[103][106] ;
79G4double G4GoudsmitSaundersonMscModel::TCSP[103][106] ;
80G4double G4GoudsmitSaundersonMscModel::FTCSP[103][106] ;
81
82//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam)
84 : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(GeV),isInitialized(false)
85{
86 fr=0.02,rangeinit=0.,masslimite=0.6*MeV;
87 particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm;
88 tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig;
89 tlimitminfix=1.e-6*mm;stepmin=tlimitminfix;lambdalimit=1.*mm;smallstep=1.e10;
90 theManager=G4LossTableManager::Instance();
91 inside=false;insideskin=false;
92 samplez=false;
93
94 GSTable = new G4GoudsmitSaundersonTable();
95
96 if(ener[0] < 0.0){
97 G4cout << "### G4GoudsmitSaundersonMscModel loading ELSEPA data" << G4endl;
98 LoadELSEPAXSections();
99 }
100}
101//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel()
103{
104 delete GSTable;
105}
106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107void G4GoudsmitSaundersonMscModel::Initialise(const G4ParticleDefinition* p,
108 const G4DataVector&)
109{
110 skindepth=skin*stepmin;
111 SetParticle(p);
112 if(isInitialized) return;
113 fParticleChange = GetParticleChangeForMSC();
114 InitialiseSafetyHelper();
115 isInitialized=true;
116}
117//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
118
119G4double G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
120 G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
121{
122 //Build cross section table : Taken from Ref.7
123 G4double cs=0.0;
124 G4double kinEnergy = kineticEnergy;
125 if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
126 if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
127
128 //value0=Lambda0;value1=Lambda1
129 G4double value0,value1;
130 CalculateIntegrals(p,Z,kinEnergy,value0,value1);
131
132 if(value1 > 0.0) cs = 1./value1;
133
134 return cs;
135}
136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137
138void G4GoudsmitSaundersonMscModel::SampleScattering(const G4DynamicParticle* dynParticle,
139 G4double safety)
140{
141 G4double kineticEnergy = dynParticle->GetKineticEnergy();
142 if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)) return ;
143
144 G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
145 G4double phi1,phi2,cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
146 G4double q1,Gamma,Eta,delta,nu,nu0,nu1,nu2,nu_interm;
147
148 ///////////////////////////////////////////
149 // Effective energy and path-length from Eq. 4.7.15+16 of Ref.4
150 G4double eloss = theManager->GetEnergy(particle,tPathLength,currentCouple);
151 if(eloss>0.5*kineticEnergy)return;
152 G4double ee = kineticEnergy - 0.5*eloss;
153 G4double ttau = ee/electron_mass_c2;
154 G4double ttau2 = ttau*ttau;
155 G4double epsilonpp= eloss/ee;
156 G4double temp2 = 0.166666*(4+ttau*(6+ttau*(7+ttau*(4+ttau))))*(epsilonpp/(ttau+1)/(ttau+2))*(epsilonpp/(ttau+1)/(ttau+2));
157 G4double cst1=epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72);
158
159 kineticEnergy *= (1 - cst1);
160 tPathLength *= (1 - temp2);
161 ///////////////////////////////////////////
162 // additivity rule for mixture xsection calculation
163 const G4Material* mat = currentCouple->GetMaterial();
164 G4int nelm = mat->GetNumberOfElements();
165 const G4ElementVector* theElementVector = mat->GetElementVector();
166 const G4double* theFraction = mat->GetFractionVector();
167 G4double atomPerVolume = mat->GetTotNbOfAtomsPerVolume();
168 G4double scrA,llambda0,llambda1;
169 scrA=0.0;llambda0 =0.;llambda1=0.;
170 for(G4int i=0;i<nelm;i++)
171 {
172 G4double l0,l1;
173 CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,l0,l1);
174 llambda0 += (theFraction[i]/l0);
175 llambda1 += (theFraction[i]/l1);
176 }
177 if(llambda0>DBL_MIN)llambda0 =1./llambda0;
178 if(llambda1>DBL_MIN)llambda1 =1./llambda1;
179 G4double g1=llambda0/llambda1;
180 G4double x1,x0;
181
182 x0=g1/2.;
183 do
184 {
185 x1 = x0-(x0*((1.+x0)*std::log(1.+1./x0)-1.0)-g1/2.)/( (1.+2.*x0)*std::log(1.+1./x0)-2.0);// x1=x0-f(x0)/f'(x0)
186 delta = std::abs( x1 - x0 );
187 x0 = x1; // new approximation becomes the old approximation for the next iteration
188 } while (delta > 1e-10);
189 scrA = x1;
190
191 G4double us=0.0,vs=0.0,ws=1.0,x_coord=0.0,y_coord=0.0,z_coord=1.0;
192 G4double lambdan=0.;
193 if(llambda0>0.)lambdan=atomPerVolume*tPathLength/llambda0;
194 if((lambdan<=1.0e-12)||(lambdan>1.0e+5))return;
195 G4bool noscatt=false;
196 G4bool singlescatt=false;
197 G4bool mscatt=false;
198
199 G4double epsilon1=G4UniformRand();
200 if(epsilon1<(exp(-lambdan)))noscatt=true;// no scattering
201 else if(epsilon1<((1.+lambdan)*exp(-lambdan)))//single scattering
202 {singlescatt=true;
203 ws=G4UniformRand();
204 ws= 1.-2.*scrA*ws/(1.-ws + scrA);
205 G4double phi0=twopi*G4UniformRand();
206 us=sqrt(1.-ws*ws)*cos(phi0);
207 vs=sqrt(1.-ws*ws)*sin(phi0);
208 G4double rr=G4UniformRand();
209 x_coord=(rr*us);
210 y_coord=(rr*vs);
211 z_coord=((1.-rr)+rr*ws);
212 }
213 else
214 {mscatt=true;
215 // Ref.2 subsection 4.4 "The best solution found"
216 // Sample first substep scattering angle
217 SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
218 phi1 = twopi*G4UniformRand();
219 cosPhi1 = cos(phi1);
220 sinPhi1 = sin(phi1);
221
222 // Sample second substep scattering angle
223 SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
224 phi2 = twopi*G4UniformRand();
225 cosPhi2 = cos(phi2);
226 sinPhi2 = sin(phi2);
227
228 // Scattering direction
229 us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
230 vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
231 ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2;
232 }
233 G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
234 G4ThreeVector newDirection(us,vs,ws);
235 newDirection.rotateUz(oldDirection);
236 fParticleChange->ProposeMomentumDirection(newDirection);
237
238 if((safety > tlimitminfix)&&(latDisplasment))
239 {
240 // Scattering coordinates
241 if(mscatt)
242 {
243 if(scrA<DBL_MIN)scrA=DBL_MIN;
244 q1 = 2.*scrA*((1. + scrA)*log(1. + 1./scrA) - 1.);
245 if(q1<DBL_MIN)q1=DBL_MIN;
246 Gamma = 6.*scrA*(1. + scrA)*((1. + 2.*scrA)*log(1. + 1./scrA) - 2.)/q1;
247 Eta = atomPerVolume*tPathLength/llambda1;
248 delta = 0.90824829 - Eta*(0.102062073-Gamma*0.026374715);
249
250 nu = G4UniformRand();
251 nu = std::sqrt(nu);
252 nu0 = (1.0 - nu)/2.;
253 nu1 = nu*delta;
254 nu2 = nu*(1.0-delta);
255 nu_interm = 1.0 - nu0 - nu1 - nu2;
256 x_coord=(nu1*sinTheta1*cosPhi1+nu2*sinTheta2*(cosPhi1*cosPhi2-cosTheta1*sinPhi1*sinPhi2)+nu_interm*us);
257 y_coord=(nu1*sinTheta1*sinPhi1+nu2*sinTheta2*(sinPhi1*cosPhi2+cosTheta1*cosPhi1*sinPhi2)+nu_interm*vs);
258 z_coord=(nu0+nu1*cosTheta1+nu2*cosTheta2+ nu_interm*ws) ;
259 }
260 G4double r=sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
261
262 G4double check= 1.- tPathLength/zPathLength;
263 if(check<=0.) return;
264 else if(r>check) {r=check;x_coord /=r;y_coord /=r;z_coord /=r;}
265
266 r *=tPathLength;
267
268 if(r > tlimitminfix) {
269
270 G4ThreeVector latDirection = G4ThreeVector(x_coord*tPathLength,y_coord*tPathLength,z_coord*tPathLength);
271 latDirection.rotateUz(oldDirection);
272
273 ComputeDisplacement(fParticleChange, latDirection, r, safety);
274 }
275 }
276}
277
278//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
279void G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA,
280 G4double &cost, G4double &sint)
281{
282 G4double u,Qn1,r1,tet;
283 G4double xi=0.;
284 Qn1=2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
285
286 if((lambdan<1.)||(Qn1<0.001))//plural scatt. or small angle scatt.
287 {
288 G4double xi1,lambdai=0.;
289 G4int i=0;
290 do {xi1=G4UniformRand();
291 lambdai -=std::log(xi1);
292 xi +=2.*scrA*xi1/(1.-xi1 + scrA);
293 i++;
294 }while((lambdai<lambdan)&&(i<30));
295
296 }
297 else {
298 if(Qn1>0.5)xi=2.*G4UniformRand();//isotropic distribution
299 else{// procedure described by Benedito in Ref.1
300 do{r1=G4UniformRand();
301 u=GSTable->SampleTheta(lambdan,scrA,G4UniformRand());
302 xi = 2.*u;
303 tet=acos(1.-xi);
304 }while(tet*r1*r1>sin(tet));
305 }
306 }
307
308
309 if(xi<0.)xi=0.;
310 if(xi>2.)xi=2.;
311 cost=(1. - xi);
312 sint=sqrt(xi*(2.-xi));
313
314}
315
316//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
317// Cubic spline log-log interpolation of Lambda0 and Lambda1
318// Screening parameter calculated according to Eq. 37 of Ref.1
319void G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z,
320 G4double kinEnergy,G4double &Lam0,
321 G4double &Lam1)
322{
323 //////// BEGIN OF: LAMBDA CALCULATION ////////////////////////////////
324 G4double TCSEForThisAtom[106],FTCSEForThisAtom[106],TCSPForThisAtom[106],FTCSPForThisAtom[106];
325 G4double summ00=0.0;
326 G4double summ10=0.0;
327 G4double InterpolatedValue=0.0;
328
329
330 G4int iZ = G4int(Z);
331 if(iZ > 103) iZ = 103;
332 for(G4int i=0;i<106;i++)
333 {
334 TCSEForThisAtom[i]=TCSE[iZ-1][i];FTCSEForThisAtom[i]=FTCSE[iZ-1][i];
335 TCSPForThisAtom[i]=TCSP[iZ-1][i];FTCSPForThisAtom[i]=FTCSP[iZ-1][i];
336 }
337
338 G4double kineticE = kinEnergy;
339 if(kineticE<lowKEnergy)kineticE=lowKEnergy;
340 if(kineticE>highKEnergy)kineticE=highKEnergy;
341 kineticE /= eV;
342
343 if(p==G4Electron::Electron())
344 {
345 MyValue= new G4DataInterpolation(ener,TCSEForThisAtom,106,0.0,0);
346 InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
347 delete MyValue;
348 summ00 = std::exp(InterpolatedValue);
349 MyValue= new G4DataInterpolation(ener,FTCSEForThisAtom,106,0.0,0);
350 InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
351 delete MyValue;
352 summ10 = std::exp(InterpolatedValue);
353 }
354 if(p==G4Positron::Positron())
355 {
356 MyValue= new G4DataInterpolation(ener,TCSPForThisAtom,106,0.0,0);
357 InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
358 delete MyValue;
359 summ00 = std::exp(InterpolatedValue);
360 MyValue= new G4DataInterpolation(ener,FTCSPForThisAtom,106,0.0,0);
361 InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
362 delete MyValue;
363 summ10 = std::exp(InterpolatedValue);
364 }
365
366 summ00 *=barn;
367 summ10 *=barn;
368
369 Lam0=1./((1.+1./Z)*summ00);
370 Lam1=1./((1.+1./Z)*summ10);
371
372}
373
374//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
375//t->g->t step transformations taken from Ref.6
376G4double G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track,
377 G4PhysicsTable* theTable,
378 G4double currentMinimalStep)
379{
380 tPathLength = currentMinimalStep;
381 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
382 G4StepStatus stepStatus = sp->GetStepStatus();
383
384 const G4DynamicParticle* dp = track.GetDynamicParticle();
385
386 if(stepStatus == fUndefined) {
387 inside = false;
388 insideskin = false;
389 tlimit = geombig;
390 SetParticle( dp->GetDefinition() );
391 }
392
393 theLambdaTable = theTable;
394 currentCouple = track.GetMaterialCutsCouple();
395 currentMaterialIndex = currentCouple->GetIndex();
396 currentKinEnergy = dp->GetKineticEnergy();
397 currentRange =
398 theManager->GetRangeFromRestricteDEDX(particle,currentKinEnergy,currentCouple);
399
400 lambda1 = GetLambda(currentKinEnergy);
401
402 // stop here if small range particle
403 if(inside) return tPathLength;
404
405 if(tPathLength > currentRange) tPathLength = currentRange;
406
407 G4double presafety = sp->GetSafety();
408
409 // far from geometry boundary
410 if(currentRange < presafety)
411 {
412 inside = true;
413 return tPathLength;
414 }
415
416 // standard version
417 //
418 if (steppingAlgorithm == fUseDistanceToBoundary)
419 {
420 //compute geomlimit and presafety
421 G4double geomlimit = ComputeGeomLimit(track, presafety, tPathLength);
422
423 // is far from boundary
424 if(currentRange <= presafety)
425 {
426 inside = true;
427 return tPathLength;
428 }
429
430 smallstep += 1.;
431 insideskin = false;
432
433 if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
434 {
435 rangeinit = currentRange;
436 if(stepStatus == fUndefined) smallstep = 1.e10;
437 else smallstep = 1.;
438
439 // constraint from the geometry
440 if((geomlimit < geombig) && (geomlimit > geommin))
441 {
442 if(stepStatus == fGeomBoundary)
443 tgeom = geomlimit/facgeom;
444 else
445 tgeom = 2.*geomlimit/facgeom;
446 }
447 else
448 tgeom = geombig;
449
450 //define stepmin here (it depends on lambda!)
451 //rough estimation of lambda_elastic/lambda_transport
452 G4double rat = currentKinEnergy/MeV ;
453 rat = 1.e-3/(rat*(10.+rat)) ;
454 //stepmin ~ lambda_elastic
455 stepmin = rat*lambda1;
456 skindepth = skin*stepmin;
457
458 //define tlimitmin
459 tlimitmin = 10.*stepmin;
460 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
461
462 }
463
464 //step limit
465 tlimit = facrange*rangeinit;
466 if(tlimit < facsafety*presafety)
467 tlimit = facsafety*presafety;
468
469 //lower limit for tlimit
470 if(tlimit < tlimitmin) tlimit = tlimitmin;
471
472 if(tlimit > tgeom) tlimit = tgeom;
473
474 // shortcut
475 if((tPathLength < tlimit) && (tPathLength < presafety) &&
476 (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
477 return tPathLength;
478
479 // step reduction near to boundary
480 if(smallstep < skin)
481 {
482 tlimit = stepmin;
483 insideskin = true;
484 }
485 else if(geomlimit < geombig)
486 {
487 if(geomlimit > skindepth)
488 {
489 if(tlimit > geomlimit-0.999*skindepth)
490 tlimit = geomlimit-0.999*skindepth;
491 }
492 else
493 {
494 insideskin = true;
495 if(tlimit > stepmin) tlimit = stepmin;
496 }
497 }
498
499 if(tlimit < stepmin) tlimit = stepmin;
500
501 if(tPathLength > tlimit) tPathLength = tlimit ;
502
503 }
504 // for 'normal' simulation with or without magnetic field
505 // there no small step/single scattering at boundaries
506 else if(steppingAlgorithm == fUseSafety)
507 {
508 // compute presafety again if presafety <= 0 and no boundary
509 // i.e. when it is needed for optimization purposes
510 if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
511 presafety = ComputeSafety(sp->GetPosition(),tPathLength);
512
513 // is far from boundary
514 if(currentRange < presafety)
515 {
516 inside = true;
517 return tPathLength;
518 }
519
520 if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
521 {
522 rangeinit = currentRange;
523 fr = facrange;
524 // 9.1 like stepping for e+/e- only (not for muons,hadrons)
525 if(mass < masslimite)
526 {
527 if(lambda1 > currentRange)
528 rangeinit = lambda1;
529 if(lambda1 > lambdalimit)
530 fr *= 0.75+0.25*lambda1/lambdalimit;
531 }
532
533 //lower limit for tlimit
534 G4double rat = currentKinEnergy/MeV ;
535 rat = 1.e-3/(rat*(10.+rat)) ;
536 tlimitmin = 10.*lambda1*rat;
537 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
538 }
539 //step limit
540 tlimit = fr*rangeinit;
541
542 if(tlimit < facsafety*presafety)
543 tlimit = facsafety*presafety;
544
545 //lower limit for tlimit
546 if(tlimit < tlimitmin) tlimit = tlimitmin;
547
548 if(tPathLength > tlimit) tPathLength = tlimit;
549 }
550
551 // version similar to 7.1 (needed for some experiments)
552 else
553 {
554 if (stepStatus == fGeomBoundary)
555 {
556 if (currentRange > lambda1) tlimit = facrange*currentRange;
557 else tlimit = facrange*lambda1;
558
559 if(tlimit < tlimitmin) tlimit = tlimitmin;
560 if(tPathLength > tlimit) tPathLength = tlimit;
561 }
562 }
563 return tPathLength ;
564}
565
566//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
567G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double)
568{
569 par1 = -1. ;
570 par2 = par3 = 0. ;
571
572 // do the true -> geom transformation
573 zPathLength = tPathLength;
574
575 // z = t for very small tPathLength
576 if(tPathLength < tlimitminfix) return zPathLength;
577
578 // this correction needed to run MSC with eIoni and eBrem inactivated
579 // and makes no harm for a normal run
580 if(tPathLength > currentRange)
581 tPathLength = currentRange ;
582
583 G4double tau = tPathLength/lambda1 ;
584
585 if ((tau <= tausmall) || insideskin) {
586 zPathLength = tPathLength;
587 if(zPathLength > lambda1) zPathLength = lambda1;
588 return zPathLength;
589 }
590
591 G4double zmean = tPathLength;
592 if (tPathLength < currentRange*dtrl) {
593 if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
594 else zmean = lambda1*(1.-exp(-tau));
595 } else if(currentKinEnergy < mass) {
596 par1 = 1./currentRange ;
597 par2 = 1./(par1*lambda1) ;
598 par3 = 1.+par2 ;
599 if(tPathLength < currentRange)
600 zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
601 else
602 zmean = 1./(par1*par3) ;
603 } else {
604 G4double T1 = theManager->GetEnergy(particle,currentRange-tPathLength,currentCouple);
605
606 lambda11 = GetLambda(T1);
607
608 par1 = (lambda1-lambda11)/(lambda1*tPathLength) ;
609 par2 = 1./(par1*lambda1) ;
610 par3 = 1.+par2 ;
611 zmean = (1.-exp(par3*log(lambda11/lambda1)))/(par1*par3) ;
612 }
613
614 zPathLength = zmean ;
615 // sample z
616 if(samplez)
617 {
618 const G4double ztmax = 0.99, onethird = 1./3. ;
619 G4double zt = zmean/tPathLength ;
620
621 if (tPathLength > stepmin && zt < ztmax)
622 {
623 G4double u,cz1;
624 if(zt >= onethird)
625 {
626 G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
627 cz1 = 1.+cz ;
628 G4double u0 = cz/cz1 ;
629 G4double grej ;
630 do {
631 u = exp(log(G4UniformRand())/cz1) ;
632 grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
633 } while (grej < G4UniformRand()) ;
634 }
635 else
636 {
637 cz1 = 1./zt-1.;
638 u = 1.-exp(log(G4UniformRand())/cz1) ;
639 }
640 zPathLength = tPathLength*u ;
641 }
642 }
643 if(zPathLength > lambda1) zPathLength = lambda1;
644
645 return zPathLength;
646}
647
648//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
649
650G4double G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength)
651{
652 // step defined other than transportation
653 if(geomStepLength == zPathLength && tPathLength <= currentRange)
654 return tPathLength;
655
656 // t = z for very small step
657 zPathLength = geomStepLength;
658 tPathLength = geomStepLength;
659 if(geomStepLength < tlimitminfix) return tPathLength;
660
661 // recalculation
662 if((geomStepLength > lambda1*tausmall) && !insideskin)
663 {
664 if(par1 < 0.)
665 tPathLength = -lambda1*log(1.-geomStepLength/lambda1) ;
666 else
667 {
668 if(par1*par3*geomStepLength < 1.)
669 tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
670 else
671 tPathLength = currentRange;
672 }
673 }
674 if(tPathLength < geomStepLength) tPathLength = geomStepLength;
675
676 return tPathLength;
677}
678
679//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
680void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()
681{
682 ///////////////////////////////////////
683 //Total & first transport x sections of e-/e+ from ELSEPA code
684 G4String filename = "XSECTIONS.dat";
685
686 char* path = getenv("G4LEDATA");
687 if (!path)
688 {
689 G4String excep = "G4GoudsmitSaundersonTable: G4LEDATA environment variable not set";
690 G4Exception(excep);
691 }
692
693 G4String pathString(path);
694 G4String dirFile = pathString + "/msc_GS/" + filename;
695 FILE *infile;
696 infile = fopen(dirFile,"r");
697 if (infile == 0)
698 {
699 G4String excep = "G4GoudsmitSaunderson - data files: " + dirFile + " not found";
700 G4Exception(excep);
701 }
702
703 // Read parameters from tables and take logarithms
704 G4float aRead;
705 for(G4int i=0 ; i<106 ;i++){
706 fscanf(infile,"%f\t",&aRead);
707 if(aRead > 0.0) aRead = std::log(aRead);
708 else aRead = 0.0;
709 ener[i]=aRead;
710 }
711 for(G4int j=0;j<103;j++){
712 for(G4int i=0;i<106;i++){
713 fscanf(infile,"%f\t",&aRead);
714 if(aRead > 0.0) aRead = std::log(aRead);
715 else aRead = 0.0;
716 TCSE[j][i]=aRead;
717 }
718 }
719 for(G4int j=0;j<103;j++){
720 for(G4int i=0;i<106;i++){
721 fscanf(infile,"%f\t",&aRead);
722 if(aRead > 0.0) aRead = std::log(aRead);
723 else aRead = 0.0;
724 FTCSE[j][i]=aRead;
725 }
726 }
727 for(G4int j=0;j<103;j++){
728 for(G4int i=0;i<106;i++){
729 fscanf(infile,"%f\t",&aRead);
730 if(aRead > 0.0) aRead = std::log(aRead);
731 else aRead = 0.0;
732 TCSP[j][i]=aRead;
733 }
734 }
735 for(G4int j=0;j<103;j++){
736 for(G4int i=0;i<106;i++){
737 fscanf(infile,"%f\t",&aRead);
738 if(aRead > 0.0) aRead = std::log(aRead);
739 else aRead = 0.0;
740 FTCSP[j][i]=aRead;
741 }
742 }
743
744 fclose(infile);
745 //End loading XSections and Energies
746
747}
748
749//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.