ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ScreenedNuclearRecoil.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ScreenedNuclearRecoil.cc
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 //
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 <stdio.h>
86 
87 #include "globals.hh"
88 
90 
92  "G4ScreenedNuclearRecoil.cc,v 1.57 2008/05/07 11:51:26 marcus Exp GEANT4 tag ";
93 }
94 
95 #include "G4ParticleTypes.hh"
96 #include "G4ParticleTable.hh"
97 #include "G4IonTable.hh"
98 #include "G4VParticleChange.hh"
100 #include "G4DataVector.hh"
101 #include "G4Track.hh"
102 #include "G4Step.hh"
103 
104 #include "G4Material.hh"
105 #include "G4Element.hh"
106 #include "G4Isotope.hh"
107 #include "G4MaterialCutsCouple.hh"
108 #include "G4ElementVector.hh"
109 #include "G4IsotopeVector.hh"
110 
111 #include "G4EmProcessSubType.hh"
112 
113 #include "G4ParticleDefinition.hh"
114 #include "G4DynamicParticle.hh"
115 #include "G4ProcessManager.hh"
116 #include "G4StableIsotopes.hh"
117 #include "G4LindhardPartition.hh"
118 
119 #include "Randomize.hh"
120 
121 #include <iostream>
122 #include <iomanip>
123 
124 #include "c2_factory.hh"
125 static c2_factory<G4double> c2; // this makes a lot of notation shorter
127 
129 {
130  screeningData.clear();
131  MFPTables.clear();
132 }
133 
134 const G4double G4ScreenedCoulombCrossSection::massmap[nMassMapElements+1]={
135  0, 1.007940, 4.002602, 6.941000, 9.012182, 10.811000, 12.010700,
136  14.006700, 15.999400, 18.998403, 20.179700, 22.989770, 24.305000, 26.981538,
137  28.085500,
138  30.973761, 32.065000, 35.453000, 39.948000, 39.098300, 40.078000, 44.955910,
139  47.867000,
140  50.941500, 51.996100, 54.938049, 55.845000, 58.933200, 58.693400, 63.546000,
141  65.409000,
142  69.723000, 72.640000, 74.921600, 78.960000, 79.904000, 83.798000, 85.467800,
143  87.620000,
144  88.905850, 91.224000, 92.906380, 95.940000, 98.000000, 101.070000, 102.905500,
145  106.420000,
146  107.868200, 112.411000, 114.818000, 118.710000, 121.760000, 127.600000,
147  126.904470, 131.293000,
148  132.905450, 137.327000, 138.905500, 140.116000, 140.907650, 144.240000,
149  145.000000, 150.360000,
150  151.964000, 157.250000, 158.925340, 162.500000, 164.930320, 167.259000,
151  168.934210, 173.040000,
152  174.967000, 178.490000, 180.947900, 183.840000, 186.207000, 190.230000,
153  192.217000, 195.078000,
154  196.966550, 200.590000, 204.383300, 207.200000, 208.980380, 209.000000,
155  210.000000, 222.000000,
156  223.000000, 226.000000, 227.000000, 232.038100, 231.035880, 238.028910,
157  237.000000, 244.000000,
158  243.000000, 247.000000, 247.000000, 251.000000, 252.000000, 257.000000,
159  258.000000, 259.000000,
160  262.000000, 261.000000, 262.000000, 266.000000, 264.000000, 277.000000,
161  268.000000, 281.000000,
162  272.000000, 285.000000, 282.500000, 289.000000, 287.500000, 292.000000};
163 
166  const G4MaterialCutsCouple* couple)
167 {
168  // Select randomly an element within the material, according to number
169  // density only
170  const G4Material* material = couple->GetMaterial();
171  G4int nMatElements = material->GetNumberOfElements();
172  const G4ElementVector* elementVector = material->GetElementVector();
173  const G4Element *element=0;
175 
176  // Special case: the material consists of one element
177  if (nMatElements == 1)
178  {
179  element= (*elementVector)[0];
180  }
181  else
182  {
183  // Composite material
184  G4double random = G4UniformRand() * material->GetTotNbOfAtomsPerVolume();
185  G4double nsum=0.0;
186  const G4double *atomDensities=material->GetVecNbOfAtomsPerVolume();
187 
188  for (G4int k=0 ; k < nMatElements ; k++ )
189  {
190  nsum+=atomDensities[k];
191  element= (*elementVector)[k];
192  if (nsum >= random) break;
193  }
194  }
195 
196  G4int N=0;
197  G4int Z=(G4int)std::floor(element->GetZ()+0.5);
198 
199  G4int nIsotopes=element->GetNumberOfIsotopes();
200  if(!nIsotopes) {
201  if(Z<=92) {
202  // we have no detailed material isotopic info available,
203  // so use G4StableIsotopes table up to Z=92
204  static G4StableIsotopes theIso;
205  // get a stable isotope table for default results
206  nIsotopes=theIso.GetNumberOfIsotopes(Z);
207  G4double random = 100.0*G4UniformRand();
208  // values are expressed as percent, sum is 100
209  G4int tablestart=theIso.GetFirstIsotope(Z);
210  G4double asum=0.0;
211  for(G4int i=0; i<nIsotopes; i++) {
212  asum+=theIso.GetAbundance(i+tablestart);
213  N=theIso.GetIsotopeNucleonCount(i+tablestart);
214  if(asum >= random) break;
215  }
216  } else {
217  // too heavy for stable isotope table, just use mean mass
218  N=(G4int)std::floor(element->GetN()+0.5);
219  }
220  } else {
221  G4int i;
222  const G4IsotopeVector *isoV=element->GetIsotopeVector();
223  G4double random = G4UniformRand();
224  G4double *abundance=element->GetRelativeAbundanceVector();
225  G4double asum=0.0;
226  for(i=0; i<nIsotopes; i++) {
227  asum+=abundance[i];
228  N=(*isoV)[i]->GetN();
229  if(asum >= random) break;
230  }
231  }
232 
233  // get the official definition of this nucleus, to get the correct
234  // value of A note that GetIon is very slow, so we will cache ones
235  // we have already found ourselves.
236  ParticleCache::iterator p=targetMap.find(Z*1000+N);
237  if (p != targetMap.end()) {
238  target=(*p).second;
239  } else{
240  target=G4IonTable::GetIonTable()->GetIon(Z, N, 0.0);
241  targetMap[Z*1000+N]=target;
242  }
243  return target;
244 }
245 
247 {
248  const G4int nmfpvals=200;
249 
250  std::vector<G4double> evals(nmfpvals), mfpvals(nmfpvals);
251 
252  // sum up inverse MFPs per element for each material
253  const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
254  if (materialTable == 0) { return; }
255  //G4Exception("G4ScreenedCoulombCrossSection::BuildMFPTables
256  //- no MaterialTable found)");
257 
259 
260  for (G4int matidx=0; matidx < nMaterials; matidx++) {
261 
262  const G4Material* material= (*materialTable)[matidx];
263  const G4ElementVector &elementVector =
264  *(material->GetElementVector());
265  const G4int nMatElements = material->GetNumberOfElements();
266 
267  const G4Element *element=0;
268  const G4double *atomDensities=material->GetVecNbOfAtomsPerVolume();
269 
270  G4double emin=0, emax=0;
271  // find innermost range of cross section functions
272  for (G4int kel=0 ; kel < nMatElements ; kel++ )
273  {
274  element=elementVector[kel];
275  G4int Z=(G4int)std::floor(element->GetZ()+0.5);
276  const G4_c2_function &ifunc=sigmaMap[Z];
277  if(!kel || ifunc.xmin() > emin) emin=ifunc.xmin();
278  if(!kel || ifunc.xmax() < emax) emax=ifunc.xmax();
279  }
280 
281  G4double logint=std::log(emax/emin) / (nmfpvals-1) ;
282  // logarithmic increment for tables
283 
284  // compute energy scale for interpolator. Force exact values at
285  // both ends to avoid range errors
286  for (G4int i=1; i<nmfpvals-1; i++) evals[i]=emin*std::exp(logint*i);
287  evals.front()=emin;
288  evals.back()=emax;
289 
290  // zero out the inverse mfp sums to start
291  for (G4int eidx=0; eidx < nmfpvals; eidx++) mfpvals[eidx] = 0.0;
292 
293  // sum inverse mfp for each element in this material and for each
294  // energy
295  for (G4int kel=0 ; kel < nMatElements ; kel++ )
296  {
297  element=elementVector[kel];
298  G4int Z=(G4int)std::floor(element->GetZ()+0.5);
299  const G4_c2_function &sigma=sigmaMap[Z];
300  G4double ndens = atomDensities[kel];
301  // compute atom fraction for this element in this material
302 
303  for (G4int eidx=0; eidx < nmfpvals; eidx++) {
304  mfpvals[eidx] += ndens*sigma(evals[eidx]);
305  }
306  }
307 
308  // convert inverse mfp to regular mfp
309  for (G4int eidx=0; eidx < nmfpvals; eidx++) {
310  mfpvals[eidx] = 1.0/mfpvals[eidx];
311  }
312  // and make a new interpolating function out of the sum
313  MFPTables[matidx] = c2.log_log_interpolating_function().load(evals,
314  mfpvals,true,0,true,0);
315  }
316 }
317 
319 G4ScreenedNuclearRecoil(const G4String& processName,
320  const G4String &ScreeningKey,
321  G4bool GenerateRecoils,
322  G4double RecoilCutoff, G4double PhysicsCutoff) :
323  G4VDiscreteProcess(processName, fElectromagnetic),
324  screeningKey(ScreeningKey),
325  generateRecoils(GenerateRecoils), avoidReactions(1),
326  recoilCutoff(RecoilCutoff), physicsCutoff(PhysicsCutoff),
327  hardeningFraction(0.0), hardeningFactor(1.0),
328  externalCrossSectionConstructor(0),
329  NIELPartitionFunction(new G4LindhardRobinsonPartition)
330 {
331  // for now, point to class instance of this. Doing it by creating a new
332  // one fails
333  // to correctly update NIEL
334  // not even this is needed... done in G4VProcess().
335  // pParticleChange=&aParticleChange;
336  processMaxEnergy=50000.0*MeV;
337  highEnergyLimit=100.0*MeV;
339  registerDepositedEnergy=1; // by default, don't hide NIEL
340  MFPScale=1.0;
341  // SetVerboseLevel(2);
345 }
346 
348 {
349 
350  std::map<G4int, G4ScreenedCoulombCrossSection*>::iterator xt=
351  crossSectionHandlers.begin();
352  for(;xt != crossSectionHandlers.end(); xt++) {
353  delete (*xt).second;
354  }
355  crossSectionHandlers.clear();
356 }
357 
359 {
360  // I don't think I like deleting the processes here... they are better
361  // abandoned
362  // if the creator doesn't get rid of them
363  // std::vector<G4ScreenedCollisionStage *>::iterator stage=
364  //collisionStages.begin();
365  //for(; stage != collisionStages.end(); stage++) delete (*stage);
366 
367  collisionStages.clear();
368 }
369 
371  const G4VNIELPartition *part)
372 {
375 }
376 
379 {
380  if(!NIELPartitionFunction) {
382  } else {
384  material, energy);
385  IonizingLoss+=energy*(1-part);
386  NIEL += energy*part;
387  }
388 }
389 
391 {
392  ResetTables();
393 }
394 
395 // returns true if it appears the nuclei collided, and we are interested
396 // in checking
398  G4double A, G4double a1, G4double apsis) {
399  return avoidReactions && (apsis < (1.1*(std::pow(A,1.0/3.0)+
400  std::pow(a1,1.0/3.0)) + 1.4)*fermi);
401  // nuclei are within 1.4 fm (reduced pion Compton wavelength) of each
402  // other at apsis,
403  // this is hadronic, skip it
404 }
405 
413  return xc;
414 }
415 
417  G4double,
418  G4ForceCondition* cond)
419 {
420  const G4DynamicParticle* incoming = track.GetDynamicParticle();
421  G4double energy = incoming->GetKineticEnergy();
422  G4double a1=incoming->GetDefinition()->GetPDGMass()/amu_c2;
423 
424  G4double meanFreePath;
425  *cond=NotForced;
426 
427  if (energy < lowEnergyLimit || energy < recoilCutoff*a1) {
428  *cond=Forced;
429  return 1.0*nm;
430 /* catch and stop slow particles to collect their NIEL! */
431  } else if (energy > processMaxEnergy*a1) {
432  return DBL_MAX; // infinite mean free path
433  } else if (energy > highEnergyLimit*a1) energy=highEnergyLimit*a1;
434 /* constant MFP at high energy */
435 
436  G4double fz1=incoming->GetDefinition()->GetPDGCharge();
437  G4int z1=(G4int)(fz1/eplus + 0.5);
438 
439  std::map<G4int, G4ScreenedCoulombCrossSection*>::iterator xh=
440  crossSectionHandlers.find(z1);
442 
443  if (xh==crossSectionHandlers.end()) {
445  xs->LoadData(screeningKey, z1, a1, physicsCutoff);
446  xs->BuildMFPTables();
447  } else xs=(*xh).second;
448 
449  const G4MaterialCutsCouple* materialCouple =
450  track.GetMaterialCutsCouple();
451  size_t materialIndex = materialCouple->GetMaterial()->GetIndex();
452 
453  const G4_c2_function &mfp=*(*xs)[materialIndex];
454 
455  // make absolutely certain we don't get an out-of-range energy
456  meanFreePath = mfp(std::min(std::max(energy, mfp.xmin()), mfp.xmax()));
457 
458  // G4cout << "MFP: " << meanFreePath << " index " << materialIndex
459  //<< " energy " << energy << " MFPScale " << MFPScale << G4endl;
460 
461  return meanFreePath*MFPScale;
462 }
463 
465  const G4Track& aTrack, const G4Step& aStep)
466 {
467  validCollision=1;
468  pParticleChange->Initialize(aTrack);
469  NIEL=0.0; // default is no NIEL deposited
470  IonizingLoss=0.0;
471 
472  // do universal setup
473 
474  const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
475  G4ParticleDefinition *baseParticle=aTrack.GetDefinition();
476 
477  G4double fz1=baseParticle->GetPDGCharge()/eplus;
478  G4int z1=(G4int)(fz1+0.5);
479  G4double a1=baseParticle->GetPDGMass()/amu_c2;
480  G4double incidentEnergy = incidentParticle->GetKineticEnergy();
481 
482  // Select randomly one element and (possibly) isotope in the
483  // current material.
484  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
485 
486  const G4Material* mat = couple->GetMaterial();
487 
488  G4double P=0.0; // the impact parameter of this collision
489 
490  if(incidentEnergy < GetRecoilCutoff()*a1) {
491  // check energy sanity on entry
492  DepositEnergy(z1, baseParticle->GetPDGMass()/amu_c2, mat,
493  incidentEnergy);
495  // stop the particle and bail out
496  validCollision=0;
497  } else {
498 
499  G4double numberDensity=mat->GetTotNbOfAtomsPerVolume();
500  G4double lattice=0.5/std::pow(numberDensity,1.0/3.0);
501  // typical lattice half-spacing
503  G4double sigopi=1.0/(pi*numberDensity*length);
504  // this is sigma0/pi
505 
506  // compute the impact parameter very early, so if is rejected
507  // as too far away, little effort is wasted
508  // this is the TRIM method for determining an impact parameter
509  // based on the flight path
510  // this gives a cumulative distribution of
511  // N(P)= 1-exp(-pi P^2 n l)
512  // which says the probability of NOT hitting a disk of area
513  // sigma= pi P^2 =exp(-sigma N l)
514  // which may be reasonable
515  if(sigopi < lattice*lattice) {
516  // normal long-flight approximation
517  P = std::sqrt(-std::log(G4UniformRand()) *sigopi);
518  } else {
519  // short-flight limit
520  P = std::sqrt(G4UniformRand())*lattice;
521  }
522 
523  G4double fraction=GetHardeningFraction();
524  if(fraction && G4UniformRand() < fraction) {
525  // pick out some events, and increase the central cross
526  // section by reducing the impact parameter
527  P /= std::sqrt(GetHardeningFactor());
528  }
529 
530 
531  // check if we are far enough away that the energy transfer
532  // must be below cutoff,
533  // and leave everything alone if so, saving a lot of time.
534  if(P*P > sigopi) {
535  if(GetVerboseLevel() > 1)
536  printf("ScreenedNuclear impact reject: length=%.3f P=%.4f limit=%.4f\n",
537  length/angstrom, P/angstrom,std::sqrt(sigopi)/angstrom);
538  // no collision, don't follow up with anything
539  validCollision=0;
540  }
541  }
542 
543  // find out what we hit, and record it in our kinematics block.
545  kinematics.a1=a1;
546 
547  if(validCollision) {
550  G4ParticleDefinition *recoilIon=
551  xsect->SelectRandomUnweightedTarget(couple);
552  kinematics.crossSection=xsect;
553  kinematics.recoilIon=recoilIon;
555  kinematics.a2=recoilIon->GetPDGMass()/amu_c2;
556  } else {
559  kinematics.a2=0;
560  }
561 
562  std::vector<G4ScreenedCollisionStage *>::iterator stage=
563  collisionStages.begin();
564 
565  for(; stage != collisionStages.end(); stage++)
566  (*stage)->DoCollisionStep(this,aTrack, aStep);
567 
571  //MHM G4cout << "depositing energy, total = "
572  //<< IonizingLoss+NIEL << " NIEL = " << NIEL << G4endl;
573  }
574 
575  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
576 }
577 
579 // instantiate all the needed functions statically, so no allocation is
580 // done at run time
581 // we will be solving x^2 - x phi(x*au)/eps - beta^2 == 0.0
582 // or, for easier scaling, x'^2 - x' au phi(x')/eps - beta^2 au^2
583 // note that only the last of these gets deleted, since it owns the rest
584  phifunc(c2.const_plugin_function()),
585  xovereps(c2.linear(0., 0., 0.)),
586  // will fill this in with the right slope at run time
587  diff(c2.quadratic(0., 0., 0., 1.)-xovereps*phifunc)
588 {
589 }
590 
592  G4ScreenedNuclearRecoil *master,
593  const G4ScreeningTables *screen, G4double eps, G4double beta)
594 {
595  G4double au=screen->au;
596  G4CoulombKinematicsInfo &kin=master->GetKinematics();
597  G4double A=kin.a2;
598  G4double a1=kin.a1;
599 
600  G4double xx0; // first estimate of closest approach
601  if(eps < 5.0) {
602  G4double y=std::log(eps);
603  G4double mlrho4=((((3.517e-4*y+1.401e-2)*y+2.393e-1)*y+2.734)*y+2.220);
604  G4double rho4=std::exp(-mlrho4); // W&M eq. 18
605  G4double bb2=0.5*beta*beta;
606  xx0=std::sqrt(bb2+std::sqrt(bb2*bb2+rho4)); // W&M eq. 17
607  } else {
608  G4double ee=1.0/(2.0*eps);
609  xx0=ee+std::sqrt(ee*ee+beta*beta); // W&M eq. 15 (Rutherford value)
610  if(master->CheckNuclearCollision(A, a1, xx0*au)) return 0;
611  // nuclei too close
612 
613  }
614 
615  // we will be solving x^2 - x phi(x*au)/eps - beta^2 == 0.0
616  // or, for easier scaling, x'^2 - x' au phi(x')/eps - beta^2 au^2
617  xovereps.reset(0., 0.0, au/eps); // slope of x*au/eps term
618  phifunc.set_function(&(screen->EMphiData.get()));
619  // install interpolating table
620  G4double xx1, phip, phip2;
621  G4int root_error;
622  xx1=diff->find_root(phifunc.xmin(), std::min(10*xx0*au,phifunc.xmax()),
623  std::min(xx0*au, phifunc.xmax()), beta*beta*au*au,
624  &root_error, &phip, &phip2)/au;
625 
626  if(root_error) {
627  G4cout << "Screened Coulomb Root Finder Error" << G4endl;
628  G4cout << "au " << au << " A " << A << " a1 " << a1
629  << " xx1 " << xx1 << " eps " << eps
630  << " beta " << beta << G4endl;
631  G4cout << " xmin " << phifunc.xmin() << " xmax "
632  << std::min(10*xx0*au,phifunc.xmax()) ;
633  G4cout << " f(xmin) " << phifunc(phifunc.xmin())
634  << " f(xmax) "
635  << phifunc(std::min(10*xx0*au,phifunc.xmax())) ;
636  G4cout << " xstart " << std::min(xx0*au, phifunc.xmax())
637  << " target " << beta*beta*au*au ;
638  G4cout << G4endl;
639  throw c2_exception("Failed root find");
640  }
641 
642  // phiprime is scaled by one factor of au because phi is evaluated
643  // at (xx0*au),
644  G4double phiprime=phip*au;
645 
646  //lambda0 is from W&M 19
647  G4double lambda0=1.0/std::sqrt(0.5+beta*beta/(2.0*xx1*xx1)
648  -phiprime/(2.0*eps));
649 
650  // compute the 6-term Lobatto integral alpha (per W&M 21, with
651  // different coefficients)
652  // this is probably completely un-needed but gives the highest
653  // quality results,
654  G4double alpha=(1.0+ lambda0)/30.0;
655  G4double xvals[]={0.98302349, 0.84652241, 0.53235309, 0.18347974};
656  G4double weights[]={0.03472124, 0.14769029, 0.23485003, 0.18602489};
657  for(G4int k=0; k<4; k++) {
658  G4double x, ff;
659  x=xx1/xvals[k];
660  ff=1.0/std::sqrt(1.0-phifunc(x*au)/(x*eps)-beta*beta/(x*x));
661  alpha+=weights[k]*ff;
662  }
663 
665  // throws an exception if used without setting again
666 
667  G4double thetac1=pi*beta*alpha/xx1;
668  // complement of CM scattering angle
669  G4double sintheta=std::sin(thetac1); //note sin(pi-theta)=sin(theta)
670  G4double costheta=-std::cos(thetac1); // note cos(pi-theta)=-cos(theta)
671  // G4double psi=std::atan2(sintheta, costheta+a1/A);
672  // lab scattering angle (M&T 3rd eq. 8.69)
673 
674  // numerics note: because we checked above for reasonable values
675  // of beta which give real recoils,
676  // we don't have to look too closely for theta -> 0 here
677  // (which would cause sin(theta)
678  // and 1-cos(theta) to both vanish and make the atan2 ill behaved).
679  G4double zeta=std::atan2(sintheta, 1-costheta);
680  // lab recoil angle (M&T 3rd eq. 8.73)
681  G4double coszeta=std::cos(zeta);
682  G4double sinzeta=std::sin(zeta);
683 
684  kin.sinTheta=sintheta;
685  kin.cosTheta=costheta;
686  kin.sinZeta=sinzeta;
687  kin.cosZeta=coszeta;
688  return 1; // all OK, collision is valid
689 }
690 
692  G4ScreenedNuclearRecoil *master,
693  const G4Track& aTrack, const G4Step&) {
694 
695  if(!master->GetValidCollision()) return;
696 
697  G4ParticleChange &aParticleChange=master->GetParticleChange();
698  G4CoulombKinematicsInfo &kin=master->GetKinematics();
699 
700  const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
701  G4ParticleDefinition *baseParticle=aTrack.GetDefinition();
702 
703  G4double incidentEnergy = incidentParticle->GetKineticEnergy();
704 
705  // this adjustment to a1 gives the right results for soft
706  // (constant gamma)
707  // relativistic collisions. Hard collisions are wrong anyway, since the
708  // Coulombic and hadronic terms interfere and cannot be added.
709  G4double gamma=(1.0+incidentEnergy/baseParticle->GetPDGMass());
710  G4double a1=kin.a1*gamma; // relativistic gamma correction
711 
712  G4ParticleDefinition *recoilIon=kin.recoilIon;
713  G4double A=recoilIon->GetPDGMass()/amu_c2;
714  G4int Z=(G4int)((recoilIon->GetPDGCharge()/eplus)+0.5);
715 
716  G4double Ec = incidentEnergy*(A/(A+a1));
717  // energy in CM frame (non-relativistic!)
718  const G4ScreeningTables *screen=kin.crossSection->GetScreening(Z);
719  G4double au=screen->au; // screening length
720 
721  G4double beta = kin.impactParameter/au;
722  // dimensionless impact parameter
723  G4double eps = Ec/(screen->z1*Z*elm_coupling/au);
724  // dimensionless energy
725 
726  G4bool ok=DoScreeningComputation(master, screen, eps, beta);
727  if(!ok) {
728  master->SetValidCollision(0); // flag bad collision
729  return; // just bail out without setting valid flag
730  }
731 
732  G4double eRecoil=4*incidentEnergy*a1*A*kin.cosZeta*kin.cosZeta
733  /((a1+A)*(a1+A));
734  kin.eRecoil=eRecoil;
735 
736  if(incidentEnergy-eRecoil < master->GetRecoilCutoff()*a1) {
737  aParticleChange.ProposeEnergy(0.0);
738  master->DepositEnergy(int(screen->z1), a1, kin.targetMaterial,
739  incidentEnergy-eRecoil);
740  }
741 
742  if(master->GetEnableRecoils() &&
743  eRecoil > master->GetRecoilCutoff() * kin.a2) {
744  kin.recoilIon=recoilIon;
745  } else {
746  kin.recoilIon=0; // this flags no recoil to be generated
747  master->DepositEnergy(Z, A, kin.targetMaterial, eRecoil) ;
748  }
749 }
750 
752  const G4Track& aTrack, const G4Step&) {
753 
754  if(!master->GetValidCollision()) return;
755 
756  G4CoulombKinematicsInfo &kin=master->GetKinematics();
757  G4ParticleChange &aParticleChange=master->GetParticleChange();
758 
759  const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
760  G4double incidentEnergy = incidentParticle->GetKineticEnergy();
761  G4double eRecoil=kin.eRecoil;
762 
763  G4double azimuth=G4UniformRand()*(2.0*pi);
764  G4double sa=std::sin(azimuth);
765  G4double ca=std::cos(azimuth);
766 
767  G4ThreeVector recoilMomentumDirection(kin.sinZeta*ca,
768  kin.sinZeta*sa, kin.cosZeta);
769  G4ParticleMomentum incidentDirection =
770  incidentParticle->GetMomentumDirection();
771  recoilMomentumDirection=
772  recoilMomentumDirection.rotateUz(incidentDirection);
773  G4ThreeVector recoilMomentum=
774  recoilMomentumDirection*std::sqrt(2.0*eRecoil*kin.a2*amu_c2);
775 
776  if(aParticleChange.GetEnergy() != 0.0) {
777  // DoKinematics hasn't stopped it!
778  G4ThreeVector beamMomentum=
779  incidentParticle->GetMomentum()-recoilMomentum;
780  aParticleChange.ProposeMomentumDirection(beamMomentum.unit()) ;
781  aParticleChange.ProposeEnergy(incidentEnergy-eRecoil);
782  }
783 
784  if(kin.recoilIon) {
785  G4DynamicParticle* recoil =
786  new G4DynamicParticle (kin.recoilIon,
787  recoilMomentumDirection,eRecoil) ;
788 
789  aParticleChange.SetNumberOfSecondaries(1);
790  aParticleChange.AddSecondary(recoil);
791  }
792 }
793 
795 IsApplicable(const G4ParticleDefinition& aParticleType)
796 {
797  return aParticleType == *(G4Proton::Proton()) ||
798  aParticleType.GetParticleType() == "nucleus" ||
799  aParticleType.GetParticleType() == "static_nucleus";
800 }
801 
802 void
805 {
806  G4String nam = aParticleType.GetParticleName();
807  if(nam == "GenericIon" || nam == "proton"
808  || nam == "deuteron" || nam == "triton"
809  || nam == "alpha" || nam == "He3") {
810  G4cout << G4endl << GetProcessName() << ": for " << nam
811  << " SubType= " << GetProcessSubType()
812  << " maxEnergy(MeV)= " << processMaxEnergy/MeV << G4endl;
813  }
814 }
815 
816 void
819 {
820 }
821 
822 // This used to be the file mhmScreenedNuclearRecoil_native.cc
823 // it has been included here to collect this file into a smaller
824 // number of packages
825 
826 #include "G4DataVector.hh"
827 #include "G4Material.hh"
828 #include "G4Element.hh"
829 #include "G4Isotope.hh"
830 #include "G4MaterialCutsCouple.hh"
831 #include "G4ElementVector.hh"
832 #include <vector>
833 
835  G4double rMax, G4double *auval)
836 {
837  static const size_t ncoef=4;
838  static G4double scales[ncoef]={-3.2, -0.9432, -0.4028, -0.2016};
839  static G4double coefs[ncoef]={0.1818,0.5099,0.2802,0.0281};
840 
841  G4double au=
842  0.8854*angstrom*0.529/(std::pow(z1, 0.23)+std::pow(z2,0.23));
843  std::vector<G4double> r(npoints), phi(npoints);
844 
845  for(size_t i=0; i<npoints; i++) {
846  G4double rr=(float)i/(float)(npoints-1);
847  r[i]=rr*rr*rMax;
848  // use quadratic r scale to make sampling fine near the center
849  G4double sum=0.0;
850  for(size_t j=0; j<ncoef; j++)
851  sum+=coefs[j]*std::exp(scales[j]*r[i]/au);
852  phi[i]=sum;
853  }
854 
855  // compute the derivative at the origin for the spline
856  G4double phiprime0=0.0;
857  for(size_t j=0; j<ncoef; j++)
858  phiprime0+=scales[j]*coefs[j]*std::exp(scales[j]*r[0]/au);
859  phiprime0*=(1.0/au); // put back in natural units;
860 
861  *auval=au;
862  return c2.lin_log_interpolating_function().load(r, phi, false,
863  phiprime0,true,0);
864 }
865 
867  G4double rMax, G4double *auval)
868 {
869  static const size_t ncoef=3;
870  static G4double scales[ncoef]={-6.0, -1.2, -0.3};
871  static G4double coefs[ncoef]={0.10, 0.55, 0.35};
872 
873  G4double au=0.8853*0.529*angstrom/std::sqrt(std::pow(z1, 0.6667)
874  +std::pow(z2,0.6667));
875  std::vector<G4double> r(npoints), phi(npoints);
876 
877  for(size_t i=0; i<npoints; i++) {
878  G4double rr=(float)i/(float)(npoints-1);
879  r[i]=rr*rr*rMax;
880  // use quadratic r scale to make sampling fine near the center
881  G4double sum=0.0;
882  for(size_t j=0; j<ncoef; j++)
883  sum+=coefs[j]*std::exp(scales[j]*r[i]/au);
884  phi[i]=sum;
885  }
886 
887  // compute the derivative at the origin for the spline
888  G4double phiprime0=0.0;
889  for(size_t j=0; j<ncoef; j++)
890  phiprime0+=scales[j]*coefs[j]*std::exp(scales[j]*r[0]/au);
891  phiprime0*=(1.0/au); // put back in natural units;
892 
893  *auval=au;
894  return c2.lin_log_interpolating_function().load(r, phi, false,
895  phiprime0,true,0);
896 }
897 
899  G4double rMax, G4double *auval)
900 {
901 //from Loftager, Besenbacher, Jensen & Sorensen
902 //PhysRev A20, 1443++, 1979
903  G4double au=0.8853*0.529*angstrom/std::sqrt(std::pow(z1, 0.6667)
904  +std::pow(z2,0.6667));
905  std::vector<G4double> r(npoints), phi(npoints);
906 
907  for(size_t i=0; i<npoints; i++) {
908  G4double rr=(float)i/(float)(npoints-1);
909  r[i]=rr*rr*rMax;
910  // use quadratic r scale to make sampling fine near the center
911 
912  G4double y=std::sqrt(9.67*r[i]/au);
913  G4double ysq=y*y;
914  G4double phipoly=1+y+0.3344*ysq+0.0485*y*ysq+0.002647*ysq*ysq;
915  phi[i]=phipoly*std::exp(-y);
916  // G4cout << r[i] << " " << phi[i] << G4endl;
917  }
918 
919  // compute the derivative at the origin for the spline
920  G4double logphiprime0=(9.67/2.0)*(2*0.3344-1.0);
921  // #avoid 0/0 on first element
922  logphiprime0 *= (1.0/au); // #put back in natural units
923 
924  *auval=au;
925  return c2.lin_log_interpolating_function().load(r, phi, false,
926  logphiprime0*phi[0],
927  true,0);
928 }
929 
931  G4double rMax, G4double *auval)
932 {
933 // hybrid of LJ and ZBL, uses LJ if x < 0.25*auniv, ZBL if x > 1.5*auniv, and
935 // is very near the point where the functions naturally cross.
936  G4double auzbl, aulj;
937 
938  c2p zbl=ZBLScreening(z1, z2, npoints, rMax, &auzbl);
939  c2p lj=LJScreening(z1, z2, npoints, rMax, &aulj);
940 
941  G4double au=(auzbl+aulj)*0.5;
942  lj->set_domain(lj->xmin(), 0.25*au);
943  zbl->set_domain(1.5*au,zbl->xmax());
944 
945  c2p conn=
946  c2.connector_function(lj->xmax(), lj, zbl->xmin(), zbl, true,0);
948  c2p keepit(pw);
949  pw.append_function(lj);
950  pw.append_function(conn);
951  pw.append_function(zbl);
952 
953  *auval=au;
954  keepit.release_for_return();
955  return pw;
956 }
957 
959 }
960 
966 }
967 
968 std::vector<G4String>
970  std::vector<G4String> keys;
971  // find the available screening keys
972  std::map<std::string, ScreeningFunc>::const_iterator sfunciter=phiMap.begin();
973  for(; sfunciter != phiMap.end(); sfunciter++)
974  keys.push_back((*sfunciter).first);
975  return keys;
976 }
977 
978 static inline G4double cm_energy(G4double a1, G4double a2, G4double t0) {
979  // "relativistically correct energy in CM frame"
980  G4double m1=a1*amu_c2, mass2=a2*amu_c2;
981  G4double mc2=(m1+mass2);
982  G4double f=2.0*mass2*t0/(mc2*mc2);
983  // old way: return (f < 1e-6) ? 0.5*mc2*f : mc2*(std::sqrt(1.0+f)-1.0);
984  // formally equivalent to previous, but numerically stable for all
985  // f without conditional
986  // uses identity (sqrt(1+x) - 1)(sqrt(1+x) + 1) = x
987  return mc2*f/(std::sqrt(1.0+f)+1.0);
988 }
989 
990 static inline G4double thetac(G4double m1, G4double mass2, G4double eratio) {
991  G4double s2th2=eratio*( (m1+mass2)*(m1+mass2)/(4.0*m1*mass2) );
992  G4double sth2=std::sqrt(s2th2);
993  return 2.0*std::asin(sth2);
994 }
995 
997  G4int z1, G4double a1,
998  G4double recoilCutoff)
999 {
1000  static const size_t sigLen=200;
1001  // since sigma doesn't matter much, a very coarse table will do
1002  G4DataVector energies(sigLen);
1003  G4DataVector data(sigLen);
1004 
1005  a1=standardmass(z1);
1006  // use standardized values for mass for building tables
1007 
1008  const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
1010 
1011  for (G4int im=0; im<nMaterials; im++)
1012  {
1013  const G4Material* material= (*materialTable)[im];
1014  const G4ElementVector* elementVector = material->GetElementVector();
1015  const G4int nMatElements = material->GetNumberOfElements();
1016 
1017  for (G4int iEl=0; iEl<nMatElements; iEl++)
1018  {
1019  G4Element* element = (*elementVector)[iEl];
1020  G4int Z = (G4int) element->GetZ();
1021  G4double a2=element->GetA()*(mole/gram);
1022 
1023  if(sigmaMap.find(Z)!=sigmaMap.end()) continue;
1024  // we've already got this element
1025 
1026  // find the screening function generator we need
1027  std::map<std::string, ScreeningFunc>::iterator sfunciter=
1028  phiMap.find(screeningKey);
1029  if(sfunciter==phiMap.end()) {
1031  ed << "No such screening key <"
1032  << screeningKey << ">";
1033  G4Exception("G4NativeScreenedCoulombCrossSection::LoadData",
1034  "em0003",FatalException,ed);
1035  }
1036  ScreeningFunc sfunc=(*sfunciter).second;
1037 
1038  G4double au;
1039  G4_c2_ptr screen=sfunc(z1, Z, 200, 50.0*angstrom, &au);
1040  // generate the screening data
1041  G4ScreeningTables st;
1042 
1043  st.EMphiData=screen; //save our phi table
1044  st.z1=z1; st.m1=a1; st.z2=Z; st.m2=a2; st.emin=recoilCutoff;
1045  st.au=au;
1046 
1047  // now comes the hard part... build the total cross section
1048  // tables from the phi table
1049  // based on (pi-thetac) = pi*beta*alpha/x0, but noting that
1050  // alpha is very nearly unity, always
1051  // so just solve it wth alpha=1, which makes the solution
1052  // much easier
1053  // this function returns an approximation to
1054  // (beta/x0)^2=phi(x0)/(eps*x0)-1 ~ ((pi-thetac)/pi)^2
1055  // Since we don't need exact sigma values, this is good enough
1056  // (within a factor of 2 almost always)
1057  // this rearranges to phi(x0)/(x0*eps) =
1058  // 2*theta/pi - theta^2/pi^2
1059 
1060  c2_linear_p<G4double> &c2eps=c2.linear(0.0, 0.0, 1.0);
1061  // will store an appropriate eps inside this in loop
1062  G4_c2_ptr phiau=screen(c2.linear(0.0, 0.0, au));
1063  G4_c2_ptr x0func(phiau/c2eps);
1064  // this will be phi(x)/(x*eps) when c2eps is correctly set
1065  x0func->set_domain(1e-6*angstrom/au, 0.9999*screen->xmax()/au);
1066  // needed for inverse function
1067  // use the c2_inverse_function interface for the root finder
1068  // it is more efficient for an ordered
1069  // computation of values.
1070  G4_c2_ptr x0_solution(c2.inverse_function(x0func));
1071 
1072  G4double m1c2=a1*amu_c2;
1073  G4double escale=z1*Z*elm_coupling/au;
1074  // energy at screening distance
1075  G4double emax=m1c2;
1076  // model is doubtful in very relativistic range
1077  G4double eratkin=0.9999*(4*a1*a2)/((a1+a2)*(a1+a2));
1078  // #maximum kinematic ratio possible at 180 degrees
1079  G4double cmfact0=st.emin/cm_energy(a1, a2, st.emin);
1080  G4double l1=std::log(emax);
1081  G4double l0=std::log(st.emin*cmfact0/eratkin);
1082 
1083  if(verbosity >=1)
1084  G4cout << "Native Screening: " << screeningKey << " "
1085  << z1 << " " << a1 << " " <<
1086  Z << " " << a2 << " " << recoilCutoff << G4endl;
1087 
1088  for(size_t idx=0; idx<sigLen; idx++) {
1089  G4double ee=std::exp(idx*((l1-l0)/sigLen)+l0);
1090  G4double gamma=1.0+ee/m1c2;
1091  G4double eratio=(cmfact0*st.emin)/ee;
1092  // factor by which ee needs to be reduced to get emin
1093  G4double theta=thetac(gamma*a1, a2, eratio);
1094 
1095  G4double eps=cm_energy(a1, a2, ee)/escale;
1096  // #make sure lab energy is converted to CM for these
1097  // calculations
1098  c2eps.reset(0.0, 0.0, eps);
1099  // set correct slope in this function
1100 
1101  G4double q=theta/pi;
1102  // G4cout << ee << " " << m1c2 << " " << gamma << " "
1103  // << eps << " " << theta << " " << q << G4endl;
1104  // old way using root finder
1105  // G4double x0= x0func->find_root(1e-6*angstrom/au,
1106  // 0.9999*screen.xmax()/au, 1.0, 2*q-q*q);
1107  // new way using c2_inverse_function which caches
1108  // useful information so should be a bit faster
1109  // since we are scanning this in strict order.
1110  G4double x0=0;
1111  try {
1112  x0=x0_solution(2*q-q*q);
1113  } catch(c2_exception&) {
1114  G4Exception("G4ScreenedNuclearRecoil::LoadData",
1115  "em0003",FatalException,
1116  "failure in inverse solution to generate MFP tables");
1117  }
1118  G4double betasquared=x0*x0 - x0*phiau(x0)/eps;
1119  G4double sigma=pi*betasquared*au*au;
1120  energies[idx]=ee;
1121  data[idx]=sigma;
1122  }
1123  screeningData[Z]=st;
1124  sigmaMap[Z] =
1125  c2.log_log_interpolating_function().load(energies, data,
1126  true,0,true,0);
1127  }
1128  }
1129 }