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