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