| Geant4 Cross Reference |
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 /// \file electromagnetic/TestEm7/src/G4ScreenedNuclearRecoil.cc
27 /// \brief Implementation of the G4ScreenedNuclearRecoil class
28 //
29 //
30 //
31 // Class Description
32 // Process for screened electromagnetic nuclear elastic scattering;
33 // Physics comes from:
34 // Marcus H. Mendenhall and Robert A. Weller,
35 // "Algorithms for the rapid computation of classical cross
36 // sections for screened Coulomb collisions "
37 // Nuclear Instruments and Methods in Physics Research B58 (1991) 11-17
38 // The only input required is a screening function phi(r/a) which is the ratio
39 // of the actual interatomic potential for two atoms with atomic
40 // numbers Z1 and Z2,
41 // to the unscreened potential Z1*Z2*e^2/r where e^2 is elm_coupling in
42 // Geant4 units
43 //
44 // First version, April 2004, Marcus H. Mendenhall, Vanderbilt University
45 //
46 // 5 May, 2004, Marcus Mendenhall
47 // Added an option for enhancing hard collisions statistically, to allow
48 // backscattering calculations to be carried out with much improved event rates,
49 // without distorting the multiple-scattering broadening too much.
50 // the method SetCrossSectionHardening(G4double fraction, G4double
51 // HardeningFactor)
52 // sets what fraction of the events will be randomly hardened,
53 // and the factor by which the impact area is reduced for such selected events.
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 from 0.0 to 1.0*nanometer
60 // to avoid new verbose warning about 0 MFP in 4.6.2p02
61 //
62 // 17 December, 2004, Marcus Mendenhall
63 // added code to permit screening out overly close collisions which are
64 // expected to be hadronic, not Coulombic
65 //
66 // 19 December, 2004, Marcus Mendenhall
67 // massive rewrite to add modular physics stages and plug-in cross section table
68 // computation. This allows one to select (e.g.) between the normal external
69 // python process and an embedded python interpreter (which is much faster)
70 // for generating the tables.
71 // It also allows one to switch between sub-sampled scattering (event biasing)
72 // and normal scattering, and between non-relativistic kinematics and
73 // relativistic kinematic approximations, without having a class for every
74 // combination. Further, one can add extra stages to the scattering, which can
75 // implement various book-keeping processes.
76 //
77 // January 2007, Marcus Mendenhall
78 // Reorganized heavily for inclusion in Geant4 Core. All modules merged into
79 // one source and header, all historic code removed.
80 //
81 // Class Description - End
82
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84
85 #include "G4ScreenedNuclearRecoil.hh"
86
87 #include "globals.hh"
88
89 #include <stdio.h>
90
91 const char* G4ScreenedCoulombCrossSectionInfo::CVSFileVers()
92 {
93 return "G4ScreenedNuclearRecoil.cc,v 1.57 2008/05/07 11:51:26 marcus Exp GEANT4 tag ";
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 a lot of notation shorter
123 typedef c2_ptr<G4double> c2p;
124
125 G4ScreenedCoulombCrossSection::~G4ScreenedCoulombCrossSection()
126 {
127 screeningData.clear();
128 MFPTables.clear();
129 }
130
131 const G4double G4ScreenedCoulombCrossSection::massmap[nMassMapElements + 1] = {
132 0, 1.007940, 4.002602, 6.941000, 9.012182, 10.811000, 12.010700, 14.006700,
133 15.999400, 18.998403, 20.179700, 22.989770, 24.305000, 26.981538, 28.085500, 30.973761,
134 32.065000, 35.453000, 39.948000, 39.098300, 40.078000, 44.955910, 47.867000, 50.941500,
135 51.996100, 54.938049, 55.845000, 58.933200, 58.693400, 63.546000, 65.409000, 69.723000,
136 72.640000, 74.921600, 78.960000, 79.904000, 83.798000, 85.467800, 87.620000, 88.905850,
137 91.224000, 92.906380, 95.940000, 98.000000, 101.070000, 102.905500, 106.420000, 107.868200,
138 112.411000, 114.818000, 118.710000, 121.760000, 127.600000, 126.904470, 131.293000, 132.905450,
139 137.327000, 138.905500, 140.116000, 140.907650, 144.240000, 145.000000, 150.360000, 151.964000,
140 157.250000, 158.925340, 162.500000, 164.930320, 167.259000, 168.934210, 173.040000, 174.967000,
141 178.490000, 180.947900, 183.840000, 186.207000, 190.230000, 192.217000, 195.078000, 196.966550,
142 200.590000, 204.383300, 207.200000, 208.980380, 209.000000, 210.000000, 222.000000, 223.000000,
143 226.000000, 227.000000, 232.038100, 231.035880, 238.028910, 237.000000, 244.000000, 243.000000,
144 247.000000, 247.000000, 251.000000, 252.000000, 257.000000, 258.000000, 259.000000, 262.000000,
145 261.000000, 262.000000, 266.000000, 264.000000, 277.000000, 268.000000, 281.000000, 272.000000,
146 285.000000, 282.500000, 289.000000, 287.500000, 292.000000};
147
148 G4ParticleDefinition*
149 G4ScreenedCoulombCrossSection::SelectRandomUnweightedTarget(const G4MaterialCutsCouple* couple)
150 {
151 // Select randomly an element within the material, according to number
152 // density only
153 const G4Material* material = couple->GetMaterial();
154 G4int nMatElements = material->GetNumberOfElements();
155 const G4ElementVector* elementVector = material->GetElementVector();
156 const G4Element* element = 0;
157 G4ParticleDefinition* target = 0;
158
159 // Special case: the material consists of one element
160 if (nMatElements == 1) {
161 element = (*elementVector)[0];
162 }
163 else {
164 // Composite material
165 G4double random = G4UniformRand() * material->GetTotNbOfAtomsPerVolume();
166 G4double nsum = 0.0;
167 const G4double* atomDensities = material->GetVecNbOfAtomsPerVolume();
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->GetNumberOfIsotopes();
180 if (0 < nIsotopes) {
181 if (Z <= 92) {
182 // we have no detailed material isotopic info available,
183 // so use G4StableIsotopes table up to Z=92
184 static G4StableIsotopes theIso;
185 // get a stable isotope table for default results
186 nIsotopes = theIso.GetNumberOfIsotopes(Z);
187 G4double random = 100.0 * G4UniformRand();
188 // values are expressed as percent, sum is 100
189 G4int tablestart = theIso.GetFirstIsotope(Z);
190 G4double asum = 0.0;
191 for (G4int i = 0; i < nIsotopes; i++) {
192 asum += theIso.GetAbundance(i + tablestart);
193 N = theIso.GetIsotopeNucleonCount(i + tablestart);
194 if (asum >= random) break;
195 }
196 }
197 else {
198 // too heavy for stable isotope table, just use mean mass
199 N = (G4int)std::floor(element->GetN() + 0.5);
200 }
201 }
202 else {
203 G4int i;
204 const G4IsotopeVector* isoV = element->GetIsotopeVector();
205 G4double random = G4UniformRand();
206 G4double* abundance = element->GetRelativeAbundanceVector();
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 nucleus, to get the correct
216 // value of A note that GetIon is very slow, so we will cache ones
217 // we have already found ourselves.
218 ParticleCache::iterator p = targetMap.find(Z * 1000 + N);
219 if (p != targetMap.end()) {
220 target = (*p).second;
221 }
222 else {
223 target = G4IonTable::GetIonTable()->GetIon(Z, N, 0.0);
224 targetMap[Z * 1000 + N] = target;
225 }
226 return target;
227 }
228
229 void G4ScreenedCoulombCrossSection::BuildMFPTables()
230 {
231 const G4int nmfpvals = 200;
232
233 std::vector<G4double> evals(nmfpvals), mfpvals(nmfpvals);
234
235 // sum up inverse MFPs per element for each material
236 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
237 if (materialTable == 0) {
238 return;
239 }
240 // G4Exception("G4ScreenedCoulombCrossSection::BuildMFPTables
241 //- no MaterialTable found)");
242
243 G4int nMaterials = G4Material::GetNumberOfMaterials();
244
245 for (G4int matidx = 0; matidx < nMaterials; matidx++) {
246 const G4Material* material = (*materialTable)[matidx];
247 const G4ElementVector& elementVector = *(material->GetElementVector());
248 const G4int nMatElements = material->GetNumberOfElements();
249
250 const G4Element* element = 0;
251 const G4double* atomDensities = material->GetVecNbOfAtomsPerVolume();
252
253 G4double emin = 0, emax = 0;
254 // find innermost range of cross section functions
255 for (G4int kel = 0; kel < nMatElements; kel++) {
256 element = elementVector[kel];
257 G4int Z = (G4int)std::floor(element->GetZ() + 0.5);
258 const G4_c2_function& ifunc = sigmaMap[Z];
259 if (!kel || ifunc.xmin() > emin) emin = ifunc.xmin();
260 if (!kel || ifunc.xmax() < emax) emax = ifunc.xmax();
261 }
262
263 G4double logint = std::log(emax / emin) / (nmfpvals - 1);
264 // logarithmic increment for tables
265
266 // compute energy scale for interpolator. Force exact values at
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 this material and for each
278 // energy
279 for (G4int kel = 0; kel < nMatElements; kel++) {
280 element = elementVector[kel];
281 G4int Z = (G4int)std::floor(element->GetZ() + 0.5);
282 const G4_c2_function& sigma = sigmaMap[Z];
283 G4double ndens = atomDensities[kel];
284 // compute atom fraction for this element in this material
285
286 for (G4int eidx = 0; eidx < nmfpvals; eidx++) {
287 mfpvals[eidx] += ndens * sigma(evals[eidx]);
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 out of the sum
296 MFPTables[matidx] = c2.log_log_interpolating_function().load(evals, mfpvals, true, 0, true, 0);
297 }
298 }
299
300 G4ScreenedNuclearRecoil::G4ScreenedNuclearRecoil(const G4String& processName,
301 const G4String& ScreeningKey,
302 G4bool GenerateRecoils, G4double RecoilCutoff,
303 G4double PhysicsCutoff)
304 : G4VDiscreteProcess(processName, fElectromagnetic),
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 G4LindhardRobinsonPartition)
314 {
315 // for now, point to class instance of this. Doing it by creating a new
316 // one fails
317 // to correctly update NIEL
318 // not even this is needed... done in G4VProcess().
319 // pParticleChange=&aParticleChange;
320 processMaxEnergy = 50000.0 * MeV;
321 highEnergyLimit = 100.0 * MeV;
322 lowEnergyLimit = physicsCutoff;
323 registerDepositedEnergy = 1; // by default, don't hide NIEL
324 MFPScale = 1.0;
325 // SetVerboseLevel(2);
326 AddStage(new G4ScreenedCoulombClassicalKinematics);
327 AddStage(new G4SingleScatter);
328 SetProcessSubType(fCoulombScattering);
329 }
330
331 void G4ScreenedNuclearRecoil::ResetTables()
332 {
333 std::map<G4int, G4ScreenedCoulombCrossSection*>::iterator xt = crossSectionHandlers.begin();
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 processes here... they are better
343 // abandoned
344 // if the creator doesn't get rid of them
345 // std::vector<G4ScreenedCollisionStage *>::iterator stage=
346 // collisionStages.begin();
347 // for(; stage != collisionStages.end(); stage++) delete (*stage);
348
349 collisionStages.clear();
350 }
351
352 void G4ScreenedNuclearRecoil::SetNIELPartitionFunction(const G4VNIELPartition* part)
353 {
354 if (NIELPartitionFunction) delete NIELPartitionFunction;
355 NIELPartitionFunction = part;
356 }
357
358 void G4ScreenedNuclearRecoil::DepositEnergy(G4int z1, G4double a1, const G4Material* material,
359 G4double energy)
360 {
361 if (!NIELPartitionFunction) {
362 IonizingLoss += energy;
363 }
364 else {
365 G4double part = NIELPartitionFunction->PartitionNIEL(z1, a1, material, energy);
366 IonizingLoss += energy * (1 - part);
367 NIEL += energy * part;
368 }
369 }
370
371 G4ScreenedNuclearRecoil::~G4ScreenedNuclearRecoil()
372 {
373 ResetTables();
374 }
375
376 // returns true if it appears the nuclei collided, and we are interested
377 // in checking
378 G4bool G4ScreenedNuclearRecoil::CheckNuclearCollision(G4double A, G4double a1, G4double apsis)
379 {
380 return avoidReactions
381 && (apsis < (1.1 * (std::pow(A, 1.0 / 3.0) + std::pow(a1, 1.0 / 3.0)) + 1.4) * fermi);
382 // nuclei are within 1.4 fm (reduced pion Compton wavelength) of each
383 // other at apsis,
384 // this is hadronic, skip it
385 }
386
387 G4ScreenedCoulombCrossSection* G4ScreenedNuclearRecoil::GetNewCrossSectionHandler(void)
388 {
389 G4ScreenedCoulombCrossSection* xc;
390 if (!externalCrossSectionConstructor)
391 xc = new G4NativeScreenedCoulombCrossSection;
392 else
393 xc = externalCrossSectionConstructor->create();
394 xc->SetVerbosity(verboseLevel);
395 return xc;
396 }
397
398 G4double G4ScreenedNuclearRecoil::GetMeanFreePath(const G4Track& track, G4double,
399 G4ForceCondition* cond)
400 {
401 const G4DynamicParticle* incoming = track.GetDynamicParticle();
402 G4double energy = incoming->GetKineticEnergy();
403 G4double a1 = incoming->GetDefinition()->GetPDGMass() / amu_c2;
404
405 G4double meanFreePath;
406 *cond = NotForced;
407
408 if (energy < lowEnergyLimit || energy < recoilCutoff * a1) {
409 *cond = Forced;
410 return 1.0 * nm;
411 /* catch and stop slow particles to collect their NIEL! */
412 }
413 else if (energy > processMaxEnergy * a1) {
414 return DBL_MAX; // infinite mean free path
415 }
416 else if (energy > highEnergyLimit * a1)
417 energy = highEnergyLimit * a1;
418 /* constant MFP at high energy */
419
420 G4double fz1 = incoming->GetDefinition()->GetPDGCharge();
421 G4int z1 = (G4int)(fz1 / eplus + 0.5);
422
423 std::map<G4int, G4ScreenedCoulombCrossSection*>::iterator xh = crossSectionHandlers.find(z1);
424 G4ScreenedCoulombCrossSection* xs;
425
426 if (xh == crossSectionHandlers.end()) {
427 xs = crossSectionHandlers[z1] = GetNewCrossSectionHandler();
428 xs->LoadData(screeningKey, z1, a1, physicsCutoff);
429 xs->BuildMFPTables();
430 }
431 else
432 xs = (*xh).second;
433
434 const G4MaterialCutsCouple* materialCouple = track.GetMaterialCutsCouple();
435 size_t materialIndex = materialCouple->GetMaterial()->GetIndex();
436
437 const G4_c2_function& mfp = *(*xs)[materialIndex];
438
439 // make absolutely certain we don't get an out-of-range energy
440 meanFreePath = mfp(std::min(std::max(energy, mfp.xmin()), mfp.xmax()));
441
442 // G4cout << "MFP: " << meanFreePath << " index " << materialIndex
443 //<< " energy " << energy << " MFPScale " << MFPScale << G4endl;
444
445 return meanFreePath * MFPScale;
446 }
447
448 G4VParticleChange* G4ScreenedNuclearRecoil::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
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 = aTrack.GetDynamicParticle();
458 G4ParticleDefinition* baseParticle = aTrack.GetDefinition();
459
460 G4double fz1 = baseParticle->GetPDGCharge() / eplus;
461 G4int z1 = (G4int)(fz1 + 0.5);
462 G4double a1 = baseParticle->GetPDGMass() / amu_c2;
463 G4double incidentEnergy = incidentParticle->GetKineticEnergy();
464
465 // Select randomly one element and (possibly) isotope in the
466 // current material.
467 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
468
469 const G4Material* mat = couple->GetMaterial();
470
471 G4double P = 0.0; // the impact parameter of this collision
472
473 if (incidentEnergy < GetRecoilCutoff() * a1) {
474 // check energy sanity on entry
475 DepositEnergy(z1, baseParticle->GetPDGMass() / amu_c2, mat, incidentEnergy);
476 GetParticleChange().ProposeEnergy(0.0);
477 // stop the particle and bail out
478 validCollision = 0;
479 }
480 else {
481 G4double numberDensity = mat->GetTotNbOfAtomsPerVolume();
482 G4double lattice = 0.5 / std::pow(numberDensity, 1.0 / 3.0);
483 // typical lattice half-spacing
484 G4double length = GetCurrentInteractionLength();
485 G4double sigopi = 1.0 / (pi * numberDensity * length);
486 // this is sigma0/pi
487
488 // compute the impact parameter very early, so if is rejected
489 // as too far away, little effort is wasted
490 // this is the TRIM method for determining an impact parameter
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 hitting a disk of area
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()) * sigopi);
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 the central cross
509 // section by reducing the impact parameter
510 P /= std::sqrt(GetHardeningFactor());
511 }
512
513 // check if we are far enough away that the energy transfer
514 // must be below cutoff,
515 // and leave everything alone if so, saving a lot of time.
516 if (P * P > sigopi) {
517 if (GetVerboseLevel() > 1)
518 printf("ScreenedNuclear impact reject: length=%.3f P=%.4f limit=%.4f\n", length / angstrom,
519 P / angstrom, std::sqrt(sigopi) / angstrom);
520 // no collision, don't follow up with anything
521 validCollision = 0;
522 }
523 }
524
525 // find out what we hit, and record it in our kinematics block.
526 kinematics.targetMaterial = mat;
527 kinematics.a1 = a1;
528
529 if (validCollision) {
530 G4ScreenedCoulombCrossSection* xsect = GetCrossSectionHandlers()[z1];
531 G4ParticleDefinition* recoilIon = xsect->SelectRandomUnweightedTarget(couple);
532 kinematics.crossSection = xsect;
533 kinematics.recoilIon = recoilIon;
534 kinematics.impactParameter = P;
535 kinematics.a2 = recoilIon->GetPDGMass() / amu_c2;
536 }
537 else {
538 kinematics.recoilIon = 0;
539 kinematics.impactParameter = 0;
540 kinematics.a2 = 0;
541 }
542
543 std::vector<G4ScreenedCollisionStage*>::iterator stage = collisionStages.begin();
544
545 for (; stage != collisionStages.end(); stage++)
546 (*stage)->DoCollisionStep(this, aTrack, aStep);
547
548 if (registerDepositedEnergy) {
549 pParticleChange->ProposeLocalEnergyDeposit(IonizingLoss + NIEL);
550 pParticleChange->ProposeNonIonizingEnergyDeposit(NIEL);
551 // MHM G4cout << "depositing energy, total = "
552 //<< IonizingLoss+NIEL << " NIEL = " << NIEL << G4endl;
553 }
554
555 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
556 }
557
558 G4ScreenedCoulombClassicalKinematics::G4ScreenedCoulombClassicalKinematics()
559 : // instantiate all the needed functions statically, so no allocation is
560 // done at run time
561 // we will be solving x^2 - x phi(x*au)/eps - beta^2 == 0.0
562 // or, for easier scaling, x'^2 - x' au phi(x')/eps - beta^2 au^2
563 // note that only the last of these gets deleted, since it owns the rest
564 phifunc(c2.const_plugin_function()),
565 xovereps(c2.linear(0., 0., 0.)),
566 // will fill this in with the right slope at run time
567 diff(c2.quadratic(0., 0., 0., 1.) - xovereps * phifunc)
568 {}
569
570 G4bool G4ScreenedCoulombClassicalKinematics::DoScreeningComputation(G4ScreenedNuclearRecoil* master,
571 const G4ScreeningTables* screen,
572 G4double eps, G4double beta)
573 {
574 G4double au = screen->au;
575 G4CoulombKinematicsInfo& kin = master->GetKinematics();
576 G4double A = kin.a2;
577 G4double a1 = kin.a1;
578
579 G4double xx0; // first estimate of closest approach
580 if (eps < 5.0) {
581 G4double y = std::log(eps);
582 G4double mlrho4 = ((((3.517e-4 * y + 1.401e-2) * y + 2.393e-1) * y + 2.734) * y + 2.220);
583 G4double rho4 = std::exp(-mlrho4); // W&M eq. 18
584 G4double bb2 = 0.5 * beta * beta;
585 xx0 = std::sqrt(bb2 + std::sqrt(bb2 * bb2 + rho4)); // W&M eq. 17
586 }
587 else {
588 G4double ee = 1.0 / (2.0 * eps);
589 xx0 = ee + std::sqrt(ee * ee + beta * beta); // W&M eq. 15 (Rutherford value)
590 if (master->CheckNuclearCollision(A, a1, xx0 * au)) return 0;
591 // nuclei too close
592 }
593
594 // we will be solving x^2 - x phi(x*au)/eps - beta^2 == 0.0
595 // or, for easier scaling, x'^2 - x' au phi(x')/eps - beta^2 au^2
596 xovereps.reset(0., 0.0, au / eps); // slope of x*au/eps term
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::min(10 * xx0 * au, phifunc.xmax()),
602 std::min(xx0 * au, phifunc.xmax()), beta * beta * au * au, &root_error,
603 &phip, &phip2)
604 / au;
605
606 if (root_error) {
607 G4cout << "Screened Coulomb Root Finder Error" << G4endl;
608 G4cout << "au " << au << " A " << A << " a1 " << a1 << " xx1 " << xx1 << " eps " << eps
609 << " beta " << beta << G4endl;
610 G4cout << " xmin " << phifunc.xmin() << " xmax " << std::min(10 * xx0 * au, phifunc.xmax());
611 G4cout << " f(xmin) " << phifunc(phifunc.xmin()) << " f(xmax) "
612 << phifunc(std::min(10 * xx0 * au, phifunc.xmax()));
613 G4cout << " xstart " << std::min(xx0 * au, phifunc.xmax()) << " target "
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 because phi is evaluated
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 * xx1 * xx1) - phiprime / (2.0 * eps));
626
627 // compute the 6-term Lobatto integral alpha (per W&M 21, with
628 // different coefficients)
629 // this is probably completely un-needed but gives the highest
630 // quality results,
631 G4double alpha = (1.0 + lambda0) / 30.0;
632 G4double xvals[] = {0.98302349, 0.84652241, 0.53235309, 0.18347974};
633 G4double weights[] = {0.03472124, 0.14769029, 0.23485003, 0.18602489};
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) / (x * eps) - beta * beta / (x * x));
638 alpha += weights[k] * ff;
639 }
640
641 phifunc.unset_function();
642 // throws an exception if used without setting again
643
644 G4double thetac1 = pi * beta * alpha / xx1;
645 // complement of CM scattering angle
646 G4double sintheta = std::sin(thetac1); // note sin(pi-theta)=sin(theta)
647 G4double costheta = -std::cos(thetac1); // note cos(pi-theta)=-cos(theta)
648 // G4double psi=std::atan2(sintheta, costheta+a1/A);
649 // lab scattering angle (M&T 3rd eq. 8.69)
650
651 // numerics note: because we checked above for reasonable values
652 // of beta which give real recoils,
653 // we don't have to look too closely for theta -> 0 here
654 // (which would cause sin(theta)
655 // and 1-cos(theta) to both vanish and make the atan2 ill behaved).
656 G4double zeta = std::atan2(sintheta, 1 - costheta);
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::DoCollisionStep(G4ScreenedNuclearRecoil* master,
669 const G4Track& aTrack, const G4Step&)
670 {
671 if (!master->GetValidCollision()) return;
672
673 G4ParticleChange& aParticleChange = master->GetParticleChange();
674 G4CoulombKinematicsInfo& kin = master->GetKinematics();
675
676 const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
677 G4ParticleDefinition* baseParticle = aTrack.GetDefinition();
678
679 G4double incidentEnergy = incidentParticle->GetKineticEnergy();
680
681 // this adjustment to a1 gives the right results for soft
682 // (constant gamma)
683 // relativistic collisions. Hard collisions are wrong anyway, since the
684 // Coulombic and hadronic terms interfere and cannot be added.
685 G4double gamma = (1.0 + incidentEnergy / baseParticle->GetPDGMass());
686 G4double a1 = kin.a1 * gamma; // relativistic gamma correction
687
688 G4ParticleDefinition* recoilIon = kin.recoilIon;
689 G4double A = recoilIon->GetPDGMass() / amu_c2;
690 G4int Z = (G4int)((recoilIon->GetPDGCharge() / eplus) + 0.5);
691
692 G4double Ec = incidentEnergy * (A / (A + a1));
693 // energy in CM frame (non-relativistic!)
694 const G4ScreeningTables* screen = kin.crossSection->GetScreening(Z);
695 G4double au = screen->au; // screening length
696
697 G4double beta = kin.impactParameter / au;
698 // dimensionless impact parameter
699 G4double eps = Ec / (screen->z1 * Z * elm_coupling / au);
700 // dimensionless energy
701
702 G4bool ok = DoScreeningComputation(master, screen, eps, beta);
703 if (!ok) {
704 master->SetValidCollision(0); // flag bad collision
705 return; // just bail out without setting valid flag
706 }
707
708 G4double eRecoil =
709 4 * incidentEnergy * a1 * A * kin.cosZeta * kin.cosZeta / ((a1 + A) * (a1 + A));
710 kin.eRecoil = eRecoil;
711
712 if (incidentEnergy - eRecoil < master->GetRecoilCutoff() * a1) {
713 aParticleChange.ProposeEnergy(0.0);
714 master->DepositEnergy(int(screen->z1), a1, kin.targetMaterial, incidentEnergy - eRecoil);
715 }
716
717 if (master->GetEnableRecoils() && eRecoil > master->GetRecoilCutoff() * kin.a2) {
718 kin.recoilIon = recoilIon;
719 }
720 else {
721 kin.recoilIon = 0; // this flags no recoil to be generated
722 master->DepositEnergy(Z, A, kin.targetMaterial, eRecoil);
723 }
724 }
725
726 void G4SingleScatter::DoCollisionStep(G4ScreenedNuclearRecoil* master, const G4Track& aTrack,
727 const G4Step&)
728 {
729 if (!master->GetValidCollision()) return;
730
731 G4CoulombKinematicsInfo& kin = master->GetKinematics();
732 G4ParticleChange& aParticleChange = master->GetParticleChange();
733
734 const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
735 G4double incidentEnergy = incidentParticle->GetKineticEnergy();
736 G4double eRecoil = kin.eRecoil;
737
738 G4double azimuth = G4UniformRand() * (2.0 * pi);
739 G4double sa = std::sin(azimuth);
740 G4double ca = std::cos(azimuth);
741
742 G4ThreeVector recoilMomentumDirection(kin.sinZeta * ca, kin.sinZeta * sa, kin.cosZeta);
743 G4ParticleMomentum incidentDirection = incidentParticle->GetMomentumDirection();
744 recoilMomentumDirection = recoilMomentumDirection.rotateUz(incidentDirection);
745 G4ThreeVector recoilMomentum =
746 recoilMomentumDirection * std::sqrt(2.0 * eRecoil * kin.a2 * amu_c2);
747
748 if (aParticleChange.GetEnergy() != 0.0) {
749 // DoKinematics hasn't stopped it!
750 G4ThreeVector beamMomentum = incidentParticle->GetMomentum() - recoilMomentum;
751 aParticleChange.ProposeMomentumDirection(beamMomentum.unit());
752 aParticleChange.ProposeEnergy(incidentEnergy - eRecoil);
753 }
754
755 if (kin.recoilIon) {
756 G4DynamicParticle* recoil =
757 new G4DynamicParticle(kin.recoilIon, recoilMomentumDirection, eRecoil);
758
759 aParticleChange.SetNumberOfSecondaries(1);
760 aParticleChange.AddSecondary(recoil);
761 }
762 }
763
764 G4bool G4ScreenedNuclearRecoil::IsApplicable(const G4ParticleDefinition& aParticleType)
765 {
766 return aParticleType == *(G4Proton::Proton()) || aParticleType.GetParticleType() == "nucleus"
767 || aParticleType.GetParticleType() == "static_nucleus";
768 }
769
770 void G4ScreenedNuclearRecoil::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
771 {
772 G4String nam = aParticleType.GetParticleName();
773 if (nam == "GenericIon" || nam == "proton" || nam == "deuteron" || nam == "triton"
774 || nam == "alpha" || nam == "He3")
775 {
776 G4cout << G4endl << GetProcessName() << ": for " << nam
777 << " SubType= " << GetProcessSubType()
778 << " maxEnergy(MeV)= " << processMaxEnergy / MeV << G4endl;
779 }
780 }
781
782 void G4ScreenedNuclearRecoil::DumpPhysicsTable(const G4ParticleDefinition&) {}
783
784 // This used to be the file mhmScreenedNuclearRecoil_native.cc
785 // it has been included here to collect this file into a smaller
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 z2, size_t npoints, G4double rMax, G4double* auval)
798 {
799 static const size_t ncoef = 4;
800 static G4double scales[ncoef] = {-3.2, -0.9432, -0.4028, -0.2016};
801 static G4double coefs[ncoef] = {0.1818, 0.5099, 0.2802, 0.0281};
802
803 G4double au = 0.8854 * angstrom * 0.529 / (std::pow(z1, 0.23) + std::pow(z2, 0.23));
804 std::vector<G4double> r(npoints), phi(npoints);
805
806 for (size_t i = 0; i < npoints; i++) {
807 G4double rr = (float)i / (float)(npoints - 1);
808 r[i] = rr * rr * rMax;
809 // use quadratic r scale to make sampling fine near the center
810 G4double sum = 0.0;
811 for (size_t j = 0; j < ncoef; j++)
812 sum += coefs[j] * std::exp(scales[j] * r[i] / au);
813 phi[i] = sum;
814 }
815
816 // compute the derivative at the origin for the spline
817 G4double phiprime0 = 0.0;
818 for (size_t j = 0; j < ncoef; j++)
819 phiprime0 += scales[j] * coefs[j] * std::exp(scales[j] * r[0] / au);
820 phiprime0 *= (1.0 / au); // put back in natural units;
821
822 *auval = au;
823 return c2.lin_log_interpolating_function().load(r, phi, false, phiprime0, true, 0);
824 }
825
826 G4_c2_function& MoliereScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval)
827 {
828 static const size_t ncoef = 3;
829 static G4double scales[ncoef] = {-6.0, -1.2, -0.3};
830 static G4double coefs[ncoef] = {0.10, 0.55, 0.35};
831
832 G4double au = 0.8853 * 0.529 * angstrom / std::sqrt(std::pow(z1, 0.6667) + std::pow(z2, 0.6667));
833 std::vector<G4double> r(npoints), phi(npoints);
834
835 for (size_t i = 0; i < npoints; i++) {
836 G4double rr = (float)i / (float)(npoints - 1);
837 r[i] = rr * rr * rMax;
838 // use quadratic r scale to make sampling fine near the center
839 G4double sum = 0.0;
840 for (size_t j = 0; j < ncoef; j++)
841 sum += coefs[j] * std::exp(scales[j] * r[i] / au);
842 phi[i] = sum;
843 }
844
845 // compute the derivative at the origin for the spline
846 G4double phiprime0 = 0.0;
847 for (size_t j = 0; j < ncoef; j++)
848 phiprime0 += scales[j] * coefs[j] * std::exp(scales[j] * r[0] / au);
849 phiprime0 *= (1.0 / au); // put back in natural units;
850
851 *auval = au;
852 return c2.lin_log_interpolating_function().load(r, phi, false, phiprime0, true, 0);
853 }
854
855 G4_c2_function& LJScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval)
856 {
857 // from Loftager, Besenbacher, Jensen & Sorensen
858 // PhysRev A20, 1443++, 1979
859 G4double au = 0.8853 * 0.529 * angstrom / std::sqrt(std::pow(z1, 0.6667) + std::pow(z2, 0.6667));
860 std::vector<G4double> r(npoints), phi(npoints);
861
862 for (size_t i = 0; i < npoints; i++) {
863 G4double rr = (float)i / (float)(npoints - 1);
864 r[i] = rr * rr * rMax;
865 // use quadratic r scale to make sampling fine near the center
866
867 G4double y = std::sqrt(9.67 * r[i] / au);
868 G4double ysq = y * y;
869 G4double phipoly = 1 + y + 0.3344 * ysq + 0.0485 * y * ysq + 0.002647 * ysq * ysq;
870 phi[i] = phipoly * std::exp(-y);
871 // G4cout << r[i] << " " << phi[i] << G4endl;
872 }
873
874 // compute the derivative at the origin for the spline
875 G4double logphiprime0 = (9.67 / 2.0) * (2 * 0.3344 - 1.0);
876 // #avoid 0/0 on first element
877 logphiprime0 *= (1.0 / au); // #put back in natural units
878
879 *auval = au;
880 return c2.lin_log_interpolating_function().load(r, phi, false, logphiprime0 * phi[0], true, 0);
881 }
882
883 G4_c2_function& LJZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval)
884 {
885 // hybrid of LJ and ZBL, uses LJ if x < 0.25*auniv, ZBL if x > 1.5*auniv, and
886 /// connector in between. These numbers are selected so the switchover
887 // is very near the point where the functions naturally cross.
888 G4double auzbl, aulj;
889
890 c2p zbl = ZBLScreening(z1, z2, npoints, rMax, &auzbl);
891 c2p lj = LJScreening(z1, z2, npoints, rMax, &aulj);
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(), lj, zbl->xmin(), zbl, true, 0);
898 c2_piecewise_function_p<G4double>& pw = c2.piecewise_function();
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::~G4NativeScreenedCoulombCrossSection() {}
910
911 G4NativeScreenedCoulombCrossSection::G4NativeScreenedCoulombCrossSection()
912 {
913 AddScreeningFunction("zbl", ZBLScreening);
914 AddScreeningFunction("lj", LJScreening);
915 AddScreeningFunction("mol", MoliereScreening);
916 AddScreeningFunction("ljzbl", LJZBLScreening);
917 }
918
919 std::vector<G4String> G4NativeScreenedCoulombCrossSection::GetScreeningKeys() const
920 {
921 std::vector<G4String> keys;
922 // find the available screening keys
923 std::map<std::string, ScreeningFunc>::const_iterator sfunciter = phiMap.begin();
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, G4double a2, G4double t0)
930 {
931 // "relativistically correct energy in CM frame"
932 G4double m1 = a1 * amu_c2, mass2 = a2 * amu_c2;
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 : mc2*(std::sqrt(1.0+f)-1.0);
936 // formally equivalent to previous, but numerically stable for all
937 // f without conditional
938 // uses identity (sqrt(1+x) - 1)(sqrt(1+x) + 1) = x
939 return mc2 * f / (std::sqrt(1.0 + f) + 1.0);
940 }
941
942 static inline G4double thetac(G4double m1, G4double mass2, G4double eratio)
943 {
944 G4double s2th2 = eratio * ((m1 + mass2) * (m1 + mass2) / (4.0 * m1 * mass2));
945 G4double sth2 = std::sqrt(s2th2);
946 return 2.0 * std::asin(sth2);
947 }
948
949 void G4NativeScreenedCoulombCrossSection::LoadData(G4String screeningKey, G4int z1, G4double a1,
950 G4double recoilCutoff)
951 {
952 static const size_t sigLen = 200;
953 // since sigma doesn't matter much, a very coarse table will do
954 G4DataVector energies(sigLen);
955 G4DataVector data(sigLen);
956
957 a1 = standardmass(z1);
958 // use standardized values for mass for building tables
959
960 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
961 G4int nMaterials = G4Material::GetNumberOfMaterials();
962
963 for (G4int im = 0; im < nMaterials; im++) {
964 const G4Material* material = (*materialTable)[im];
965 const G4ElementVector* elementVector = material->GetElementVector();
966 const G4int nMatElements = material->GetNumberOfElements();
967
968 for (G4int iEl = 0; iEl < nMatElements; iEl++) {
969 const G4Element* element = (*elementVector)[iEl];
970 G4int Z = element->GetZasInt();
971 G4double a2 = element->GetA() * (mole / gram);
972
973 if (sigmaMap.find(Z) != sigmaMap.end()) continue;
974 // we've already got this element
975
976 // find the screening function generator we need
977 std::map<std::string, ScreeningFunc>::iterator sfunciter = phiMap.find(screeningKey);
978 if (sfunciter == phiMap.end()) {
979 G4ExceptionDescription ed;
980 ed << "No such screening key <" << screeningKey << ">";
981 G4Exception("G4NativeScreenedCoulombCrossSection::LoadData", "em0003", FatalException, ed);
982 }
983 ScreeningFunc sfunc = (*sfunciter).second;
984
985 G4double au;
986 G4_c2_ptr screen = sfunc(z1, Z, 200, 50.0 * angstrom, &au);
987 // generate the screening data
988 G4ScreeningTables st;
989
990 st.EMphiData = screen; // save our phi table
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 total cross section
999 // tables from the phi table
1000 // based on (pi-thetac) = pi*beta*alpha/x0, but noting that
1001 // alpha is very nearly unity, always
1002 // so just solve it wth alpha=1, which makes the solution
1003 // much easier
1004 // this function returns an approximation to
1005 // (beta/x0)^2=phi(x0)/(eps*x0)-1 ~ ((pi-thetac)/pi)^2
1006 // Since we don't need exact sigma values, this is good enough
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.linear(0.0, 0.0, 1.0);
1012 // will store an appropriate eps inside this in loop
1013 G4_c2_ptr phiau = screen(c2.linear(0.0, 0.0, au));
1014 G4_c2_ptr x0func(phiau / c2eps);
1015 // this will be phi(x)/(x*eps) when c2eps is correctly set
1016 x0func->set_domain(1e-6 * angstrom / au, 0.9999 * screen->xmax() / au);
1017 // needed for inverse function
1018 // use the c2_inverse_function interface for the root finder
1019 // it is more efficient for an ordered
1020 // computation of values.
1021 G4_c2_ptr x0_solution(c2.inverse_function(x0func));
1022
1023 G4double m1c2 = a1 * amu_c2;
1024 G4double escale = z1 * Z * elm_coupling / au;
1025 // energy at screening distance
1026 G4double emax = m1c2;
1027 // model is doubtful in very relativistic range
1028 G4double eratkin = 0.9999 * (4 * a1 * a2) / ((a1 + a2) * (a1 + a2));
1029 // #maximum kinematic ratio possible at 180 degrees
1030 G4double cmfact0 = st.emin / cm_energy(a1, a2, st.emin);
1031 G4double l1 = std::log(emax);
1032 G4double l0 = std::log(st.emin * cmfact0 / eratkin);
1033
1034 if (verbosity >= 1)
1035 G4cout << "Native Screening: " << screeningKey << " " << z1 << " " << a1 << " " << Z << " "
1036 << a2 << " " << recoilCutoff << G4endl;
1037
1038 for (size_t idx = 0; idx < sigLen; idx++) {
1039 G4double ee = std::exp(idx * ((l1 - l0) / sigLen) + l0);
1040 G4double gamma = 1.0 + ee / m1c2;
1041 G4double eratio = (cmfact0 * st.emin) / ee;
1042 // factor by which ee needs to be reduced to get emin
1043 G4double theta = thetac(gamma * a1, a2, eratio);
1044
1045 G4double eps = cm_energy(a1, a2, ee) / escale;
1046 // #make sure lab energy is converted to CM for these
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 << " " << gamma << " "
1053 // << eps << " " << theta << " " << q << G4endl;
1054 // old way using root finder
1055 // G4double x0= x0func->find_root(1e-6*angstrom/au,
1056 // 0.9999*screen.xmax()/au, 1.0, 2*q-q*q);
1057 // new way using c2_inverse_function which caches
1058 // useful information so should be a bit faster
1059 // since we are scanning this in strict order.
1060 G4double x0 = 0;
1061 try {
1062 x0 = x0_solution(2 * q - q * q);
1063 }
1064 catch (c2_exception&) {
1065 G4Exception("G4ScreenedNuclearRecoil::LoadData", "em0003", FatalException,
1066 "failure in inverse solution to generate MFP tables");
1067 }
1068 G4double betasquared = x0 * x0 - x0 * phiau(x0) / eps;
1069 G4double sigma = pi * betasquared * au * au;
1070 energies[idx] = ee;
1071 data[idx] = sigma;
1072 }
1073 screeningData[Z] = st;
1074 sigmaMap[Z] = c2.log_log_interpolating_function().load(energies, data, true, 0, true, 0);
1075 }
1076 }
1077 }
1078