47 fpMolWaterDensity =
nullptr;
50 kineticEnergyCorrection[0]=0.;
51 kineticEnergyCorrection[1]=0.;
52 kineticEnergyCorrection[2]=0.;
53 kineticEnergyCorrection[3]=0.;
65 G4cout <<
"Miller & Green excitation model is constructed " <<
G4endl;
86 G4cout <<
"Calling G4DNAMillerGreenExcitationModel::Initialise()" <<
G4endl;
93 hydrogenDef = instance->
GetIon(
"hydrogen");
95 alphaPlusDef = instance->
GetIon(
"alpha+");
96 heliumDef = instance->
GetIon(
"helium");
107 lowEnergyLimit[proton] = 10. * eV;
108 highEnergyLimit[proton] = 500. * keV;
110 kineticEnergyCorrection[0] = 1.;
111 slaterEffectiveCharge[0][0] = 0.;
112 slaterEffectiveCharge[1][0] = 0.;
113 slaterEffectiveCharge[2][0] = 0.;
114 sCoefficient[0][0] = 0.;
115 sCoefficient[1][0] = 0.;
116 sCoefficient[2][0] = 0.;
119 lowEnergyLimit[hydrogen] = 10. * eV;
120 highEnergyLimit[hydrogen] = 500. * keV;
122 kineticEnergyCorrection[0] = 1.;
123 slaterEffectiveCharge[0][0] = 0.;
124 slaterEffectiveCharge[1][0] = 0.;
125 slaterEffectiveCharge[2][0] = 0.;
126 sCoefficient[0][0] = 0.;
127 sCoefficient[1][0] = 0.;
128 sCoefficient[2][0] = 0.;
131 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
132 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
134 kineticEnergyCorrection[1] = 0.9382723/3.727417;
135 slaterEffectiveCharge[0][1]=0.;
136 slaterEffectiveCharge[1][1]=0.;
137 slaterEffectiveCharge[2][1]=0.;
138 sCoefficient[0][1]=0.;
139 sCoefficient[1][1]=0.;
140 sCoefficient[2][1]=0.;
143 lowEnergyLimit[alphaPlus] = 1. * keV;
144 highEnergyLimit[alphaPlus] = 400. * MeV;
146 kineticEnergyCorrection[2] = 0.9382723/3.727417;
147 slaterEffectiveCharge[0][2]=2.0;
150 slaterEffectiveCharge[1][2]=2.00;
151 slaterEffectiveCharge[2][2]=2.00;
153 sCoefficient[0][2]=0.7;
154 sCoefficient[1][2]=0.15;
155 sCoefficient[2][2]=0.15;
158 lowEnergyLimit[helium] = 1. * keV;
159 highEnergyLimit[helium] = 400. * MeV;
161 kineticEnergyCorrection[3] = 0.9382723/3.727417;
162 slaterEffectiveCharge[0][3]=1.7;
163 slaterEffectiveCharge[1][3]=1.15;
164 slaterEffectiveCharge[2][3]=1.15;
165 sCoefficient[0][3]=0.5;
166 sCoefficient[1][3]=0.25;
167 sCoefficient[2][3]=0.25;
171 if (particle==protonDef)
177 if (particle==hydrogenDef)
183 if (particle==alphaPlusPlusDef)
189 if (particle==alphaPlusDef)
195 if (particle==heliumDef)
208 G4cout <<
"Miller & Green excitation model is initialized " <<
G4endl
219 if (isInitialised) {
return; }
221 isInitialised =
true;
233 if (verboseLevel > 3)
234 G4cout <<
"Calling CrossSectionPerVolume() of G4DNAMillerGreenExcitationModel" <<
G4endl;
239 particleDefinition != protonDef
241 particleDefinition != hydrogenDef
243 particleDefinition != alphaPlusPlusDef
245 particleDefinition != alphaPlusDef
247 particleDefinition != heliumDef
260 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
261 pos1 = lowEnergyLimit.find(particleName);
263 if (pos1 != lowEnergyLimit.end())
265 lowLim = pos1->second;
268 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
269 pos2 = highEnergyLimit.find(particleName);
271 if (pos2 != highEnergyLimit.end())
273 highLim = pos2->second;
276 if (k >= lowLim && k <= highLim)
278 crossSection = Sum(k,particleDefinition);
331 if (verboseLevel > 2)
333 G4cout <<
"__________________________________" <<
G4endl;
334 G4cout <<
"G4DNAMillerGreenExcitationModel - XS INFO START" <<
G4endl;
336 G4cout <<
"Cross section per water molecule (cm^2)=" << crossSection/cm/cm <<
G4endl;
337 G4cout <<
"Cross section per water molecule (cm^-1)=" << crossSection*waterDensity/(1./cm) <<
G4endl;
339 G4cout <<
"G4DNAMillerGreenExcitationModel - XS INFO END" <<
G4endl;
342 return crossSection*waterDensity;
354 if (verboseLevel > 3)
355 G4cout <<
"Calling SampleSecondaries() of G4DNAMillerGreenExcitationModel" <<
G4endl;
362 const G4double excitation[]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
363 G4double excitationEnergy = excitation[level];
367 if (!statCode) newEnergy = particleEnergy0 - excitationEnergy;
369 else newEnergy = particleEnergy0;
379 level, theIncomingTrack);
392 return PartialCrossSection(kineticEnergy, level, particleDefinition);
414 const G4double sigma0(1.E+8 * barn);
416 const G4double aj[]={876.*eV, 2084.* eV, 1373.*eV, 692.*eV, 900.*eV};
417 const G4double jj[]={19820.*eV, 23490.*eV, 27770.*eV, 30830.*eV, 33080.*eV};
418 const G4double omegaj[]={0.85, 0.88, 0.88, 0.78, 0.78};
421 const G4double Eliq[5]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
423 G4int particleTypeIndex = 0;
425 if (particleDefinition == protonDef) particleTypeIndex=0;
426 if (particleDefinition == hydrogenDef) particleTypeIndex=0;
427 if (particleDefinition == alphaPlusPlusDef) particleTypeIndex=1;
428 if (particleDefinition == alphaPlusDef) particleTypeIndex=2;
429 if (particleDefinition == heliumDef) particleTypeIndex=3;
432 tCorrected = k * kineticEnergyCorrection[particleTypeIndex];
435 if (tCorrected < Eliq[excitationLevel])
return 0;
441 numerator = gpow->
powA(z * aj[excitationLevel], omegaj[excitationLevel]) *
442 gpow->
powA(tCorrected - Eliq[excitationLevel], nu);
446 if (particleDefinition == hydrogenDef)
447 numerator = gpow->
powA(z * 0.75*aj[excitationLevel], omegaj[excitationLevel]) *
448 gpow->
powA(tCorrected - Eliq[excitationLevel], nu);
452 power = omegaj[excitationLevel] + nu;
455 denominator = gpow->
powA(jj[excitationLevel], power) + gpow->
powA(tCorrected, power);
459 zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, Eliq[excitationLevel], slaterEffectiveCharge[0][particleTypeIndex], 1.) +
460 sCoefficient[1][particleTypeIndex] * S_2s(k, Eliq[excitationLevel], slaterEffectiveCharge[1][particleTypeIndex], 2.) +
461 sCoefficient[2][particleTypeIndex] * S_2p(k, Eliq[excitationLevel], slaterEffectiveCharge[2][particleTypeIndex], 2.) );
463 if (particleDefinition == hydrogenDef) zEff = 1.;
465 G4double cross = sigma0 * zEff * zEff * numerator / denominator;
477 std::deque<G4double> values;
479 if ( particle == alphaPlusPlusDef ||
480 particle == protonDef||
481 particle == hydrogenDef ||
482 particle == alphaPlusDef ||
483 particle == heliumDef
489 G4double partial = PartialCrossSection(k,i,particle);
490 values.push_front(partial);
501 if (values[i] > value)
return i;
558 for (
G4int i=0; i<nLevels; i++)
560 totalCrossSection += PartialCrossSection(k,i,particle);
562 return totalCrossSection;
575 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
576 G4double value = 1. -
G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
592 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
593 G4double value = 1. -
G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
609 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
610 G4double value = 1. -
G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
625 G4double tElectron = 0.511/3728. * t;
629 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (_slaterEffectiveCharge/shellNumber);
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4GLOB_DLL std::ostream G4cout
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
static G4DNAGenericIonsManager * Instance()
G4ParticleDefinition * GetIon(const G4String &name)
~G4DNAMillerGreenExcitationModel() override
G4ParticleChangeForGamma * fParticleChangeForGamma
G4double GetPartialCrossSection(const G4Material *, G4int, const G4ParticleDefinition *, G4double) override
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4DNAMillerGreenExcitationModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="DNAMillerGreenExcitationModel")
G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax) override
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
const std::vector< G4double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
static G4DNAMolecularMaterial * Instance()
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
std::size_t GetIndex() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4int GetLeptonNumber() const
G4double GetPDGCharge() const
const G4String & GetParticleName() const
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
static G4Proton * ProtonDefinition()
void SetHighEnergyLimit(G4double)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void SetLowEnergyLimit(G4double)
const G4Track * GetCurrentTrack() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)