| 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 /// \file electromagnetic/TestEm7/src/G4Screen
27 /// \brief Implementation of the G4ScreenedNuc
28 //
29 //
30 //
31 // Class Description
32 // Process for screened electromagnetic nuclea
33 // Physics comes from:
34 // Marcus H. Mendenhall and Robert A. Weller,
35 // "Algorithms for the rapid computation o
36 // sections for screened Coulomb collision
37 // Nuclear Instruments and Methods in Phy
38 // The only input required is a screening func
39 // of the actual interatomic potential for two
40 // numbers Z1 and Z2,
41 // to the unscreened potential Z1*Z2*e^2/r whe
42 // Geant4 units
43 //
44 // First version, April 2004, Marcus H. Menden
45 //
46 // 5 May, 2004, Marcus Mendenhall
47 // Added an option for enhancing hard collisio
48 // backscattering calculations to be carried o
49 // without distorting the multiple-scattering
50 // the method SetCrossSectionHardening(G4doubl
51 // Hardeni
52 // sets what fraction of the events will be ra
53 // and the factor by which the impact area is
54 //
55 // 21 November, 2004, Marcus Mendenhall
56 // added static_nucleus to IsApplicable
57 //
58 // 7 December, 2004, Marcus Mendenhall
59 // changed mean free path of stopping particle
60 // to avoid new verbose warning about 0 MFP in
61 //
62 // 17 December, 2004, Marcus Mendenhall
63 // added code to permit screening out overly c
64 // expected to be hadronic, not Coulombic
65 //
66 // 19 December, 2004, Marcus Mendenhall
67 // massive rewrite to add modular physics stag
68 // computation. This allows one to select (e.
69 // python process and an embedded python inter
70 // for generating the tables.
71 // It also allows one to switch between sub-sa
72 // and normal scattering, and between non-rela
73 // relativistic kinematic approximations, with
74 // combination. Further, one can add extra sta
75 // implement various book-keeping processes.
76 //
77 // January 2007, Marcus Mendenhall
78 // Reorganized heavily for inclusion in Geant4
79 // one source and header, all historic code re
80 //
81 // Class Description - End
82
83 //....oooOO0OOooo........oooOO0OOooo........oo
84
85 #include "G4ScreenedNuclearRecoil.hh"
86
87 #include "globals.hh"
88
89 #include <stdio.h>
90
91 const char* G4ScreenedCoulombCrossSectionInfo:
92 {
93 return "G4ScreenedNuclearRecoil.cc,v 1.57 20
94 }
95
96 #include "c2_factory.hh"
97
98 #include "G4DataVector.hh"
99 #include "G4DynamicParticle.hh"
100 #include "G4Element.hh"
101 #include "G4ElementVector.hh"
102 #include "G4EmProcessSubType.hh"
103 #include "G4IonTable.hh"
104 #include "G4Isotope.hh"
105 #include "G4IsotopeVector.hh"
106 #include "G4LindhardPartition.hh"
107 #include "G4Material.hh"
108 #include "G4MaterialCutsCouple.hh"
109 #include "G4ParticleChangeForLoss.hh"
110 #include "G4ParticleDefinition.hh"
111 #include "G4ParticleTable.hh"
112 #include "G4ParticleTypes.hh"
113 #include "G4ProcessManager.hh"
114 #include "G4StableIsotopes.hh"
115 #include "G4Step.hh"
116 #include "G4Track.hh"
117 #include "G4VParticleChange.hh"
118 #include "Randomize.hh"
119
120 #include <iomanip>
121 #include <iostream>
122 static c2_factory<G4double> c2; // this makes
123 typedef c2_ptr<G4double> c2p;
124
125 G4ScreenedCoulombCrossSection::~G4ScreenedCoul
126 {
127 screeningData.clear();
128 MFPTables.clear();
129 }
130
131 const G4double G4ScreenedCoulombCrossSection::
132 0, 1.007940, 4.002602, 6.941000
133 15.999400, 18.998403, 20.179700, 22.98977
134 32.065000, 35.453000, 39.948000, 39.09830
135 51.996100, 54.938049, 55.845000, 58.93320
136 72.640000, 74.921600, 78.960000, 79.90400
137 91.224000, 92.906380, 95.940000, 98.00000
138 112.411000, 114.818000, 118.710000, 121.7600
139 137.327000, 138.905500, 140.116000, 140.9076
140 157.250000, 158.925340, 162.500000, 164.9303
141 178.490000, 180.947900, 183.840000, 186.2070
142 200.590000, 204.383300, 207.200000, 208.9803
143 226.000000, 227.000000, 232.038100, 231.0358
144 247.000000, 247.000000, 251.000000, 252.0000
145 261.000000, 262.000000, 266.000000, 264.0000
146 285.000000, 282.500000, 289.000000, 287.5000
147
148 G4ParticleDefinition*
149 G4ScreenedCoulombCrossSection::SelectRandomUnw
150 {
151 // Select randomly an element within the mat
152 // density only
153 const G4Material* material = couple->GetMate
154 G4int nMatElements = material->GetNumberOfEl
155 const G4ElementVector* elementVector = mater
156 const G4Element* element = 0;
157 G4ParticleDefinition* target = 0;
158
159 // Special case: the material consists of on
160 if (nMatElements == 1) {
161 element = (*elementVector)[0];
162 }
163 else {
164 // Composite material
165 G4double random = G4UniformRand() * materi
166 G4double nsum = 0.0;
167 const G4double* atomDensities = material->
168
169 for (G4int k = 0; k < nMatElements; k++) {
170 nsum += atomDensities[k];
171 element = (*elementVector)[k];
172 if (nsum >= random) break;
173 }
174 }
175
176 G4int N = 0;
177 G4int Z = element->GetZasInt();
178
179 G4int nIsotopes = element->GetNumberOfIsotop
180 if (0 < nIsotopes) {
181 if (Z <= 92) {
182 // we have no detailed material isotopic
183 // so use G4StableIsotopes table up to Z
184 static G4StableIsotopes theIso;
185 // get a stable isotope table for defaul
186 nIsotopes = theIso.GetNumberOfIsotopes(Z
187 G4double random = 100.0 * G4UniformRand(
188 // values are expressed as percent, sum
189 G4int tablestart = theIso.GetFirstIsotop
190 G4double asum = 0.0;
191 for (G4int i = 0; i < nIsotopes; i++) {
192 asum += theIso.GetAbundance(i + tables
193 N = theIso.GetIsotopeNucleonCount(i +
194 if (asum >= random) break;
195 }
196 }
197 else {
198 // too heavy for stable isotope table, j
199 N = (G4int)std::floor(element->GetN() +
200 }
201 }
202 else {
203 G4int i;
204 const G4IsotopeVector* isoV = element->Get
205 G4double random = G4UniformRand();
206 G4double* abundance = element->GetRelative
207 G4double asum = 0.0;
208 for (i = 0; i < nIsotopes; i++) {
209 asum += abundance[i];
210 N = (*isoV)[i]->GetN();
211 if (asum >= random) break;
212 }
213 }
214
215 // get the official definition of this nucle
216 // value of A note that GetIon is very slow,
217 // we have already found ourselves.
218 ParticleCache::iterator p = targetMap.find(Z
219 if (p != targetMap.end()) {
220 target = (*p).second;
221 }
222 else {
223 target = G4IonTable::GetIonTable()->GetIon
224 targetMap[Z * 1000 + N] = target;
225 }
226 return target;
227 }
228
229 void G4ScreenedCoulombCrossSection::BuildMFPTa
230 {
231 const G4int nmfpvals = 200;
232
233 std::vector<G4double> evals(nmfpvals), mfpva
234
235 // sum up inverse MFPs per element for each
236 const G4MaterialTable* materialTable = G4Mat
237 if (materialTable == 0) {
238 return;
239 }
240 // G4Exception("G4ScreenedCoulombCrossSectio
241 //- no MaterialTable found)");
242
243 G4int nMaterials = G4Material::GetNumberOfMa
244
245 for (G4int matidx = 0; matidx < nMaterials;
246 const G4Material* material = (*materialTab
247 const G4ElementVector& elementVector = *(m
248 const G4int nMatElements = material->GetNu
249
250 const G4Element* element = 0;
251 const G4double* atomDensities = material->
252
253 G4double emin = 0, emax = 0;
254 // find innermost range of cross section f
255 for (G4int kel = 0; kel < nMatElements; ke
256 element = elementVector[kel];
257 G4int Z = (G4int)std::floor(element->Get
258 const G4_c2_function& ifunc = sigmaMap[Z
259 if (!kel || ifunc.xmin() > emin) emin =
260 if (!kel || ifunc.xmax() < emax) emax =
261 }
262
263 G4double logint = std::log(emax / emin) /
264 // logarithmic increment for tables
265
266 // compute energy scale for interpolator.
267 // both ends to avoid range errors
268 for (G4int i = 1; i < nmfpvals - 1; i++)
269 evals[i] = emin * std::exp(logint * i);
270 evals.front() = emin;
271 evals.back() = emax;
272
273 // zero out the inverse mfp sums to start
274 for (G4int eidx = 0; eidx < nmfpvals; eidx
275 mfpvals[eidx] = 0.0;
276
277 // sum inverse mfp for each element in thi
278 // energy
279 for (G4int kel = 0; kel < nMatElements; ke
280 element = elementVector[kel];
281 G4int Z = (G4int)std::floor(element->Get
282 const G4_c2_function& sigma = sigmaMap[Z
283 G4double ndens = atomDensities[kel];
284 // compute atom fraction for this elemen
285
286 for (G4int eidx = 0; eidx < nmfpvals; ei
287 mfpvals[eidx] += ndens * sigma(evals[e
288 }
289 }
290
291 // convert inverse mfp to regular mfp
292 for (G4int eidx = 0; eidx < nmfpvals; eidx
293 mfpvals[eidx] = 1.0 / mfpvals[eidx];
294 }
295 // and make a new interpolating function o
296 MFPTables[matidx] = c2.log_log_interpolati
297 }
298 }
299
300 G4ScreenedNuclearRecoil::G4ScreenedNuclearReco
301
302
303
304 : G4VDiscreteProcess(processName, fElectroma
305 screeningKey(ScreeningKey),
306 generateRecoils(GenerateRecoils),
307 avoidReactions(1),
308 recoilCutoff(RecoilCutoff),
309 physicsCutoff(PhysicsCutoff),
310 hardeningFraction(0.0),
311 hardeningFactor(1.0),
312 externalCrossSectionConstructor(0),
313 NIELPartitionFunction(new G4LindhardRobins
314 {
315 // for now, point to class instance of this.
316 // one fails
317 // to correctly update NIEL
318 // not even this is needed... done in G4VPro
319 // pParticleChange=&aParticleChange;
320 processMaxEnergy = 50000.0 * MeV;
321 highEnergyLimit = 100.0 * MeV;
322 lowEnergyLimit = physicsCutoff;
323 registerDepositedEnergy = 1; // by default,
324 MFPScale = 1.0;
325 // SetVerboseLevel(2);
326 AddStage(new G4ScreenedCoulombClassicalKinem
327 AddStage(new G4SingleScatter);
328 SetProcessSubType(fCoulombScattering);
329 }
330
331 void G4ScreenedNuclearRecoil::ResetTables()
332 {
333 std::map<G4int, G4ScreenedCoulombCrossSectio
334 for (; xt != crossSectionHandlers.end(); xt+
335 delete (*xt).second;
336 }
337 crossSectionHandlers.clear();
338 }
339
340 void G4ScreenedNuclearRecoil::ClearStages()
341 {
342 // I don't think I like deleting the process
343 // abandoned
344 // if the creator doesn't get rid of them
345 // std::vector<G4ScreenedCollisionStage *>::
346 // collisionStages.begin();
347 // for(; stage != collisionStages.end(); sta
348
349 collisionStages.clear();
350 }
351
352 void G4ScreenedNuclearRecoil::SetNIELPartition
353 {
354 if (NIELPartitionFunction) delete NIELPartit
355 NIELPartitionFunction = part;
356 }
357
358 void G4ScreenedNuclearRecoil::DepositEnergy(G4
359 G4
360 {
361 if (!NIELPartitionFunction) {
362 IonizingLoss += energy;
363 }
364 else {
365 G4double part = NIELPartitionFunction->Par
366 IonizingLoss += energy * (1 - part);
367 NIEL += energy * part;
368 }
369 }
370
371 G4ScreenedNuclearRecoil::~G4ScreenedNuclearRec
372 {
373 ResetTables();
374 }
375
376 // returns true if it appears the nuclei colli
377 // in checking
378 G4bool G4ScreenedNuclearRecoil::CheckNuclearCo
379 {
380 return avoidReactions
381 && (apsis < (1.1 * (std::pow(A, 1.0 /
382 // nuclei are within 1.4 fm (reduced pion Co
383 // other at apsis,
384 // this is hadronic, skip it
385 }
386
387 G4ScreenedCoulombCrossSection* G4ScreenedNucle
388 {
389 G4ScreenedCoulombCrossSection* xc;
390 if (!externalCrossSectionConstructor)
391 xc = new G4NativeScreenedCoulombCrossSecti
392 else
393 xc = externalCrossSectionConstructor->crea
394 xc->SetVerbosity(verboseLevel);
395 return xc;
396 }
397
398 G4double G4ScreenedNuclearRecoil::GetMeanFreeP
399
400 {
401 const G4DynamicParticle* incoming = track.Ge
402 G4double energy = incoming->GetKineticEnergy
403 G4double a1 = incoming->GetDefinition()->Get
404
405 G4double meanFreePath;
406 *cond = NotForced;
407
408 if (energy < lowEnergyLimit || energy < reco
409 *cond = Forced;
410 return 1.0 * nm;
411 /* catch and stop slow particles to collec
412 }
413 else if (energy > processMaxEnergy * a1) {
414 return DBL_MAX; // infinite mean free pat
415 }
416 else if (energy > highEnergyLimit * a1)
417 energy = highEnergyLimit * a1;
418 /* constant MFP at high energy */
419
420 G4double fz1 = incoming->GetDefinition()->Ge
421 G4int z1 = (G4int)(fz1 / eplus + 0.5);
422
423 std::map<G4int, G4ScreenedCoulombCrossSectio
424 G4ScreenedCoulombCrossSection* xs;
425
426 if (xh == crossSectionHandlers.end()) {
427 xs = crossSectionHandlers[z1] = GetNewCros
428 xs->LoadData(screeningKey, z1, a1, physics
429 xs->BuildMFPTables();
430 }
431 else
432 xs = (*xh).second;
433
434 const G4MaterialCutsCouple* materialCouple =
435 size_t materialIndex = materialCouple->GetMa
436
437 const G4_c2_function& mfp = *(*xs)[materialI
438
439 // make absolutely certain we don't get an o
440 meanFreePath = mfp(std::min(std::max(energy,
441
442 // G4cout << "MFP: " << meanFreePath << " in
443 //<< " energy " << energy << " MFPScale " <<
444
445 return meanFreePath * MFPScale;
446 }
447
448 G4VParticleChange* G4ScreenedNuclearRecoil::Po
449 {
450 validCollision = 1;
451 pParticleChange->Initialize(aTrack);
452 NIEL = 0.0; // default is no NIEL deposited
453 IonizingLoss = 0.0;
454
455 // do universal setup
456
457 const G4DynamicParticle* incidentParticle =
458 G4ParticleDefinition* baseParticle = aTrack.
459
460 G4double fz1 = baseParticle->GetPDGCharge()
461 G4int z1 = (G4int)(fz1 + 0.5);
462 G4double a1 = baseParticle->GetPDGMass() / a
463 G4double incidentEnergy = incidentParticle->
464
465 // Select randomly one element and (possibly
466 // current material.
467 const G4MaterialCutsCouple* couple = aTrack.
468
469 const G4Material* mat = couple->GetMaterial(
470
471 G4double P = 0.0; // the impact parameter o
472
473 if (incidentEnergy < GetRecoilCutoff() * a1)
474 // check energy sanity on entry
475 DepositEnergy(z1, baseParticle->GetPDGMass
476 GetParticleChange().ProposeEnergy(0.0);
477 // stop the particle and bail out
478 validCollision = 0;
479 }
480 else {
481 G4double numberDensity = mat->GetTotNbOfAt
482 G4double lattice = 0.5 / std::pow(numberDe
483 // typical lattice half-spacing
484 G4double length = GetCurrentInteractionLen
485 G4double sigopi = 1.0 / (pi * numberDensit
486 // this is sigma0/pi
487
488 // compute the impact parameter very early
489 // as too far away, little effort is waste
490 // this is the TRIM method for determining
491 // based on the flight path
492 // this gives a cumulative distribution of
493 // N(P)= 1-exp(-pi P^2 n l)
494 // which says the probability of NOT hitti
495 // sigma= pi P^2 =exp(-sigma N l)
496 // which may be reasonable
497 if (sigopi < lattice * lattice) {
498 // normal long-flight approximation
499 P = std::sqrt(-std::log(G4UniformRand())
500 }
501 else {
502 // short-flight limit
503 P = std::sqrt(G4UniformRand()) * lattice
504 }
505
506 G4double fraction = GetHardeningFraction()
507 if (fraction && G4UniformRand() < fraction
508 // pick out some events, and increase th
509 // section by reducing the impact parame
510 P /= std::sqrt(GetHardeningFactor());
511 }
512
513 // check if we are far enough away that th
514 // must be below cutoff,
515 // and leave everything alone if so, savin
516 if (P * P > sigopi) {
517 if (GetVerboseLevel() > 1)
518 printf("ScreenedNuclear impact reject:
519 P / angstrom, std::sqrt(sigopi)
520 // no collision, don't follow up with an
521 validCollision = 0;
522 }
523 }
524
525 // find out what we hit, and record it in ou
526 kinematics.targetMaterial = mat;
527 kinematics.a1 = a1;
528
529 if (validCollision) {
530 G4ScreenedCoulombCrossSection* xsect = Get
531 G4ParticleDefinition* recoilIon = xsect->S
532 kinematics.crossSection = xsect;
533 kinematics.recoilIon = recoilIon;
534 kinematics.impactParameter = P;
535 kinematics.a2 = recoilIon->GetPDGMass() /
536 }
537 else {
538 kinematics.recoilIon = 0;
539 kinematics.impactParameter = 0;
540 kinematics.a2 = 0;
541 }
542
543 std::vector<G4ScreenedCollisionStage*>::iter
544
545 for (; stage != collisionStages.end(); stage
546 (*stage)->DoCollisionStep(this, aTrack, aS
547
548 if (registerDepositedEnergy) {
549 pParticleChange->ProposeLocalEnergyDeposit
550 pParticleChange->ProposeNonIonizingEnergyD
551 // MHM G4cout << "depositing energy, total
552 //<< IonizingLoss+NIEL << " NIEL = " << NI
553 }
554
555 return G4VDiscreteProcess::PostStepDoIt(aTra
556 }
557
558 G4ScreenedCoulombClassicalKinematics::G4Screen
559 : // instantiate all the needed functions s
560 // done at run time
561 // we will be solving x^2 - x phi(x*au)/e
562 // or, for easier scaling, x'^2 - x' au p
563 // note that only the last of these gets
564 phifunc(c2.const_plugin_function()),
565 xovereps(c2.linear(0., 0., 0.)),
566 // will fill this in with the right slope
567 diff(c2.quadratic(0., 0., 0., 1.) - xovere
568 {}
569
570 G4bool G4ScreenedCoulombClassicalKinematics::D
571
572
573 {
574 G4double au = screen->au;
575 G4CoulombKinematicsInfo& kin = master->GetKi
576 G4double A = kin.a2;
577 G4double a1 = kin.a1;
578
579 G4double xx0; // first estimate of closest
580 if (eps < 5.0) {
581 G4double y = std::log(eps);
582 G4double mlrho4 = ((((3.517e-4 * y + 1.401
583 G4double rho4 = std::exp(-mlrho4); // W&M
584 G4double bb2 = 0.5 * beta * beta;
585 xx0 = std::sqrt(bb2 + std::sqrt(bb2 * bb2
586 }
587 else {
588 G4double ee = 1.0 / (2.0 * eps);
589 xx0 = ee + std::sqrt(ee * ee + beta * beta
590 if (master->CheckNuclearCollision(A, a1, x
591 // nuclei too close
592 }
593
594 // we will be solving x^2 - x phi(x*au)/eps
595 // or, for easier scaling, x'^2 - x' au phi(
596 xovereps.reset(0., 0.0, au / eps); // slope
597 phifunc.set_function(&(screen->EMphiData.get
598 // install interpolating table
599 G4double xx1, phip, phip2;
600 G4int root_error;
601 xx1 = diff->find_root(phifunc.xmin(), std::m
602 std::min(xx0 * au, phi
603 &phip, &phip2)
604 / au;
605
606 if (root_error) {
607 G4cout << "Screened Coulomb Root Finder Er
608 G4cout << "au " << au << " A " << A << " a
609 << " beta " << beta << G4endl;
610 G4cout << " xmin " << phifunc.xmin() << "
611 G4cout << " f(xmin) " << phifunc(phifunc.x
612 << phifunc(std::min(10 * xx0 * au,
613 G4cout << " xstart " << std::min(xx0 * au,
614 << beta * beta * au * au;
615 G4cout << G4endl;
616 throw c2_exception("Failed root find");
617 }
618
619 // phiprime is scaled by one factor of au be
620 // at (xx0*au),
621 G4double phiprime = phip * au;
622
623 // lambda0 is from W&M 19
624 G4double lambda0 =
625 1.0 / std::sqrt(0.5 + beta * beta / (2.0 *
626
627 // compute the 6-term Lobatto integral alpha
628 // different coefficients)
629 // this is probably completely un-needed but
630 // quality results,
631 G4double alpha = (1.0 + lambda0) / 30.0;
632 G4double xvals[] = {0.98302349, 0.84652241,
633 G4double weights[] = {0.03472124, 0.14769029
634 for (G4int k = 0; k < 4; k++) {
635 G4double x, ff;
636 x = xx1 / xvals[k];
637 ff = 1.0 / std::sqrt(1.0 - phifunc(x * au)
638 alpha += weights[k] * ff;
639 }
640
641 phifunc.unset_function();
642 // throws an exception if used without setti
643
644 G4double thetac1 = pi * beta * alpha / xx1;
645 // complement of CM scattering angle
646 G4double sintheta = std::sin(thetac1); // n
647 G4double costheta = -std::cos(thetac1); //
648 // G4double psi=std::atan2(sintheta, costhet
649 // lab scattering angle (M&T 3rd eq. 8.69)
650
651 // numerics note: because we checked above
652 // of beta which give real recoils,
653 // we don't have to look too closely for the
654 // (which would cause sin(theta)
655 // and 1-cos(theta) to both vanish and make
656 G4double zeta = std::atan2(sintheta, 1 - cos
657 // lab recoil angle (M&T 3rd eq. 8.73)
658 G4double coszeta = std::cos(zeta);
659 G4double sinzeta = std::sin(zeta);
660
661 kin.sinTheta = sintheta;
662 kin.cosTheta = costheta;
663 kin.sinZeta = sinzeta;
664 kin.cosZeta = coszeta;
665 return 1; // all OK, collision is valid
666 }
667
668 void G4ScreenedCoulombClassicalKinematics::DoC
669
670 {
671 if (!master->GetValidCollision()) return;
672
673 G4ParticleChange& aParticleChange = master->
674 G4CoulombKinematicsInfo& kin = master->GetKi
675
676 const G4DynamicParticle* incidentParticle =
677 G4ParticleDefinition* baseParticle = aTrack.
678
679 G4double incidentEnergy = incidentParticle->
680
681 // this adjustment to a1 gives the right res
682 // (constant gamma)
683 // relativistic collisions. Hard collisions
684 // Coulombic and hadronic terms interfere an
685 G4double gamma = (1.0 + incidentEnergy / bas
686 G4double a1 = kin.a1 * gamma; // relativist
687
688 G4ParticleDefinition* recoilIon = kin.recoil
689 G4double A = recoilIon->GetPDGMass() / amu_c
690 G4int Z = (G4int)((recoilIon->GetPDGCharge()
691
692 G4double Ec = incidentEnergy * (A / (A + a1)
693 // energy in CM frame (non-relativistic!)
694 const G4ScreeningTables* screen = kin.crossS
695 G4double au = screen->au; // screening leng
696
697 G4double beta = kin.impactParameter / au;
698 // dimensionless impact parameter
699 G4double eps = Ec / (screen->z1 * Z * elm_co
700 // dimensionless energy
701
702 G4bool ok = DoScreeningComputation(master, s
703 if (!ok) {
704 master->SetValidCollision(0); // flag bad
705 return; // just bail out without setting
706 }
707
708 G4double eRecoil =
709 4 * incidentEnergy * a1 * A * kin.cosZeta
710 kin.eRecoil = eRecoil;
711
712 if (incidentEnergy - eRecoil < master->GetRe
713 aParticleChange.ProposeEnergy(0.0);
714 master->DepositEnergy(int(screen->z1), a1,
715 }
716
717 if (master->GetEnableRecoils() && eRecoil >
718 kin.recoilIon = recoilIon;
719 }
720 else {
721 kin.recoilIon = 0; // this flags no recoi
722 master->DepositEnergy(Z, A, kin.targetMate
723 }
724 }
725
726 void G4SingleScatter::DoCollisionStep(G4Screen
727 const G4
728 {
729 if (!master->GetValidCollision()) return;
730
731 G4CoulombKinematicsInfo& kin = master->GetKi
732 G4ParticleChange& aParticleChange = master->
733
734 const G4DynamicParticle* incidentParticle =
735 G4double incidentEnergy = incidentParticle->
736 G4double eRecoil = kin.eRecoil;
737
738 G4double azimuth = G4UniformRand() * (2.0 *
739 G4double sa = std::sin(azimuth);
740 G4double ca = std::cos(azimuth);
741
742 G4ThreeVector recoilMomentumDirection(kin.si
743 G4ParticleMomentum incidentDirection = incid
744 recoilMomentumDirection = recoilMomentumDire
745 G4ThreeVector recoilMomentum =
746 recoilMomentumDirection * std::sqrt(2.0 *
747
748 if (aParticleChange.GetEnergy() != 0.0) {
749 // DoKinematics hasn't stopped it!
750 G4ThreeVector beamMomentum = incidentParti
751 aParticleChange.ProposeMomentumDirection(b
752 aParticleChange.ProposeEnergy(incidentEner
753 }
754
755 if (kin.recoilIon) {
756 G4DynamicParticle* recoil =
757 new G4DynamicParticle(kin.recoilIon, rec
758
759 aParticleChange.SetNumberOfSecondaries(1);
760 aParticleChange.AddSecondary(recoil);
761 }
762 }
763
764 G4bool G4ScreenedNuclearRecoil::IsApplicable(c
765 {
766 return aParticleType == *(G4Proton::Proton()
767 || aParticleType.GetParticleType() ==
768 }
769
770 void G4ScreenedNuclearRecoil::BuildPhysicsTabl
771 {
772 G4String nam = aParticleType.GetParticleName
773 if (nam == "GenericIon" || nam == "proton" |
774 || nam == "alpha" || nam == "He3")
775 {
776 G4cout << G4endl << GetProcessName() << ":
777 << " SubType= " << GetProcessSub
778 << " maxEnergy(MeV)= " << proces
779 }
780 }
781
782 void G4ScreenedNuclearRecoil::DumpPhysicsTable
783
784 // This used to be the file mhmScreenedNuclear
785 // it has been included here to collect this f
786 // number of packages
787
788 #include "G4DataVector.hh"
789 #include "G4Element.hh"
790 #include "G4ElementVector.hh"
791 #include "G4Isotope.hh"
792 #include "G4Material.hh"
793 #include "G4MaterialCutsCouple.hh"
794
795 #include <vector>
796
797 G4_c2_function& ZBLScreening(G4int z1, G4int z
798 {
799 static const size_t ncoef = 4;
800 static G4double scales[ncoef] = {-3.2, -0.94
801 static G4double coefs[ncoef] = {0.1818, 0.50
802
803 G4double au = 0.8854 * angstrom * 0.529 / (s
804 std::vector<G4double> r(npoints), phi(npoint
805
806 for (size_t i = 0; i < npoints; i++) {
807 G4double rr = (float)i / (float)(npoints -
808 r[i] = rr * rr * rMax;
809 // use quadratic r scale to make sampling
810 G4double sum = 0.0;
811 for (size_t j = 0; j < ncoef; j++)
812 sum += coefs[j] * std::exp(scales[j] * r
813 phi[i] = sum;
814 }
815
816 // compute the derivative at the origin for
817 G4double phiprime0 = 0.0;
818 for (size_t j = 0; j < ncoef; j++)
819 phiprime0 += scales[j] * coefs[j] * std::e
820 phiprime0 *= (1.0 / au); // put back in nat
821
822 *auval = au;
823 return c2.lin_log_interpolating_function().l
824 }
825
826 G4_c2_function& MoliereScreening(G4int z1, G4i
827 {
828 static const size_t ncoef = 3;
829 static G4double scales[ncoef] = {-6.0, -1.2,
830 static G4double coefs[ncoef] = {0.10, 0.55,
831
832 G4double au = 0.8853 * 0.529 * angstrom / st
833 std::vector<G4double> r(npoints), phi(npoint
834
835 for (size_t i = 0; i < npoints; i++) {
836 G4double rr = (float)i / (float)(npoints -
837 r[i] = rr * rr * rMax;
838 // use quadratic r scale to make sampling
839 G4double sum = 0.0;
840 for (size_t j = 0; j < ncoef; j++)
841 sum += coefs[j] * std::exp(scales[j] * r
842 phi[i] = sum;
843 }
844
845 // compute the derivative at the origin for
846 G4double phiprime0 = 0.0;
847 for (size_t j = 0; j < ncoef; j++)
848 phiprime0 += scales[j] * coefs[j] * std::e
849 phiprime0 *= (1.0 / au); // put back in nat
850
851 *auval = au;
852 return c2.lin_log_interpolating_function().l
853 }
854
855 G4_c2_function& LJScreening(G4int z1, G4int z2
856 {
857 // from Loftager, Besenbacher, Jensen & Sore
858 // PhysRev A20, 1443++, 1979
859 G4double au = 0.8853 * 0.529 * angstrom / st
860 std::vector<G4double> r(npoints), phi(npoint
861
862 for (size_t i = 0; i < npoints; i++) {
863 G4double rr = (float)i / (float)(npoints -
864 r[i] = rr * rr * rMax;
865 // use quadratic r scale to make sampling
866
867 G4double y = std::sqrt(9.67 * r[i] / au);
868 G4double ysq = y * y;
869 G4double phipoly = 1 + y + 0.3344 * ysq +
870 phi[i] = phipoly * std::exp(-y);
871 // G4cout << r[i] << " " << phi[i] << G4en
872 }
873
874 // compute the derivative at the origin for
875 G4double logphiprime0 = (9.67 / 2.0) * (2 *
876 // #avoid 0/0 on first element
877 logphiprime0 *= (1.0 / au); // #put back in
878
879 *auval = au;
880 return c2.lin_log_interpolating_function().l
881 }
882
883 G4_c2_function& LJZBLScreening(G4int z1, G4int
884 {
885 // hybrid of LJ and ZBL, uses LJ if x < 0.25
886 /// connector in between. These numbers are
887 // is very near the point where the function
888 G4double auzbl, aulj;
889
890 c2p zbl = ZBLScreening(z1, z2, npoints, rMax
891 c2p lj = LJScreening(z1, z2, npoints, rMax,
892
893 G4double au = (auzbl + aulj) * 0.5;
894 lj->set_domain(lj->xmin(), 0.25 * au);
895 zbl->set_domain(1.5 * au, zbl->xmax());
896
897 c2p conn = c2.connector_function(lj->xmax(),
898 c2_piecewise_function_p<G4double>& pw = c2.p
899 c2p keepit(pw);
900 pw.append_function(lj);
901 pw.append_function(conn);
902 pw.append_function(zbl);
903
904 *auval = au;
905 keepit.release_for_return();
906 return pw;
907 }
908
909 G4NativeScreenedCoulombCrossSection::~G4Native
910
911 G4NativeScreenedCoulombCrossSection::G4NativeS
912 {
913 AddScreeningFunction("zbl", ZBLScreening);
914 AddScreeningFunction("lj", LJScreening);
915 AddScreeningFunction("mol", MoliereScreening
916 AddScreeningFunction("ljzbl", LJZBLScreening
917 }
918
919 std::vector<G4String> G4NativeScreenedCoulombC
920 {
921 std::vector<G4String> keys;
922 // find the available screening keys
923 std::map<std::string, ScreeningFunc>::const_
924 for (; sfunciter != phiMap.end(); sfunciter+
925 keys.push_back((*sfunciter).first);
926 return keys;
927 }
928
929 static inline G4double cm_energy(G4double a1,
930 {
931 // "relativistically correct energy in CM fr
932 G4double m1 = a1 * amu_c2, mass2 = a2 * amu_
933 G4double mc2 = (m1 + mass2);
934 G4double f = 2.0 * mass2 * t0 / (mc2 * mc2);
935 // old way: return (f < 1e-6) ? 0.5*mc2*f :
936 // formally equivalent to previous, but nume
937 // f without conditional
938 // uses identity (sqrt(1+x) - 1)(sqrt(1+x) +
939 return mc2 * f / (std::sqrt(1.0 + f) + 1.0);
940 }
941
942 static inline G4double thetac(G4double m1, G4d
943 {
944 G4double s2th2 = eratio * ((m1 + mass2) * (m
945 G4double sth2 = std::sqrt(s2th2);
946 return 2.0 * std::asin(sth2);
947 }
948
949 void G4NativeScreenedCoulombCrossSection::Load
950
951 {
952 static const size_t sigLen = 200;
953 // since sigma doesn't matter much, a very c
954 G4DataVector energies(sigLen);
955 G4DataVector data(sigLen);
956
957 a1 = standardmass(z1);
958 // use standardized values for mass for buil
959
960 const G4MaterialTable* materialTable = G4Mat
961 G4int nMaterials = G4Material::GetNumberOfMa
962
963 for (G4int im = 0; im < nMaterials; im++) {
964 const G4Material* material = (*materialTab
965 const G4ElementVector* elementVector = mat
966 const G4int nMatElements = material->GetNu
967
968 for (G4int iEl = 0; iEl < nMatElements; iE
969 const G4Element* element = (*elementVect
970 G4int Z = element->GetZasInt();
971 G4double a2 = element->GetA() * (mole /
972
973 if (sigmaMap.find(Z) != sigmaMap.end())
974 // we've already got this element
975
976 // find the screening function generator
977 std::map<std::string, ScreeningFunc>::it
978 if (sfunciter == phiMap.end()) {
979 G4ExceptionDescription ed;
980 ed << "No such screening key <" << scr
981 G4Exception("G4NativeScreenedCoulombCr
982 }
983 ScreeningFunc sfunc = (*sfunciter).secon
984
985 G4double au;
986 G4_c2_ptr screen = sfunc(z1, Z, 200, 50.
987 // generate the screening data
988 G4ScreeningTables st;
989
990 st.EMphiData = screen; // save our phi
991 st.z1 = z1;
992 st.m1 = a1;
993 st.z2 = Z;
994 st.m2 = a2;
995 st.emin = recoilCutoff;
996 st.au = au;
997
998 // now comes the hard part... build the
999 // tables from the phi table
1000 // based on (pi-thetac) = pi*beta*alpha
1001 // alpha is very nearly unity, always
1002 // so just solve it wth alpha=1, which
1003 // much easier
1004 // this function returns an approximati
1005 // (beta/x0)^2=phi(x0)/(eps*x0)-1 ~ ((p
1006 // Since we don't need exact sigma valu
1007 // (within a factor of 2 almost always)
1008 // this rearranges to phi(x0)/(x0*eps)
1009 // 2*theta/pi - theta^2/pi^2
1010
1011 c2_linear_p<G4double>& c2eps = c2.linea
1012 // will store an appropriate eps inside
1013 G4_c2_ptr phiau = screen(c2.linear(0.0,
1014 G4_c2_ptr x0func(phiau / c2eps);
1015 // this will be phi(x)/(x*eps) when c2e
1016 x0func->set_domain(1e-6 * angstrom / au
1017 // needed for inverse function
1018 // use the c2_inverse_function interfac
1019 // it is more efficient for an ordered
1020 // computation of values.
1021 G4_c2_ptr x0_solution(c2.inverse_functi
1022
1023 G4double m1c2 = a1 * amu_c2;
1024 G4double escale = z1 * Z * elm_coupling
1025 // energy at screening distance
1026 G4double emax = m1c2;
1027 // model is doubtful in very relativist
1028 G4double eratkin = 0.9999 * (4 * a1 * a
1029 // #maximum kinematic ratio possible at
1030 G4double cmfact0 = st.emin / cm_energy(
1031 G4double l1 = std::log(emax);
1032 G4double l0 = std::log(st.emin * cmfact
1033
1034 if (verbosity >= 1)
1035 G4cout << "Native Screening: " << scr
1036 << a2 << " " << recoilCutoff <
1037
1038 for (size_t idx = 0; idx < sigLen; idx+
1039 G4double ee = std::exp(idx * ((l1 - l
1040 G4double gamma = 1.0 + ee / m1c2;
1041 G4double eratio = (cmfact0 * st.emin)
1042 // factor by which ee needs to be red
1043 G4double theta = thetac(gamma * a1, a
1044
1045 G4double eps = cm_energy(a1, a2, ee)
1046 // #make sure lab energy is converted
1047 // calculations
1048 c2eps.reset(0.0, 0.0, eps);
1049 // set correct slope in this function
1050
1051 G4double q = theta / pi;
1052 // G4cout << ee << " " << m1c2 << " "
1053 // << eps << " " << theta << " " << q
1054 // old way using root finder
1055 // G4double x0= x0func->find_root(1e-
1056 // 0.9999*screen.xmax()/au, 1.0, 2*q-
1057 // new way using c2_inverse_function
1058 // useful information so should be a
1059 // since we are scanning this in stri
1060 G4double x0 = 0;
1061 try {
1062 x0 = x0_solution(2 * q - q * q);
1063 }
1064 catch (c2_exception&) {
1065 G4Exception("G4ScreenedNuclearRecoi
1066 "failure in inverse sol
1067 }
1068 G4double betasquared = x0 * x0 - x0 *
1069 G4double sigma = pi * betasquared * a
1070 energies[idx] = ee;
1071 data[idx] = sigma;
1072 }
1073 screeningData[Z] = st;
1074 sigmaMap[Z] = c2.log_log_interpolating_
1075 }
1076 }
1077 }
1078