| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4WentzelVIModel
33 //
34 // Author: V.Ivanchenko
35 //
36 // Creation date: 09.04.2008 from G4MuMscModel
37 //
38 // Modifications:
39 // 27-05-2010 V.Ivanchenko added G4WentzelOKan
40 // compute cross sections and sam
41 //
42 //
43 // Class Description:
44 //
45 // Implementation of the model of multiple sca
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)
49 // L.Urban, CERN-OPEN-2006-077.
50
51 // -------------------------------------------
52 //
53
54 //....oooOO0OOooo........oooOO0OOooo........oo
55 //....oooOO0OOooo........oooOO0OOooo........oo
56
57 #include "G4WentzelVIModel.hh"
58 #include "G4PhysicalConstants.hh"
59 #include "G4SystemOfUnits.hh"
60 #include "Randomize.hh"
61 #include "G4ParticleChangeForMSC.hh"
62 #include "G4PhysicsTableHelper.hh"
63 #include "G4ElementVector.hh"
64 #include "G4ProductionCutsTable.hh"
65 #include "G4EmParameters.hh"
66 #include "G4Log.hh"
67 #include "G4Exp.hh"
68
69 //....oooOO0OOooo........oooOO0OOooo........oo
70
71 using namespace std;
72
73 const G4double invsqrt12 = 1./std::sqrt(12.);
74 const G4double numlimit = 0.1;
75 const G4int minNCollisions = 10;
76
77 G4WentzelVIModel::G4WentzelVIModel(G4bool comb
78 : G4VMscModel(nam),
79 singleScatteringMode(false),
80 isCombined(comb),
81 useSecondMoment(false)
82 {
83 tlimitminfix = 1.e-6*CLHEP::mm;
84 lowEnergyLimit = 1.0*CLHEP::eV;
85 SetSingleScatteringFactor(1.25);
86 wokvi = new G4WentzelOKandVIxSection(isCombi
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90
91 G4WentzelVIModel::~G4WentzelVIModel()
92 {
93 delete wokvi;
94 if(IsMaster()) {
95 delete fSecondMoments;
96 fSecondMoments = nullptr;
97 }
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oo
101
102 void G4WentzelVIModel::Initialise(const G4Part
103 const G4Data
104 {
105 // reset parameters
106 SetupParticle(p);
107 InitialiseParameters(p);
108 currentRange = 0.0;
109
110 if(isCombined) {
111 G4double tet = PolarAngleLimit();
112 if(tet <= 0.0) { cosThetaMax = 1
113 else if(tet < CLHEP::pi) { cosThetaMax = c
114 }
115 //G4cout << "G4WentzelVIModel::Initialise "
116 // << " " << this << " " << wokvi << G4end
117
118 wokvi->Initialise(p, cosThetaMax);
119 /*
120 G4cout << "G4WentzelVIModel: " << particle->
121 << " 1-cos(ThetaLimit)= " << 1 - cos
122 << " SingScatFactor= " << ssFactor
123 << G4endl;
124 */
125 currentCuts = &cuts;
126
127 // set values of some data members
128 fParticleChange = GetParticleChangeForMSC(p)
129
130 // Access to materials
131 const G4ProductionCutsTable* theCoupleTable
132 G4ProductionCutsTable::GetProductionCutsTa
133 G4int numOfCouples = (G4int)theCoupleTable->
134 nelments = 0;
135 for(G4int i=0; i<numOfCouples; ++i) {
136 G4int nelm = (G4int)theCoupleTable->GetMat
137 nelments = std::max(nelments, nelm);
138 }
139 xsecn.resize(nelments);
140 prob.resize(nelments);
141
142 // build second moment table only if transpo
143 G4PhysicsTable* table = GetCrossSectionTable
144 if(useSecondMoment && IsMaster() && nullptr
145
146 //G4cout << "### G4WentzelVIModel::Initial
147 // << table << G4endl;
148 fSecondMoments =
149 G4PhysicsTableHelper::PreparePhysicsTabl
150
151 G4bool splineFlag = true;
152 G4PhysicsVector* aVector = nullptr;
153 G4PhysicsVector* bVector = nullptr;
154 G4double emin = std::max(LowEnergyLimit(),
155 G4double emax = std::min(HighEnergyLimit()
156 if(emin < emax) {
157 std::size_t n = G4EmParameters::Instance
158 *G4lrint(std::log10(emax/emin));
159 if(n < 3) { n = 3; }
160
161 for(G4int i=0; i<numOfCouples; ++i) {
162
163 //G4cout<< "i= " << i << " Flag= " <<
164 // << G4endl;
165 if(fSecondMoments->GetFlag(i)) {
166 DefineMaterial(theCoupleTable->GetMa
167
168 delete (*fSecondMoments)[i];
169 if(nullptr == aVector) {
170 aVector = new G4PhysicsLogVector(e
171 bVector = aVector;
172 } else {
173 bVector = new G4PhysicsVector(*aVe
174 }
175 for(std::size_t j=0; j<n; ++j) {
176 G4double e = bVector->Energy(j);
177 bVector->PutValue(j, ComputeSecond
178 }
179 if(splineFlag) { bVector->FillSecond
180 (*fSecondMoments)[i] = bVector;
181 }
182 }
183 }
184 //G4cout << *fSecondMoments << G4endl;
185 }
186 }
187
188 //....oooOO0OOooo........oooOO0OOooo........oo
189
190 void G4WentzelVIModel::InitialiseLocal(const G
191 G4VEmMo
192 {
193 fSecondMoments = static_cast<G4WentzelVIMode
194 ->GetSecondMomentTable();
195 }
196
197 //....oooOO0OOooo........oooOO0OOooo........oo
198
199 void G4WentzelVIModel::DefineMaterial(const G4
200 {
201 if(cup != currentCouple) {
202 currentCouple = cup;
203 SetCurrentCouple(cup);
204 currentMaterial = cup->GetMaterial();
205 currentMaterialIndex = currentCouple->GetI
206 }
207 }
208
209 //....oooOO0OOooo........oooOO0OOooo........oo
210
211 G4double G4WentzelVIModel::ComputeCrossSection
212 const G4ParticleD
213 G4double kinEnerg
214 G4double Z, G4dou
215 G4double cutEnerg
216 {
217 G4double cross = 0.0;
218 SetupParticle(p);
219 if(kinEnergy < lowEnergyLimit) { return cros
220 if(nullptr == CurrentCouple()) {
221 G4Exception("G4WentzelVIModel::ComputeCros
222 FatalException, " G4MaterialCu
223 return 0.0;
224 }
225 DefineMaterial(CurrentCouple());
226 cosTetMaxNuc = wokvi->SetupKinematic(kinEner
227 if(cosTetMaxNuc < 1.0) {
228 G4double cut = (0.0 < fixedCut) ? fixedCu
229 G4double cost = wokvi->SetupTarget(G4lrint
230 cross = wokvi->ComputeTransportCrossSectio
231 /*
232 if(p->GetParticleName() == "e-")
233 G4cout << "G4WentzelVIModel::CS: Z= " << G
234 << " 1-cosN= " << 1 - cosTetMaxNuc
235 << " " << particle->GetParticleName
236 */
237 }
238 return cross;
239 }
240
241 //....oooOO0OOooo........oooOO0OOooo........oo
242
243 void G4WentzelVIModel::StartTracking(G4Track*
244 {
245 /*
246 G4cout << "G4WentzelVIModel::StartTracking "
247 << track->GetParticleDefinition()->GetParti
248 << " workvi: " << wokvi << G4endl;
249 */
250 SetupParticle(track->GetParticleDefinition()
251 }
252
253 //....oooOO0OOooo........oooOO0OOooo........oo
254
255 G4double G4WentzelVIModel::ComputeTruePathLeng
256 const G4Track& tr
257 G4double& current
258 {
259 G4double tlimit = currentMinimalStep;
260 const G4DynamicParticle* dp = track.GetDynam
261 const G4StepPoint* sp = track.GetStep()->Get
262 G4StepStatus stepStatus = sp->GetStepStatus(
263 singleScatteringMode = false;
264
265 //G4cout << "G4WentzelVIModel::ComputeTruePa
266 // << stepStatus << " " << track.Ge
267 // << G4endl;
268
269 // initialisation for each step, lambda may
270 preKinEnergy = dp->GetKineticEnergy();
271 effKinEnergy = preKinEnergy;
272 DefineMaterial(track.GetMaterialCutsCouple()
273 const G4double logPreKinEnergy = dp->GetLogK
274 lambdaeff = GetTransportMeanFreePath(particl
275 currentRange = GetRange(particle,preKinEnerg
276 cosTetMaxNuc = wokvi->SetupKinematic(preKinE
277
278 //G4cout << "lambdaeff= " << lambdaeff << "
279 // << " tlimit= " << tlimit << " 1-cost= " <
280
281 // extra check for abnormal situation
282 // this check needed to run MSC with eIoni a
283 if(tlimit > currentRange) { tlimit = current
284
285 // stop here if small range particle
286 if(tlimit < tlimitminfix) {
287 return ConvertTrueToGeom(tlimit, currentMi
288 }
289
290 // pre step
291 G4double presafety = sp->GetSafety();
292 // far from geometry boundary
293 if(currentRange < presafety) {
294 return ConvertTrueToGeom(tlimit, currentMi
295 }
296
297 // compute presafety again if presafety <= 0
298 // i.e. when it is needed for optimization p
299 if(stepStatus != fGeomBoundary && presafety
300 presafety = ComputeSafety(sp->GetPosition(
301 if(currentRange < presafety) {
302 return ConvertTrueToGeom(tlimit, current
303 }
304 }
305 /*
306 G4cout << "e(MeV)= " << preKinEnergy/MeV
307 << " " << particle->GetParticleName(
308 << " CurLimit(mm)= " << tlimit/mm <<"
309 << " R(mm)= " <<currentRange/mm
310 << " L0(mm^-1)= " << lambdaeff*mm
311 << G4endl;
312 */
313 // natural limit for high energy
314 G4double rlimit = std::max(facrange*currentR
315 (1.0 - cosTetMaxN
316
317 // low-energy e-
318 if(cosThetaMax > cosTetMaxNuc) {
319 rlimit = std::min(rlimit, facsafety*presaf
320 }
321
322 // cut correction
323 G4double rcut = currentCouple->GetProduction
324 //G4cout << "rcut= " << rcut << " rlimit= "
325 // << presafety << " 1-cosThetaMax= " <<1-co
326 //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc <
327 if(rcut > rlimit) { rlimit = std::min(rlimit
328
329 tlimit = std::min(tlimit, rlimit);
330 tlimit = std::max(tlimit, tlimitminfix);
331
332 // step limit in infinite media
333 tlimit = std::min(tlimit, 50*currentMaterial
334
335 //compute geomlimit and force few steps with
336 if (steppingAlgorithm == fUseDistanceToBound
337 && stepStatus == fGeomBoundary) {
338
339 G4double geomlimit = ComputeGeomLimit(trac
340 tlimit = std::min(tlimit, geomlimit/facgeo
341 }
342 /*
343 G4cout << particle->GetParticleName() << " e
344 << " L0= " << lambdaeff << " R= " <<
345 << " tlimit= " << tlimit
346 << " currentMinimalStep= " << curre
347 */
348 return ConvertTrueToGeom(tlimit, currentMini
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oo
352
353 G4double G4WentzelVIModel::ComputeGeomPathLeng
354 {
355 zPathLength = tPathLength = truelength;
356
357 // small step use only single scattering
358 cosThetaMin = 1.0;
359 ComputeTransportXSectionPerVolume(cosThetaMi
360 //G4cout << "xtsec= " << xtsec << " Nav= "
361 // << zPathLength*xtsec << G4endl;
362 if(0.0 >= lambdaeff || G4int(zPathLength*xts
363 singleScatteringMode = true;
364 lambdaeff = DBL_MAX;
365
366 } else {
367 //G4cout << "ComputeGeomPathLength: tLengt
368 // << " Leff= " << lambdaeff <<
369 // small step
370 if(tPathLength < numlimit*lambdaeff) {
371 G4double tau = tPathLength/lambdaeff;
372 zPathLength *= (1.0 - 0.5*tau + tau*tau/
373
374 // medium step
375 } else {
376 G4double e1 = 0.0;
377 if(currentRange > tPathLength) {
378 e1 = GetEnergy(particle,currentRange-t
379 }
380 effKinEnergy = 0.5*(e1 + preKinEnergy);
381 cosTetMaxNuc = wokvi->SetupKinematic(eff
382 lambdaeff = GetTransportMeanFreePath(par
383 //G4cout << " tLength= "<< tPathLength<<
384 zPathLength = lambdaeff;
385 if(tPathLength*numlimit < lambdaeff) {
386 zPathLength *= (1.0 - G4Exp(-tPathLeng
387 }
388 }
389 }
390 //G4cout << "Comp.geom: zLength= "<<zPathLen
391 // << tPathLength<< " Leff= " << lam
392 return zPathLength;
393 }
394
395 //....oooOO0OOooo........oooOO0OOooo........oo
396
397 G4double G4WentzelVIModel::ComputeTrueStepLeng
398 {
399 // initialisation of single scattering x-sec
400 /*
401 G4cout << "ComputeTrueStepLength: Step= " <<
402 << " geomL= " << zPathLength
403 << " Lambda= " << lambdaeff
404 << " 1-cosThetaMaxNuc= " << 1 - cos
405 */
406 if(singleScatteringMode) {
407 zPathLength = tPathLength = geomStepLength
408
409 } else {
410
411 // step defined by transportation
412 // change both geom and true step lengths
413 if(geomStepLength < zPathLength) {
414
415 // single scattering
416 if(G4int(geomStepLength*xtsec) < minNCol
417 zPathLength = tPathLength = geomStepLe
418 lambdaeff = DBL_MAX;
419 singleScatteringMode = true;
420
421 // multiple scattering
422 } else {
423 // small step
424 if(geomStepLength < numlimit*lambdaeff
425 G4double tau = geomStepLength/lambda
426 tPathLength = geomStepLength*(1.0 +
427
428 // energy correction for a big step
429 } else {
430 tPathLength *= geomStepLength/zPathL
431 G4double e1 = 0.0;
432 if(currentRange > tPathLength) {
433 e1 = GetEnergy(particle,currentRan
434 }
435 effKinEnergy = 0.5*(e1 + preKinEnerg
436 cosTetMaxNuc = wokvi->SetupKinematic
437 lambdaeff = GetTransportMeanFreePath
438 G4double tau = geomStepLength/lambda
439
440 if(tau < 0.999999) { tPathLength = -
441 else { tPathLength = c
442 }
443 zPathLength = geomStepLength;
444 }
445 }
446 }
447 // check of step length
448 // define threshold angle between single and
449 if(!singleScatteringMode) {
450 cosThetaMin -= ssFactor*tPathLength/lambda
451 xtsec = 0.0;
452
453 // recompute transport cross section - do
454 // anymore - cannot be applied for big ste
455 if(cosThetaMin > cosTetMaxNuc) {
456 // new computation
457 G4double cross = ComputeTransportXSectio
458 //G4cout << "%%%% cross= " << cross << "
459 // << " 1-cosTMin= " << 1.0 -
460 if(cross <= 0.0) {
461 singleScatteringMode = true;
462 tPathLength = zPathLength;
463 lambdaeff = DBL_MAX;
464 cosThetaMin = 1.0;
465 } else if(xtsec > 0.0) {
466
467 lambdaeff = 1./cross;
468 G4double tau = zPathLength*cross;
469 if(tau < numlimit) {
470 tPathLength = zPathLength*(1.0 + 0.5
471 } else if(tau < 0.999999) {
472 tPathLength = -lambdaeff*G4Log(1.0 -
473 } else {
474 tPathLength = currentRange;
475 }
476 }
477 }
478 }
479 tPathLength = std::min(tPathLength, currentR
480 /*
481 G4cout <<"Comp.true: zLength= "<<zPathLength
482 <<" Leff(mm)= "<<lambdaeff/mm<<" sig0
483 G4cout << particle->GetParticleName() << " 1
484 << " 1-cosTetMaxNuc= " << 1-cosTetMax
485 << " e(MeV)= " << preKinEnergy/MeV <<
486 << " SSmode= " << singleScatteringMod
487 */
488 return tPathLength;
489 }
490
491 //....oooOO0OOooo........oooOO0OOooo........oo
492
493 G4ThreeVector&
494 G4WentzelVIModel::SampleScattering(const G4Thr
495 G4double /*
496 {
497 fDisplacement.set(0.0,0.0,0.0);
498 //G4cout << "!##! G4WentzelVIModel::SampleSc
499 // << particle->GetParticleName() <<
500
501 // ignore scattering for zero step length an
502 if(preKinEnergy < lowEnergyLimit || tPathLen
503 { return fDisplacement; }
504
505 G4double invlambda = 0.0;
506 if(lambdaeff < DBL_MAX) { invlambda = 0.5/la
507
508 // use average kinetic energy over the step
509 G4double cut = (*currentCuts)[currentMateria
510 if(fixedCut > 0.0) { cut = fixedCut; }
511 /*
512 G4cout <<"SampleScat: E0(MeV)= "<< preKinEne
513 << " Leff= " << lambdaeff <<" sig0(
514 << " xmsc= " << tPathLength*invlamb
515 << " safety= " << safety << G4endl;
516 */
517 // step limit due msc
518 G4int nMscSteps = 1;
519 G4double x0 = tPathLength;
520 G4double z0 = x0*invlambda;
521 G4double prob2 = 0.0;
522
523 CLHEP::HepRandomEngine* rndmEngine = G4Rando
524
525 // large scattering angle case - two step ap
526 if(!singleScatteringMode) {
527 if(useSecondMoment) {
528 G4double z1 = invlambda*invlambda;
529 G4double z2 = SecondMoment(particle, cur
530 prob2 = (z2 - z1)/(1.5*z1 - z2);
531 }
532 // numerical limitation on step length for
533 if (z0 > 0.05 && z0 < 30.) {
534 x0 *= 0.5;
535 z0 *= 0.5;
536 nMscSteps = 2;
537 G4double zzz = G4Exp(-1.0/z0);
538 z0 += zzz;
539 prob2 *= (1.0 + zzz);
540 }
541 prob2 /= (1.0 + prob2);
542 }
543
544 // step limit due to single scattering
545 G4double x1 = 2*tPathLength;
546 if(0.0 < xtsec) { x1 = -G4Log(rndmEngine->fl
547
548 // no scattering case
549 if(singleScatteringMode && x1 > tPathLength)
550 { return fDisplacement; }
551
552 const G4ElementVector* theElementVector =
553 currentMaterial->GetElementVector();
554 std::size_t nelm = currentMaterial->GetNumbe
555
556 // geometry
557 G4double sint, cost, phi;
558 G4ThreeVector temp(0.0,0.0,1.0);
559
560 // current position and direction relative t
561 // because of magnetic field geometry is com
562 // end point of the step
563 G4ThreeVector dir(0.0,0.0,1.0);
564 fDisplacement.set(0.0,0.0,-zPathLength);
565
566 G4double mscfac = zPathLength/tPathLength;
567
568 // start a loop
569 G4double x2 = x0;
570 G4double step, z;
571 G4bool singleScat;
572 /*
573 G4cout << "Start of the loop x1(mm)= " <<
574 << " 1-cost1= " << 1 - cosThetaMin << " SS
575 << " xtsec= " << xtsec << " Nst= "
576 */
577 do {
578
579 //G4cout << "# x1(mm)= "<< x1<< " x2(mm)=
580 // single scattering case
581 if(singleScatteringMode && x1 > x2) {
582 fDisplacement += x2*mscfac*dir;
583 break;
584 }
585
586 // what is next single of multiple?
587 if(x1 <= x2) {
588 step = x1;
589 singleScat = true;
590 } else {
591 step = x2;
592 singleScat = false;
593 }
594
595 //G4cout << "# step(mm)= "<< step<< " sin
596
597 // new position
598 fDisplacement += step*mscfac*dir;
599
600 if(singleScat) {
601
602 // select element
603 std::size_t i = 0;
604 if(nelm > 1) {
605 G4double qsec = rndmEngine->flat()*xts
606 for (; i<nelm; ++i) { if(xsecn[i] >= q
607 }
608 G4double cosTetM =
609 wokvi->SetupTarget((*theElementVector)
610 //G4cout << "!!! " << cosThetaMin << "
611 // << prob[i] << G4endl;
612 temp = wokvi->SampleSingleScattering(cos
613
614 // direction is changed
615 temp.rotateUz(dir);
616 dir = temp;
617 //G4cout << dir << G4endl;
618
619 // new proposed step length
620 x2 -= step;
621 x1 = -G4Log(rndmEngine->flat())/xtsec;
622
623 // multiple scattering
624 } else {
625 --nMscSteps;
626 x1 -= step;
627 x2 = x0;
628
629 // sample z in interval 0 - 1
630 G4bool isFirst = true;
631 if(prob2 > 0.0 && rndmEngine->flat() < p
632 do {
633 if(isFirst) { z = -G4Log(rndmEngine->f
634 else { z = G4RandGamma::shoot(r
635 z *= z0;
636 // Loop checking, 03-Aug-2015, Vladimi
637 } while(z > 1.0);
638
639 cost = 1.0 - 2.0*z/*factCM*/;
640 if(cost > 1.0) { cost = 1.0; }
641 else if(cost < -1.0) { cost =-1.0; }
642 sint = sqrt((1.0 - cost)*(1.0 + cost));
643 phi = twopi*rndmEngine->flat();
644 G4double vx1 = sint*std::cos(phi);
645 G4double vy1 = sint*std::sin(phi);
646
647 // lateral displacement
648 if (latDisplasment) {
649 G4double rms = z0 > 0.0 ? invsqrt12*st
650 G4double r = x0*mscfac;
651 G4double dx = r*(0.5*vx1 + rms*G4RandG
652 G4double dy = r*(0.5*vy1 + rms*G4RandG
653 G4double d = r*r - dx*dx - dy*dy;
654
655 // change position
656 if(d >= 0.0) {
657 temp.set(dx, dy, std::sqrt(d) - r);
658 temp.rotateUz(dir);
659 fDisplacement += temp;
660 }
661 }
662 // change direction
663 temp.set(vx1,vy1,cost);
664 temp.rotateUz(dir);
665 dir = temp;
666 }
667 // Loop checking, 03-Aug-2015, Vladimir Iv
668 } while (0 < nMscSteps);
669
670 dir.rotateUz(oldDirection);
671
672 //G4cout<<"G4WentzelVIModel sampling is done
673 // end of sampling -------------------------
674
675 fParticleChange->ProposeMomentumDirection(di
676
677 // lateral displacement
678 fDisplacement.rotateUz(oldDirection);
679
680 /*
681 G4cout << " r(mm)= " << fDisplacement
682 << " safety= " << safety
683 << " trueStep(mm)= " << tPathL
684 << " geomStep(mm)= " << zPathL
685 << " x= " << fDisplacement.x()
686 << " y= " << fDisplacement.y()
687 << " z= " << fDisplacement.z()
688 << G4endl;
689 */
690
691 //G4cout<< "G4WentzelVIModel::SampleScatteri
692 return fDisplacement;
693 }
694
695 //....oooOO0OOooo........oooOO0OOooo........oo
696
697 G4double G4WentzelVIModel::ComputeTransportXSe
698 {
699 // prepare recomputation of x-sections
700 const G4ElementVector* theElementVector = cu
701 const G4double* theAtomNumDensityVector =
702 currentMaterial->GetVecNbOfAtomsPerVolume(
703 G4int nelm = (G4int)currentMaterial->GetNumb
704 if(nelm > nelments) {
705 nelments = nelm;
706 xsecn.resize(nelm);
707 prob.resize(nelm);
708 }
709
710 // check consistency
711 xtsec = 0.0;
712 if(cosTetMaxNuc >= cosTheta) { return 0.0; }
713
714 G4double cut = (*currentCuts)[currentMateria
715 if(fixedCut > 0.0) { cut = fixedCut; }
716
717 // loop over elements
718 G4double xs = 0.0;
719 for (G4int i=0; i<nelm; ++i) {
720 G4double costm =
721 wokvi->SetupTarget((*theElementVector)[i
722 G4double density = theAtomNumDensityVector
723
724 G4double esec = 0.0;
725 if(costm < cosTheta) {
726
727 // recompute the transport x-section
728 if(1.0 > cosTheta) {
729 xs += density*wokvi->ComputeTransportC
730 }
731 // recompute the total x-section
732 G4double nucsec = wokvi->ComputeNuclearC
733 esec = wokvi->ComputeElectronCrossSectio
734 nucsec += esec;
735 if(nucsec > 0.0) { esec /= nucsec; }
736 xtsec += nucsec*density;
737 }
738 xsecn[i] = xtsec;
739 prob[i] = esec;
740 //G4cout << i << " xs= " << xs << " xtsec
741 // << " 1-cosTheta= " << 1-cosTheta
742 // << " 1-cosTetMaxNuc2= " <<1-c
743 }
744
745 //G4cout << "ComputeXS result: xsec(1/mm)=
746 // << " txsec(1/mm)= " << xtsec <<G4
747 return xs;
748 }
749
750 //....oooOO0OOooo........oooOO0OOooo........oo
751
752 G4double G4WentzelVIModel::ComputeSecondMoment
753 G4double kinEnergy)
754 {
755 G4double xs = 0.0;
756
757 SetupParticle(p);
758 cosTetMaxNuc = wokvi->SetupKinematic(kinEner
759
760 if(cosTetMaxNuc >= 1.0) { return xs; }
761
762 const G4ElementVector* theElementVector = cu
763 const G4double* theAtomNumDensityVector =
764 currentMaterial->GetVecNbOfAtomsPerVolume(
765 std::size_t nelm = currentMaterial->GetNumbe
766
767 G4double cut = (*currentCuts)[currentMateria
768 if(fixedCut > 0.0) { cut = fixedCut; }
769
770 // loop over elements
771 for (std::size_t i=0; i<nelm; ++i) {
772 G4double costm =
773 wokvi->SetupTarget((*theElementVector)[i
774 xs += theAtomNumDensityVector[i]
775 *wokvi->ComputeSecondTransportMoment(cos
776 }
777 return xs;
778 }
779
780 //....oooOO0OOooo........oooOO0OOooo........oo
781
782 void G4WentzelVIModel::SetSingleScatteringFact
783 {
784 if(val > 0.05) {
785 ssFactor = val;
786 invssFactor = 1.0/(val - 0.05);
787 }
788 }
789
790 //....oooOO0OOooo........oooOO0OOooo........oo
791