93 "G4ScreenedNuclearRecoil.cc,v 1.57 2008/05/07 11:51:26 marcus Exp GEANT4 tag ";
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,
141 30.973761, 32.065000, 35.453000, 39.948000, 39.098300, 40.078000, 44.955910,
143 50.941500, 51.996100, 54.938049, 55.845000, 58.933200, 58.693400, 63.546000,
145 69.723000, 72.640000, 74.921600, 78.960000, 79.904000, 83.798000, 85.467800,
147 88.905850, 91.224000, 92.906380, 95.940000, 98.000000, 101.070000, 102.905500,
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};
180 if (nMatElements == 1)
182 element= (*elementVector)[0];
191 for (
G4int k=0 ; k < nMatElements ; k++ )
193 nsum+=atomDensities[k];
194 element= (*elementVector)[k];
195 if (nsum >= random)
break;
214 for(
G4int i=0; i<nIsotopes; i++) {
217 if(asum >= random)
break;
221 N=(
G4int)std::floor(element->
GetN()+0.5);
229 for(i=0; i<nIsotopes; i++) {
231 N=(*isoV)[i]->GetN();
232 if(asum >= random)
break;
239 ParticleCache::iterator p=
targetMap.find(Z*1000+N);
251 const G4int nmfpvals=200;
253 std::vector<G4double> evals(nmfpvals), mfpvals(nmfpvals);
257 if (materialTable == 0) {
return; }
263 for (
G4int matidx=0; matidx < nMaterials; matidx++) {
265 const G4Material* material= (*materialTable)[matidx];
275 for (
G4int kel=0 ; kel < nMatElements ; kel++ )
277 element=elementVector[kel];
280 if(!kel || ifunc.
xmin() > emin) emin=ifunc.
xmin();
281 if(!kel || ifunc.
xmax() < emax) emax=ifunc.
xmax();
284 G4double logint=std::log(emax/emin) / (nmfpvals-1) ;
289 for (
G4int i=1; i<nmfpvals-1; i++) evals[i]=emin*std::exp(logint*i);
294 for (
G4int eidx=0; eidx < nmfpvals; eidx++) mfpvals[eidx] = 0.0;
298 for (
G4int kel=0 ; kel < nMatElements ; kel++ )
300 element=elementVector[kel];
303 G4double ndens = atomDensities[kel];
306 for (
G4int eidx=0; eidx < nmfpvals; eidx++) {
307 mfpvals[eidx] += ndens*sigma(evals[eidx]);
312 for (
G4int eidx=0; eidx < nmfpvals; eidx++) {
313 mfpvals[eidx] = 1.0/mfpvals[eidx];
317 mfpvals,
true,0,
true,0);
327 screeningKey(ScreeningKey),
328 generateRecoils(GenerateRecoils), avoidReactions(1),
329 recoilCutoff(RecoilCutoff), physicsCutoff(PhysicsCutoff),
330 hardeningFraction(0.0), hardeningFactor(1.0),
331 externalCrossSectionConstructor(0),
353 std::map<G4int, G4ScreenedCoulombCrossSection*>::iterator xt=
403 std::pow(a1,1.0/3.0)) + 1.4)*
fermi);
442 std::map<G4int, G4ScreenedCoulombCrossSection*>::iterator xh=
450 }
else xs=(*xh).second;
503 G4double lattice=0.5/std::pow(numberDensity,1.0/3.0);
518 if(sigopi < lattice*lattice) {
539 printf(
"ScreenedNuclear impact reject: length=%.3f P=%.4f limit=%.4f\n",
565 std::vector<G4ScreenedCollisionStage *>::iterator stage=
569 (*stage)->DoCollisionStep(
this,aTrack, aStep);
587 phifunc(c2.const_plugin_function()),
588 xovereps(c2.linear(0., 0., 0.)),
590 diff(c2.quadratic(0., 0., 0., 1.)-xovereps*phifunc)
606 G4double mlrho4=((((3.517e-4*y+1.401e-2)*y+2.393e-1)*y+2.734)*y+2.220);
609 xx0=std::sqrt(bb2+std::sqrt(bb2*bb2+rho4));
612 xx0=ee+std::sqrt(ee*ee+beta*beta);
627 &root_error, &phip, &phip2)/au;
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;
640 <<
" target " << beta*beta*au*au ;
650 G4double lambda0=1.0/std::sqrt(0.5+beta*beta/(2.0*xx1*xx1)
651 -phiprime/(2.0*eps));
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++) {
663 ff=1.0/std::sqrt(1.0-
phifunc(x*au)/(x*eps)-beta*beta/(x*x));
664 alpha+=weights[k]*ff;
672 G4double sintheta=std::sin(thetac1);
673 G4double costheta=-std::cos(thetac1);
682 G4double zeta=std::atan2(sintheta, 1-costheta);
739 if(incidentEnergy-eRecoil < master->GetRecoilCutoff()*
a1) {
742 incidentEnergy-eRecoil);
774 recoilMomentumDirection=
775 recoilMomentumDirection.rotateUz(incidentDirection);
777 recoilMomentumDirection*std::sqrt(2.0*eRecoil*kin.
a2*amu_c2);
790 recoilMomentumDirection,eRecoil) ;
810 if(nam ==
"GenericIon" || nam ==
"proton"
811 || nam ==
"deuteron" || nam ==
"triton"
812 || nam ==
"alpha" || nam ==
"He3") {
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};
845 0.8854*
angstrom*0.529/(std::pow(z1, 0.23)+std::pow(z2,0.23));
846 std::vector<G4double> r(npoints), phi(npoints);
848 for(
size_t i=0; i<npoints; i++) {
849 G4double rr=(float)i/(
float)(npoints-1);
853 for(
size_t j=0; j<ncoef; j++)
854 sum+=coefs[j]*std::exp(scales[j]*r[i]/au);
860 for(
size_t j=0; j<ncoef; j++)
861 phiprime0+=scales[j]*coefs[j]*std::exp(scales[j]*r[0]/au);
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};
877 +std::pow(z2,0.6667));
878 std::vector<G4double> r(npoints), phi(npoints);
880 for(
size_t i=0; i<npoints; i++) {
881 G4double rr=(float)i/(
float)(npoints-1);
885 for(
size_t j=0; j<ncoef; j++)
886 sum+=coefs[j]*std::exp(scales[j]*r[i]/au);
892 for(
size_t j=0; j<ncoef; j++)
893 phiprime0+=scales[j]*coefs[j]*std::exp(scales[j]*r[0]/au);
907 +std::pow(z2,0.6667));
908 std::vector<G4double> r(npoints), phi(npoints);
910 for(
size_t i=0; i<npoints; i++) {
911 G4double rr=(float)i/(
float)(npoints-1);
917 G4double phipoly=1+y+0.3344*ysq+0.0485*y*ysq+0.002647*ysq*ysq;
918 phi[i]=phipoly*std::exp(-y);
923 G4double logphiprime0=(9.67/2.0)*(2*0.3344-1.0);
925 logphiprime0 *= (1.0/au);
971 std::vector<G4String>
973 std::vector<G4String> keys;
975 std::map<std::string, ScreeningFunc>::const_iterator sfunciter=
phiMap.begin();
976 for(; sfunciter !=
phiMap.end(); sfunciter++)
977 keys.push_back((*sfunciter).first);
983 G4double m1=a1*amu_c2, mass2=a2*amu_c2;
990 return mc2*f/(std::sqrt(1.0+f)+1.0);
994 G4double s2th2=eratio*( (m1+mass2)*(m1+mass2)/(4.0*m1*mass2) );
996 return 2.0*std::asin(sth2);
1003 static const size_t sigLen=200;
1012 if (materialTable == 0) {
return; }
1017 for (
G4int im=0; im<nMaterials; im++)
1019 const G4Material* material= (*materialTable)[im];
1023 for (
G4int iEl=0; iEl<nMatElements; iEl++)
1025 G4Element* element = (*elementVector)[iEl];
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",
1071 x0func->set_domain(1e-6*
angstrom/au, 0.9999*screen->
xmax()/au);
1079 G4double escale=z1*Z*elm_coupling/au;
1090 G4cout <<
"Native Screening: " << screeningKey <<
" "
1091 << z1 <<
" " << a1 <<
" " <<
1092 Z <<
" " << a2 <<
" " << recoilCutoff <<
G4endl;
1094 for(
size_t idx=0; idx<sigLen; idx++) {
1095 G4double ee=std::exp(idx*((l1-l0)/sigLen)+l0);
1104 c2eps.
reset(0.0, 0.0, eps);
1118 x0=x0_solution(2*q-q*q);
size_t GetNumberOfIsotopes() const
G4ParticleDefinition * GetDefinition() const
virtual void Initialize(const G4Track &)
static c2_factory< G4double > c2
G4_c2_const_ptr EMphiData
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
G4CoulombKinematicsInfo kinematics
G4double GetCurrentInteractionLength() const
the the interaciton length used in the last scattering.
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
void AddStage(G4ScreenedCollisionStage *stage)
virtual void LoadData(G4String screeningKey, G4int z1, G4double m1, G4double recoilCutoff)=0
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4bool GetValidCollision() const
const G4DynamicParticle * GetDynamicParticle() const
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()
G4int GetFirstIsotope(G4int Z)
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)
static G4MaterialTable * GetMaterialTable()
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 BuildMFPTables(void)
virtual ~G4NativeScreenedCoulombCrossSection()
void set_function(const c2_function< float_type > *f)
fill the container with a new function, or clear it with a null pointer
G4_c2_function & MoliereScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
const G4ElementVector * GetElementVector() const
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
const G4double * GetVecNbOfAtomsPerVolume() const
static G4double thetac(G4double m1, G4double mass2, G4double eratio)
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
std::map< G4int, G4ScreenedCoulombCrossSection * > crossSectionHandlers
const G4ThreeVector & GetMomentumDirection() const
void set_domain(float_type amin, float_type amax)
set the domain for this function.
ScreeningMap screeningData
virtual ~G4ScreenedCoulombCrossSection()
void SetValidCollision(G4bool flag)
const G4ScreeningTables * GetScreening(G4int Z)
G4ParticleDefinition * SelectRandomUnweightedTarget(const G4MaterialCutsCouple *couple)
G4NativeScreenedCoulombCrossSection()
G4double * GetRelativeAbundanceVector() const
static G4Proton * Proton()
void SetProcessSubType(G4int)
void release_for_return()
release the function without destroying it, so it can be returned from a function ...
G4double lowEnergyLimit
the energy per nucleon below which the MFP is zero
const G4String & GetParticleType() const
std::map< G4int, G4_c2_const_ptr > sigmaMap
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
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()
static const G4double A[nN]
const G4String & GetProcessName() const
virtual G4bool IsApplicable(const G4ParticleDefinition &aParticleType)
test if a prticle of type aParticleType can use this process
static G4IonTable * GetIonTable()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4int GetNumberOfIsotopes(G4int Z)
the exception class for c2_function operations.
void SetNIELPartitionFunction(const G4VNIELPartition *part)
set the pointer to a class for paritioning energy into NIEL
virtual G4ScreenedCoulombCrossSection * create()=0
G4double GetTotNbOfAtomsPerVolume() const
G4bool registerDepositedEnergy
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
void SetVerbosity(G4int v)
static log_log_interpolating_function_p< float_type > & log_log_interpolating_function()
make a *new object
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...
static c2_linear_p< float_type > & linear(float_type x0, float_type y0, float_type slope)
make a *new object
G4double standardmass(G4int z1)
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)
static const char * CVSFileVers()
const G4Material * targetMaterial
virtual ~G4ScreenedNuclearRecoil()
destructor
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
used internally by Geant4 machinery
G4VParticleChange * pParticleChange
a factory of pre-templated c2_function generators
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
void AddSecondary(G4Track *aSecondary)
virtual G4ScreenedCoulombCrossSection * GetNewCrossSectionHandler(void)
G4double GetEnergy() const
G4double GetHardeningFraction() const
get the fraction of particles which will have boosted scattering
size_t GetNumberOfElements() const
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)
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
class G4ParticleChange & GetParticleChange()
get the pointer to our ParticleChange object. for internal use, primarily.
G4ScreenedCoulombClassicalKinematics()
c2_linear_p< G4double > & xovereps
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()
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
virtual G4double PartitionNIEL(G4int z1, G4double a1, const G4Material *material, G4double energy) const =0
void unset_function()
clear our function
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()
static const double fermi
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
std::map< G4int, G4ScreenedCoulombCrossSection * > & GetCrossSectionHandlers()
static const double angstrom
Provides a factory class to avoid an infinite number of template.