source: trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

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