47 killBelowEnergy = 9*
eV;
48 lowEnergyLimit = 0 *
eV;
49 intermediateEnergyLimit = 200 *
eV;
50 highEnergyLimit = 1. *
MeV;
51 SetLowEnergyLimit(lowEnergyLimit);
52 SetHighEnergyLimit(highEnergyLimit);
64 G4cout <<
"Screened Rutherford Elastic model is constructed " <<
G4endl
66 << lowEnergyLimit /
eV <<
" eV - "
67 << highEnergyLimit /
MeV <<
" MeV"
70 fParticleChangeForGamma = 0;
85 G4cout <<
"Calling G4DNAScreenedRutherfordElasticModel::Initialise()" <<
G4endl;
89 if (LowEnergyLimit() < lowEnergyLimit)
91 G4cout <<
"G4DNAScreenedRutherfordElasticModel: low energy limit increased from " <<
92 LowEnergyLimit()/
eV <<
" eV to " << lowEnergyLimit/
eV <<
" eV" <<
G4endl;
93 SetLowEnergyLimit(lowEnergyLimit);
96 if (HighEnergyLimit() > highEnergyLimit)
98 G4cout <<
"G4DNAScreenedRutherfordElasticModel: high energy limit decreased from " <<
99 HighEnergyLimit()/
MeV <<
" MeV to " << highEnergyLimit/
MeV <<
" MeV" <<
G4endl;
100 SetHighEnergyLimit(highEnergyLimit);
105 betaCoeff.push_back(7.51525);
106 betaCoeff.push_back(-0.41912);
107 betaCoeff.push_back(7.2017E-3);
108 betaCoeff.push_back(-4.646E-5);
109 betaCoeff.push_back(1.02897E-7);
111 deltaCoeff.push_back(2.9612);
112 deltaCoeff.push_back(-0.26376);
113 deltaCoeff.push_back(4.307E-3);
114 deltaCoeff.push_back(-2.6895E-5);
115 deltaCoeff.push_back(5.83505E-8);
117 gamma035_10Coeff.push_back(-1.7013);
118 gamma035_10Coeff.push_back(-1.48284);
119 gamma035_10Coeff.push_back(0.6331);
120 gamma035_10Coeff.push_back(-0.10911);
121 gamma035_10Coeff.push_back(8.358E-3);
122 gamma035_10Coeff.push_back(-2.388E-4);
124 gamma10_100Coeff.push_back(-3.32517);
125 gamma10_100Coeff.push_back(0.10996);
126 gamma10_100Coeff.push_back(-4.5255E-3);
127 gamma10_100Coeff.push_back(5.8372E-5);
128 gamma10_100Coeff.push_back(-2.4659E-7);
130 gamma100_200Coeff.push_back(2.4775E-2);
131 gamma100_200Coeff.push_back(-2.96264E-5);
132 gamma100_200Coeff.push_back(-1.20655E-7);
138 G4cout <<
"Screened Rutherford elastic model is initialized " << G4endl
140 << LowEnergyLimit() /
eV <<
" eV - "
141 << HighEnergyLimit() /
MeV <<
" MeV"
148 if (isInitialised) {
return; }
149 fParticleChangeForGamma = GetParticleChangeForGamma();
150 isInitialised =
true;
162 if (verboseLevel > 3)
163 G4cout <<
"Calling CrossSectionPerVolume() of G4DNAScreenedRutherfordElasticModel" <<
G4endl;
171 if(waterDensity!= 0.0)
175 if (ekin < highEnergyLimit)
178 if (ekin < killBelowEnergy)
return DBL_MAX;
182 G4double crossSection = RutherfordCrossSection(ekin, z);
183 sigma =
pi * crossSection / (n * (n + 1.));
186 if (verboseLevel > 2)
188 G4cout <<
"__________________________________" <<
G4endl;
189 G4cout <<
"°°° G4DNAScreenedRutherfordElasticModel - XS INFO START" <<
G4endl;
191 G4cout <<
"°°° Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
192 G4cout <<
"°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
194 G4cout <<
"°°° G4DNAScreenedRutherfordElasticModel - XS INFO END" <<
G4endl;
216 G4double cross = z * ( z + 1) * length * length;
241 G4double numerator = (alpha_1 + beta_1 * std::log(k/
eV)) * constK * std::pow(z, 2./3.);
248 if (denominator > 0.) value = numerator / denominator;
262 if (verboseLevel > 3)
263 G4cout <<
"Calling SampleSecondaries() of G4DNAScreenedRutherfordElasticModel" <<
G4endl;
267 if (electronEnergy0 < killBelowEnergy)
269 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
270 fParticleChangeForGamma->ProposeTrackStatus(
fStopAndKill);
271 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
277 if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
279 if (electronEnergy0<intermediateEnergyLimit)
281 if (verboseLevel > 3)
G4cout <<
"---> Using Brenner & Zaider model" <<
G4endl;
282 cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0);
285 if (electronEnergy0>=intermediateEnergyLimit)
287 if (verboseLevel > 3)
G4cout <<
"---> Using Screened Rutherford model" <<
G4endl;
289 cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z);
298 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
300 xDir *= std::cos(phi);
301 yDir *= std::sin(phi);
303 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
305 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.
unit()) ;
307 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
314 G4double G4DNAScreenedRutherfordElasticModel::BrennerZaiderRandomizeCosTheta(
G4double k)
328 G4double beta = std::exp(CalculatePolynomial(k,betaCoeff));
329 G4double delta = std::exp(CalculatePolynomial(k,deltaCoeff));
334 gamma = CalculatePolynomial(k, gamma100_200Coeff);
341 gamma = std::exp(CalculatePolynomial(k, gamma10_100Coeff));
345 gamma = std::exp(CalculatePolynomial(k, gamma035_10Coeff));
351 G4double oneOverMax = 1. / (1./(4.*gamma*gamma) + beta/( (2.+2.*delta)*(2.+2.*delta) ));
362 leftDenominator = (1. + 2.*gamma - cosTheta);
363 rightDenominator = (1. + 2.*delta + cosTheta);
364 if ( (leftDenominator * rightDenominator) != 0. )
366 fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
416 G4double G4DNAScreenedRutherfordElasticModel::CalculatePolynomial(
G4double k, std::vector<G4double>& vec)
423 size_t size = vec.size();
463 fCosTheta = (1 + 2.*n - cosTheta);
464 if (fCosTheta !=0.) fCosTheta = oneOverMax / (fCosTheta*fCosTheta);
G4double GetKineticEnergy() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
const G4ThreeVector & GetMomentumDirection() const
const G4double * GetAtomicNumDensityVector() const
static G4DNAMolecularMaterial * Instance()
Hep3Vector orthogonal() const
const XML_Char int const XML_Char * value
G4DNAScreenedRutherfordElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAScreenedRutherfordElasticModel")
Hep3Vector cross(const Hep3Vector &) const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual ~G4DNAScreenedRutherfordElasticModel()