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