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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

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