source: trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.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
File size: 24.4 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: G4WentzelVIModel.cc,v 1.37 2009/10/28 10:14:13 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4WentzelVIModel
35//
36// Author: V.Ivanchenko
37//
38// Creation date: 09.04.2008 from G4MuMscModel
39//
40// Modifications:
41//
42//
43// Class Description:
44//
45// Implementation of the model of multiple scattering based on
46// G.Wentzel, Z. Phys. 40 (1927) 590.
47// H.W.Lewis, Phys Rev 78 (1950) 526.
48// J.M. Fernandez-Varea et al., NIM B73 (1993) 447.
49// L.Urban, CERN-OPEN-2006-077.
50
51// -------------------------------------------------------------------
52//
53
54//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57#include "G4WentzelVIModel.hh"
58#include "Randomize.hh"
59#include "G4LossTableManager.hh"
60#include "G4ParticleChangeForMSC.hh"
61#include "G4PhysicsTableHelper.hh"
62#include "G4ElementVector.hh"
63#include "G4ProductionCutsTable.hh"
64#include "G4PhysicsLogVector.hh"
65#include "G4Electron.hh"
66#include "G4Positron.hh"
67#include "G4Proton.hh"
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70
71G4double G4WentzelVIModel::ScreenRSquare[] = {0.0};
72G4double G4WentzelVIModel::FormFactor[] = {0.0};
73
74using namespace std;
75
76G4WentzelVIModel::G4WentzelVIModel(const G4String& nam) :
77 G4VMscModel(nam),
78 theLambdaTable(0),
79 theLambda2Table(0),
80 numlimit(0.2),
81 nbins(60),
82 nwarnings(0),
83 nwarnlimit(50),
84 currentCouple(0),
85 cosThetaMin(1.0),
86 q2Limit(TeV*TeV),
87 alpha2(fine_structure_const*fine_structure_const),
88 isInitialized(false),
89 inside(false)
90{
91 invsqrt12 = 1./sqrt(12.);
92 tlimitminfix = 1.e-6*mm;
93 theManager = G4LossTableManager::Instance();
94 fNistManager = G4NistManager::Instance();
95 theElectron = G4Electron::Electron();
96 thePositron = G4Positron::Positron();
97 theProton = G4Proton::Proton();
98 lowEnergyLimit = 0.1*keV;
99 G4double p0 = electron_mass_c2*classic_electr_radius;
100 coeff = twopi*p0*p0;
101 tkin = targetZ = mom2 = DBL_MIN;
102 ecut = etag = DBL_MAX;
103 particle = 0;
104 nelments = 5;
105 xsecn.resize(nelments);
106 prob.resize(nelments);
107
108 // Thomas-Fermi screening radii
109 // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
110
111 if(0.0 == ScreenRSquare[0]) {
112 G4double a0 = electron_mass_c2/0.88534;
113 G4double constn = 6.937e-6/(MeV*MeV);
114
115 ScreenRSquare[0] = alpha2*a0*a0;
116 for(G4int j=1; j<100; j++) {
117 G4double x = a0*fNistManager->GetZ13(j);
118 ScreenRSquare[j] = alpha2*x*x;
119 x = fNistManager->GetA27(j);
120 FormFactor[j] = constn*x*x;
121 }
122 }
123}
124
125//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126
127G4WentzelVIModel::~G4WentzelVIModel()
128{}
129
130//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131
132void G4WentzelVIModel::Initialise(const G4ParticleDefinition* p,
133 const G4DataVector& cuts)
134{
135 // reset parameters
136 SetupParticle(p);
137 tkin = targetZ = mom2 = 0.0;
138 ecut = etag = DBL_MAX;
139 currentRange = 0.0;
140 cosThetaMax = cos(PolarAngleLimit());
141 currentCuts = &cuts;
142
143 // set values of some data members
144 if(!isInitialized) {
145 isInitialized = true;
146 fParticleChange = GetParticleChangeForMSC();
147 InitialiseSafetyHelper();
148 }
149}
150
151//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152
153G4double G4WentzelVIModel::ComputeCrossSectionPerAtom(
154 const G4ParticleDefinition* p,
155 G4double kinEnergy,
156 G4double Z, G4double,
157 G4double cutEnergy, G4double)
158{
159 SetupParticle(p);
160 if(kinEnergy < lowEnergyLimit) return 0.0;
161 SetupKinematic(kinEnergy, cutEnergy);
162 SetupTarget(Z, kinEnergy);
163 G4double xsec = ComputeTransportXSectionPerAtom();
164 /*
165 G4cout << "CS: e= " << tkin << " cosEl= " << cosTetMaxElec2
166 << " cosN= " << cosTetMaxNuc2 << " xsec(bn)= " << xsec/barn
167 << " " << particle->GetParticleName() << G4endl;
168 */
169 return xsec;
170}
171
172//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
173
174G4double G4WentzelVIModel::ComputeTransportXSectionPerAtom()
175{
176 G4double xSection = 0.0;
177 G4double x, y, x1, x2, x3, x4;
178
179 // scattering off electrons
180 if(cosTetMaxElec2 < 1.0) {
181 x = (1.0 - cosTetMaxElec2)/screenZ;
182 if(x < numlimit) y = 0.5*x*x*(1.0 - 1.3333333*x + 1.5*x*x);
183 else y = log(1.0 + x) - x/(1.0 + x);
184 if(y < 0.0) {
185 nwarnings++;
186 if(nwarnings < nwarnlimit /*&& y < -1.e-10*/) {
187 G4cout << "Electron scattering <0 for L1 " << y
188 << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
189 << " Z= " << targetZ << " "
190 << particle->GetParticleName() << G4endl;
191 G4cout << " z= " << 1.0-cosTetMaxElec2 << " screenZ= " << screenZ
192 << " x= " << x << G4endl;
193 }
194 y = 0.0;
195 }
196 xSection = y;
197 }
198 /*
199 G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
200 << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
201 << " cut(MeV)= " << ecut/MeV
202 << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
203 << " zmaxN= " << (1.0 - cosTetMaxNuc2)/screenZ
204 << " costm= " << cosTetMaxNuc2 << G4endl;
205 */
206 // scattering off nucleus
207 if(cosTetMaxNuc2 < 1.0) {
208 x = 1.0 - cosTetMaxNuc2;
209 x1 = screenZ*formfactA;
210 x2 = 1.0 - x1;
211 x3 = x/screenZ;
212 x4 = formfactA*x;
213 // low-energy limit
214 if(x3 < numlimit && x1 < numlimit) {
215 y = 0.5*x3*x3*(1.0 - 1.3333333*x3 + 1.5*x3*x3 - 1.5*x1
216 + 3.0*x1*x1 + 2.666666*x3*x1)/(x2*x2*x2);
217 // high energy limit
218 } else if(x2 <= 0.0) {
219 x4 = x1*(1.0 + x3);
220 y = x3*(1.0 + 0.5*x3 - (2.0 - x1)*(1.0 + x3 + x3*x3/3.0)/x4)/(x4*x4);
221 // middle energy
222 } else {
223 y = ((1.0 + x1)*x2*log((1. + x3)/(1. + x4))
224 - x3/(1. + x3) - x4/(1. + x4))/(x2*x2);
225 }
226 //G4cout << "y= " << y << " x1= " <<x1<<" x2= " <<x2
227 // <<" x3= "<<x3<<" x4= " << x4<<G4endl;
228 if(y < 0.0) {
229 nwarnings++;
230 if(nwarnings < nwarnlimit /*&& y < -1.e-10*/) {
231 G4cout << "Nuclear scattering <0 for L1 " << y
232 << " e(MeV)= " << tkin << " Z= " << targetZ << " "
233 << particle->GetParticleName() << G4endl;
234 G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
235 << " x= " << " x1= " << x1 << " x2= " << x2
236 << " x3= " << x3 << " x4= " << x4 <<G4endl;
237 }
238 y = 0.0;
239 }
240 xSection += y*targetZ;
241 }
242 xSection *= kinFactor;
243 /*
244 G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
245 << " screenZ= " << screenZ << " formF= " << formfactA
246 << " for " << particle->GetParticleName()
247 << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
248 << " x= " << x
249 << G4endl;
250 */
251 return xSection;
252}
253
254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255
256G4double G4WentzelVIModel::ComputeTruePathLengthLimit(
257 const G4Track& track,
258 G4PhysicsTable* theTable,
259 G4double currentMinimalStep)
260{
261 G4double tlimit = currentMinimalStep;
262 const G4DynamicParticle* dp = track.GetDynamicParticle();
263 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
264 G4StepStatus stepStatus = sp->GetStepStatus();
265
266 // initialisation for 1st step
267 if(stepStatus == fUndefined) {
268 inside = false;
269 SetupParticle(dp->GetDefinition());
270 theLambdaTable = theTable;
271 }
272
273 // initialisation for each step, lambda may be computed from scratch
274 preKinEnergy = dp->GetKineticEnergy();
275 DefineMaterial(track.GetMaterialCutsCouple());
276 lambda0 = GetLambda(preKinEnergy);
277 currentRange =
278 theManager->GetRangeFromRestricteDEDX(particle,preKinEnergy,currentCouple);
279
280 // extra check for abnormal situation
281 // this check needed to run MSC with eIoni and eBrem inactivated
282 if(tlimit > currentRange) tlimit = currentRange;
283
284 // stop here if small range particle
285 if(inside) return tlimit;
286
287 // pre step
288 G4double presafety = sp->GetSafety();
289
290 // compute presafety again if presafety <= 0 and no boundary
291 // i.e. when it is needed for optimization purposes
292 if(stepStatus != fGeomBoundary && presafety < tlimitminfix)
293 presafety = ComputeSafety(sp->GetPosition(), tlimit);
294 /*
295 G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit tlimit= "
296 <<tlimit<<" safety= " << presafety
297 << " range= " <<currentRange<<G4endl;
298 */
299 // far from geometry boundary
300 if(currentRange < presafety) {
301 inside = true;
302
303 // limit mean scattering angle
304 } else {
305 G4double rlimit = facrange*lambda0;
306 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
307 if(rcut > rlimit) rlimit = std::pow(rcut*rcut*rlimit,0.33333333);
308 rlimit = std::min(rlimit, facgeom*currentMaterial->GetRadlen());
309 if(rlimit < tlimit) tlimit = rlimit;
310 }
311 /*
312 G4cout << particle->GetParticleName() << " e= " << preKinEnergy
313 << " L0= " << lambda0 << " R= " << currentRange
314 << "tlimit= " << tlimit
315 << " currentMinimalStep= " << currentMinimalStep << G4endl;
316 */
317 return tlimit;
318}
319
320//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
321
322G4double G4WentzelVIModel::ComputeGeomPathLength(G4double truelength)
323{
324 tPathLength = truelength;
325 zPathLength = tPathLength;
326 lambdaeff = lambda0;
327
328 if(lambda0 > 0.0) {
329 G4double tau = tPathLength/lambda0;
330 //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
331 // << " lambda0= " << lambda0 << " tau= " << tau << G4endl;
332 // small step
333 if(tau < numlimit) {
334 zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
335
336 // medium step
337 } else {
338 // zPathLength = lambda0*(1.0 - exp(-tPathLength/lambda0));
339 G4double e1 = 0.0;
340 if(currentRange > tPathLength) {
341 e1 = theManager->GetEnergy(particle,
342 currentRange-tPathLength,
343 currentCouple);
344 }
345 lambdaeff = GetLambda(0.5*(e1 + preKinEnergy));
346 zPathLength = lambdaeff*(1.0 - exp(-tPathLength/lambdaeff));
347 }
348 }
349 //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
350 return zPathLength;
351}
352
353//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
354
355G4double G4WentzelVIModel::ComputeTrueStepLength(G4double geomStepLength)
356{
357 // step defined other than transportation
358 if(geomStepLength == zPathLength) return tPathLength;
359
360 // step defined by transportation
361 tPathLength = geomStepLength;
362 zPathLength = geomStepLength;
363 G4double tau = zPathLength/lambdaeff;
364 tPathLength *= (1.0 + 0.5*tau + tau*tau/3.0);
365
366 if(tau > numlimit) {
367 G4double e1 = 0.0;
368 if(currentRange > tPathLength) {
369 e1 = theManager->GetEnergy(particle,
370 currentRange-tPathLength,
371 currentCouple);
372 }
373 lambdaeff = GetLambda(0.5*(e1 + preKinEnergy));
374 tau = zPathLength/lambdaeff;
375
376 if(tau < 0.999999) tPathLength = -lambdaeff*log(1.0 - tau);
377 else tPathLength = currentRange;
378
379 if(tPathLength < zPathLength) tPathLength = zPathLength;
380 }
381 if(tPathLength > currentRange) tPathLength = currentRange;
382 //G4cout<<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
383 return tPathLength;
384}
385
386//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
387
388void G4WentzelVIModel::SampleScattering(const G4DynamicParticle* dynParticle,
389 G4double safety)
390{
391 //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
392 // << particle->GetParticleName() << G4endl;
393 G4double kinEnergy = dynParticle->GetKineticEnergy();
394
395 // ignore scattering for zero step length and enegy below the limit
396 if(kinEnergy < lowEnergyLimit || tPathLength <= DBL_MIN) return;
397
398 G4double ekin = preKinEnergy;
399 if(ekin - kinEnergy > ekin*dtrl) {
400 ekin = 0.5*(preKinEnergy + kinEnergy);
401 lambdaeff = GetLambda(ekin);
402 }
403
404 G4double x1 = 0.5*tPathLength/lambdaeff;
405 G4double cut= (*currentCuts)[currentMaterialIndex];
406 /*
407 G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy<<" Eeff(MeV)= "<<ekin/MeV
408 << " L0= " << lambda0 << " Leff= " << lambdaeff
409 << " x1= " << x1 << " safety= " << safety << G4endl;
410 */
411
412 G4double xsec = 0.0;
413 G4bool largeAng = false;
414
415 // large scattering angle case
416 if(x1 > 0.5) {
417 x1 *= 0.5;
418 largeAng = true;
419
420 // normal case
421 } else {
422
423 // define threshold angle between single and multiple scattering
424 cosThetaMin = 1.0 - 3.0*x1;
425
426 // for low-energy e-,e+ no limit
427 SetupKinematic(ekin, cut);
428
429 // recompute transport cross section
430 if(cosThetaMin > cosTetMaxNuc) {
431
432 xsec = ComputeXSectionPerVolume();
433
434 if(xtsec > 0.0) x1 = 0.5*tPathLength*xtsec;
435 else x1 = 0.0;
436
437 /*
438 G4cout << "cosTetMaxNuc= " << cosTetMaxNuc
439 << " cosThetaMin= " << cosThetaMin
440 << " cosThetaMax= " << cosThetaMax
441 << " cosTetMaxElec2= " << cosTetMaxElec2 << G4endl;
442 G4cout << "Recomputed xsec(1/mm)= " << xsec << " x1= " << x1 << G4endl;
443 */
444 }
445 }
446
447 // result of central part sampling
448 G4double z;
449 do {
450 z = -x1*log(G4UniformRand());
451 } while (z > 1.0);
452
453 // cost is sampled ------------------------------
454 G4double cost = 1.0 - 2.0*z;
455 // if(cost < -1.0) cost = -1.0;
456 // else if(cost > 1.0) cost = 1.0;
457 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
458
459 G4double phi = twopi*G4UniformRand();
460
461 G4double dirx = sint*cos(phi);
462 G4double diry = sint*sin(phi);
463
464 //G4cout << "G4WentzelVIModel: step(mm)= " << tPathLength/mm
465 // << " sint= " << sint << " cost= " << cost<< G4endl;
466
467 G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
468 G4ThreeVector newDirection(dirx,diry,cost);
469 G4ThreeVector temp(0.0,0.0,1.0);
470 G4ThreeVector pos(0.0,0.0,-zPathLength);
471 G4ThreeVector dir(0.0,0.0,1.0);
472 G4bool isscat = false;
473
474 // sample MSC scattering for large angle
475 // extra central scattering for half step
476 if(largeAng) {
477 isscat = true;
478 pos.setZ(-0.5*zPathLength);
479 do {
480 z = -x1*log(G4UniformRand());
481 } while (z > 1.0);
482 cost = 1.0 - 2.0*z;
483 //if(std::fabs(cost) > 1.0) cost = 1.0;
484
485 sint = sqrt((1.0 - cost)*(1.0 + cost));
486 phi = twopi*G4UniformRand();
487
488 // position and direction for secondary scattering
489 dir.set(sint*cos(phi),sint*sin(phi),cost);
490 pos += 0.5*dir*zPathLength;
491 x1 *= 2.0;
492 }
493
494 // sample Rutherford scattering for large angle
495 if(xsec > DBL_MIN) {
496 G4double t = tPathLength;
497 G4int nelm = currentMaterial->GetNumberOfElements();
498 const G4ElementVector* theElementVector =
499 currentMaterial->GetElementVector();
500 do{
501 G4double x = -log(G4UniformRand())/xsec;
502 pos += dir*(zPathLength*std::min(x,t)/tPathLength);
503 t -= x;
504 if(t > 0.0) {
505 G4double zz1 = 1.0;
506 G4double qsec = G4UniformRand()*xsec;
507
508 // scattering off nucleus
509 G4int i = 0;
510 if(nelm > 1) {
511 for (; i<nelm; i++) {if(xsecn[i] >= qsec) break;}
512 if(i >= nelm) i = nelm - 1;
513 }
514 SetupTarget((*theElementVector)[i]->GetZ(), tkin);
515 G4double formf = formfactA;
516 G4double costm = cosTetMaxNuc2;
517 if(prob[i] > 0.0) {
518 if(G4UniformRand() <= prob[i]) {
519 formf = 0.0;
520 costm = cosTetMaxElec2;
521 }
522 }
523 if(cosThetaMin > costm) {
524
525 G4double w1 = 1. - cosThetaMin + screenZ;
526 G4double w2 = 1. - costm + screenZ;
527 G4double w3 = cosThetaMin - costm;
528 G4double grej, zz;
529 do {
530 zz = w1*w2/(w1 + G4UniformRand()*w3) - screenZ;
531 grej = 1.0/(1.0 + formf*zz);
532 } while ( G4UniformRand() > grej*grej );
533 if(zz < 0.0) zz = 0.0;
534 else if(zz > 2.0) zz = 2.0;
535 zz1 = 1.0 - zz;
536 }
537 if(zz1 < 1.0) {
538 isscat = true;
539 //G4cout << "Rutherford zz1= " << zz1 << " t= " << t << G4endl;
540 sint = sqrt((1.0 - zz1)*(1.0 + zz1));
541 //G4cout << "sint= " << sint << G4endl;
542 phi = twopi*G4UniformRand();
543 G4double vx1 = sint*cos(phi);
544 G4double vy1 = sint*sin(phi);
545 temp.set(vx1,vy1,zz1);
546 temp.rotateUz(dir);
547 dir = temp;
548 }
549 }
550 } while (t > 0.0);
551 }
552 if(isscat) newDirection.rotateUz(dir);
553 newDirection.rotateUz(oldDirection);
554
555 //G4cout << "G4WentzelVIModel sampling of scattering is done" << G4endl;
556 // end of sampling -------------------------------
557
558 fParticleChange->ProposeMomentumDirection(newDirection);
559
560 if (latDisplasment && safety > tlimitminfix) {
561 G4double rms = invsqrt12*sqrt(2.0*x1);
562 G4double dx = zPathLength*(0.5*dirx + rms*G4RandGauss::shoot(0.0,1.0));
563 G4double dy = zPathLength*(0.5*diry + rms*G4RandGauss::shoot(0.0,1.0));
564 G4double dz;
565 G4double d = (dx*dx + dy*dy)/(zPathLength*zPathLength);
566 if(d < numlimit) dz = -0.5*zPathLength*d*(1.0 + 0.25*d);
567 else if(d < 1.0) dz = -zPathLength*(1.0 - sqrt(1.0 - d));
568 else {
569 dx = dy = dz = 0.0;
570 }
571
572 temp.set(dx,dy,dz);
573 if(isscat) temp.rotateUz(dir);
574 pos += temp;
575
576 pos.rotateUz(oldDirection);
577
578 G4double r = pos.mag();
579
580 /*
581 G4cout << " r(mm)= " << r << " safety= " << safety
582 << " trueStep(mm)= " << tPathLength
583 << " geomStep(mm)= " << zPathLength
584 << G4endl;
585 */
586
587 if(r > tlimitminfix) {
588 pos /= r;
589 ComputeDisplacement(fParticleChange, pos, r, safety);
590 }
591 }
592 //G4cout << "G4WentzelVIModel::SampleScattering end" << G4endl;
593}
594
595//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
596
597G4double G4WentzelVIModel::ComputeXSectionPerVolume()
598{
599 const G4ElementVector* theElementVector =
600 currentMaterial->GetElementVector();
601 const G4double* theAtomNumDensityVector =
602 currentMaterial->GetVecNbOfAtomsPerVolume();
603 G4int nelm = currentMaterial->GetNumberOfElements();
604 if(nelm > nelments) {
605 nelments = nelm;
606 xsecn.resize(nelments);
607 prob.resize(nelments);
608 }
609
610 xtsec = 0.0;
611 G4double xs = 0.0;
612
613 for (G4int i=0; i<nelm; i++) {
614 SetupTarget((*theElementVector)[i]->GetZ(), tkin);
615 G4double density = theAtomNumDensityVector[i];
616 G4double cosnm = cosTetMaxNuc2;
617 G4double cosem = cosTetMaxElec2;
618
619 // recompute the angular limit
620 cosTetMaxNuc2 = std::max(cosnm,cosThetaMin);
621 cosTetMaxElec2 = std::max(cosem,cosThetaMin);
622 xtsec += ComputeTransportXSectionPerAtom()*density;
623 // return limit back
624 cosTetMaxElec2 = cosem;
625 cosTetMaxNuc2 = cosnm;
626
627 G4double esec = 0.0;
628 G4double nsec = 0.0;
629 G4double x1 = 1.0 - cosThetaMin + screenZ;
630 G4double f = kinFactor*density;
631
632 // scattering off electrons
633 if(cosThetaMin > cosem) {
634 esec = f*(cosThetaMin - cosem)/(x1*(1.0 - cosem + screenZ));
635 }
636
637 // scattering off nucleaus
638 if(cosThetaMin > cosnm) {
639
640 // Rutherford part
641 G4double s = screenZ*formfactA;
642 G4double z1 = 1.0 - cosnm + screenZ;
643 G4double s1 = 1.0 - s;
644 G4double d = s1/formfactA;
645
646 // check numerical limit
647 if(d < numlimit*x1) {
648 G4double x2 = x1*x1;
649 G4double z2 = z1*z1;
650 nsec = (1.0/(x1*x2) - 1.0/(z1*z2) - d*1.5*(1.0/(x2*x2) - 1.0/(z2*z2)))/
651 (3.0*formfactA*formfactA);
652 } else {
653 G4double x2 = x1 + d;
654 G4double z2 = z1 + d;
655 nsec = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1);
656 }
657 nsec *= f*targetZ;
658 }
659 nsec += esec;
660 if(nsec > 0.0) esec /= nsec;
661 xs += nsec;
662 xsecn[i] = xs;
663 prob[i] = esec;
664 //G4cout << i << " xs= " << xs << " cosThetaMin= " << cosThetaMin
665 // << " costm= " << costm << G4endl;
666 }
667
668 //G4cout << "ComputeXS result: xsec(1/mm)= " << xs
669 //<< " txsec(1/mm)= " << xtsec <<G4endl;
670 return xs;
671}
672
673//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
674
675/*
676G4double G4MuMscModel::ComputeXSectionPerVolume()
677{
678 const G4ElementVector* theElementVector =
679 currentMaterial->GetElementVector();
680 const G4double* theAtomNumDensityVector =
681 currentMaterial->GetVecNbOfAtomsPerVolume();
682 size_t nelm = currentMaterial->GetNumberOfElements();
683
684 xsece1 = 0.0;
685 xsece2 = 0.0;
686 xsecn2 = 0.0;
687 zcorr = 0.0;
688
689 G4double fac = coeff*chargeSquare*invbeta2/mom2;
690
691 for (size_t i=0; i<nelm; i++) {
692 const G4Element* elm = (*theElementVector)[i];
693 G4double Z = elm->GetZ();
694 SetupTarget(Z, tkin);
695 G4double den = fac*theAtomNumDensityVector[i]*Z;
696
697 G4double x = 1.0 - cosThetaMin;
698 G4double x1 = x + screenZ;
699 G4double x2 = 1.0/(x1*x1);
700 G4double x3 = 1.0 + x*formfactA;
701
702 //G4cout << "x= " << x << " den= " << den << " cosE= " << cosTetMaxElec << G4endl;
703 //G4cout << "cosThtaMin= " << cosThetaMin << G4endl;
704 //G4cout << "cosTetMaxNuc= " << cosTetMaxNuc << " q2Limit= " << q2Limit << G4endl;
705
706 // scattering off electrons
707 if(cosTetMaxElec < cosThetaMin) {
708
709 // flat part
710 G4double s = den*x2*x;
711 xsece1 += s;
712 zcorr += 0.5*x*s;
713
714 // Rutherford part
715 G4double z1 = 1.0 - cosTetMaxElec + screenZ;
716 G4double z2 = (cosThetaMin - cosTetMaxElec)/x1;
717 if(z2 < 0.2) s = z2*(x - 0.5*z2*(x - screenZ))/x1;
718 else s = log(1.0 + z2) - screenZ*z2/z1;
719 xsece2 += den*z2/z1;
720 zcorr += den*s;
721 }
722 den *= Z;
723
724 //G4cout << "Z= " << Z<< " cosL= " << cosTetMaxNuc << " cosMin= " << cosThetaMin << G4endl;
725 // scattering off nucleaus
726 if(cosTetMaxNuc < cosThetaMin) {
727
728 // flat part
729 G4double s = den*x2*x/(x3*x3);
730 xsece1 += s;
731 zcorr += 0.5*x*s;
732
733 // Rutherford part
734 s = screenZ*formfactA;
735 G4double w = 1.0 + 2.0*s;
736 G4double z1 = 1.0 - cosTetMaxNuc + screenZ;
737 G4double d = (1.0 - s)/formfactA;
738 G4double x4 = x1 + d;
739 G4double z4 = z1 + d;
740 G4double t1 = 1.0/(x1*z1);
741 G4double t4 = 1.0/(x4*z4);
742 G4double w1 = cosThetaMin - cosTetMaxNuc;
743 G4double w2 = log(z1*x4/(x1*z4));
744
745 den *= w;
746 xsecn2 += den*(w1*(t1 + t4) - 2.0*w2/d);
747 zcorr += den*(w*w2 - w1*(screenZ*t1 + t4/formfactA));
748 }
749 xsece[i] = xsece2;
750 xsecn[i] = xsecn2;
751 // G4cout << i << " xsece2= " << xsece2 << " xsecn2= " << xsecn2 << G4endl;
752 }
753 G4double xsec = xsece1 + xsece2 + xsecn2;
754
755 //G4cout << "xsece1= " << xsece1 << " xsece2= " << xsece2
756 //<< " xsecn2= " << xsecn2
757 // << " zsec= " << zcorr*0.5*tPathLength << G4endl;
758 zcorr *= 0.5*tPathLength;
759
760 return xsec;
761}
762*/
763
764//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
765
766void G4WentzelVIModel::ComputeMaxElectronScattering(G4double cutEnergy)
767{
768 ecut = cutEnergy;
769 G4double tmax = tkin;
770 cosTetMaxElec = 1.0;
771 if(mass > MeV) {
772 G4double ratio = electron_mass_c2/mass;
773 G4double tau = tkin/mass;
774 tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
775 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
776 cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
777 } else {
778
779 if(particle == theElectron) tmax *= 0.5;
780 G4double t = std::min(cutEnergy, tmax);
781 G4double mom21 = t*(t + 2.0*electron_mass_c2);
782 G4double t1 = tkin - t;
783 //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
784 //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
785 if(t1 > 0.0) {
786 G4double mom22 = t1*(t1 + 2.0*mass);
787 G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
788 if(ctm < 1.0) cosTetMaxElec = ctm;
789 if(ctm < -1.0) cosTetMaxElec = -1.0;
790 }
791 }
792 if(cosTetMaxElec < cosTetMaxNuc) cosTetMaxElec = cosTetMaxNuc;
793}
794
795//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.