46 fpMolWaterDensity = 0;
49 kineticEnergyCorrection[0]=0.;
50 kineticEnergyCorrection[1]=0.;
51 kineticEnergyCorrection[2]=0.;
52 kineticEnergyCorrection[3]=0.;
64 G4cout <<
"Miller & Green excitation model is constructed " <<
G4endl;
81 G4cout <<
"Calling G4DNAMillerGreenExcitationModel::Initialise()" <<
G4endl;
102 lowEnergyLimit[proton] = 10. * eV;
103 highEnergyLimit[proton] = 500. * keV;
105 kineticEnergyCorrection[0] = 1.;
106 slaterEffectiveCharge[0][0] = 0.;
107 slaterEffectiveCharge[1][0] = 0.;
108 slaterEffectiveCharge[2][0] = 0.;
109 sCoefficient[0][0] = 0.;
110 sCoefficient[1][0] = 0.;
111 sCoefficient[2][0] = 0.;
114 lowEnergyLimit[hydrogen] = 10. * eV;
115 highEnergyLimit[hydrogen] = 500. * keV;
117 kineticEnergyCorrection[0] = 1.;
118 slaterEffectiveCharge[0][0] = 0.;
119 slaterEffectiveCharge[1][0] = 0.;
120 slaterEffectiveCharge[2][0] = 0.;
121 sCoefficient[0][0] = 0.;
122 sCoefficient[1][0] = 0.;
123 sCoefficient[2][0] = 0.;
126 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
127 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
129 kineticEnergyCorrection[1] = 0.9382723/3.727417;
130 slaterEffectiveCharge[0][1]=0.;
131 slaterEffectiveCharge[1][1]=0.;
132 slaterEffectiveCharge[2][1]=0.;
133 sCoefficient[0][1]=0.;
134 sCoefficient[1][1]=0.;
135 sCoefficient[2][1]=0.;
138 lowEnergyLimit[alphaPlus] = 1. * keV;
139 highEnergyLimit[alphaPlus] = 400. * MeV;
141 kineticEnergyCorrection[2] = 0.9382723/3.727417;
142 slaterEffectiveCharge[0][2]=2.0;
145 slaterEffectiveCharge[1][2]=2.00;
146 slaterEffectiveCharge[2][2]=2.00;
148 sCoefficient[0][2]=0.7;
149 sCoefficient[1][2]=0.15;
150 sCoefficient[2][2]=0.15;
153 lowEnergyLimit[helium] = 1. * keV;
154 highEnergyLimit[helium] = 400. * MeV;
156 kineticEnergyCorrection[3] = 0.9382723/3.727417;
157 slaterEffectiveCharge[0][3]=1.7;
158 slaterEffectiveCharge[1][3]=1.15;
159 slaterEffectiveCharge[2][3]=1.15;
160 sCoefficient[0][3]=0.5;
161 sCoefficient[1][3]=0.25;
162 sCoefficient[2][3]=0.25;
166 if (particle==protonDef)
172 if (particle==hydrogenDef)
178 if (particle==alphaPlusPlusDef)
184 if (particle==alphaPlusDef)
190 if (particle==heliumDef)
203 G4cout <<
"Miller & Green excitation model is initialized " <<
G4endl
214 if (isInitialised) {
return; }
216 isInitialised =
true;
228 if (verboseLevel > 3)
229 G4cout <<
"Calling CrossSectionPerVolume() of G4DNAMillerGreenExcitationModel" <<
G4endl;
239 particleDefinition != instance->
GetIon(
"hydrogen")
241 particleDefinition != instance->
GetIon(
"alpha++")
243 particleDefinition != instance->
GetIon(
"alpha+")
245 particleDefinition != instance->
GetIon(
"helium")
256 if(waterDensity!= 0.0)
263 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
264 pos1 = lowEnergyLimit.find(particleName);
266 if (pos1 != lowEnergyLimit.end())
268 lowLim = pos1->second;
271 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
272 pos2 = highEnergyLimit.find(particleName);
274 if (pos2 != highEnergyLimit.end())
276 highLim = pos2->second;
279 if (k >= lowLim && k < highLim)
281 crossSection = Sum(k,particleDefinition);
334 if (verboseLevel > 2)
336 G4cout <<
"__________________________________" <<
G4endl;
337 G4cout <<
"°°° G4DNAMillerGreenExcitationModel - XS INFO START" <<
G4endl;
339 G4cout <<
"°°° Cross section per water molecule (cm^2)=" << crossSection/cm/cm <<
G4endl;
340 G4cout <<
"°°° Cross section per water molecule (cm^-1)=" << crossSection*waterDensity/(1./cm) <<
G4endl;
342 G4cout <<
"°°° G4DNAMillerGreenExcitationModel - XS INFO END" <<
G4endl;
347 if (verboseLevel > 2)
349 G4cout <<
"MillerGreenExcitationModel : WARNING Water density is NULL" <<
G4endl;
353 return crossSection*waterDensity;
367 if (verboseLevel > 3)
368 G4cout <<
"Calling SampleSecondaries() of G4DNAMillerGreenExcitationModel" <<
G4endl;
377 const G4double excitation[]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
378 G4double excitationEnergy = excitation[level];
380 G4double newEnergy = particleEnergy0 - excitationEnergy;
416 const G4double sigma0(1.E+8 * barn);
418 const G4double aj[]={876.*eV, 2084.* eV, 1373.*eV, 692.*eV, 900.*eV};
419 const G4double jj[]={19820.*eV, 23490.*eV, 27770.*eV, 30830.*eV, 33080.*eV};
420 const G4double omegaj[]={0.85, 0.88, 0.88, 0.78, 0.78};
423 const G4double Eliq[5]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
425 G4int particleTypeIndex = 0;
430 if (particleDefinition == instance->
GetIon(
"hydrogen")) particleTypeIndex=0;
431 if (particleDefinition == instance->
GetIon(
"alpha++")) particleTypeIndex=1;
432 if (particleDefinition == instance->
GetIon(
"alpha+")) particleTypeIndex=2;
433 if (particleDefinition == instance->
GetIon(
"helium")) particleTypeIndex=3;
436 tCorrected = k * kineticEnergyCorrection[particleTypeIndex];
439 if (tCorrected < Eliq[excitationLevel])
return 0;
445 numerator = std::pow(z * aj[excitationLevel], omegaj[excitationLevel]) *
446 std::pow(tCorrected - Eliq[excitationLevel], nu);
450 if (particleDefinition == instance->
GetIon(
"hydrogen"))
451 numerator = std::pow(z * 0.75*aj[excitationLevel], omegaj[excitationLevel]) *
452 std::pow(tCorrected - Eliq[excitationLevel], nu);
456 power = omegaj[excitationLevel] + nu;
459 denominator = std::pow(jj[excitationLevel], power) + std::pow(tCorrected, power);
463 zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, Eliq[excitationLevel], slaterEffectiveCharge[0][particleTypeIndex], 1.) +
464 sCoefficient[1][particleTypeIndex] * S_2s(k, Eliq[excitationLevel], slaterEffectiveCharge[1][particleTypeIndex], 2.) +
465 sCoefficient[2][particleTypeIndex] * S_2p(k, Eliq[excitationLevel], slaterEffectiveCharge[2][particleTypeIndex], 2.) );
467 if (particleDefinition == instance->
GetIon(
"hydrogen")) zEff = 1.;
469 G4double cross = sigma0 * zEff * zEff * numerator / denominator;
481 std::deque<double> values;
486 if ( particle == instance->
GetIon(
"alpha++") ||
488 particle == instance->
GetIon(
"hydrogen") ||
489 particle == instance->
GetIon(
"alpha+") ||
490 particle == instance->
GetIon(
"helium")
496 G4double partial = PartialCrossSection(k,i,particle);
497 values.push_front(partial);
508 if (values[i] > value)
return i;
565 for (
G4int i=0; i<nLevels; i++)
567 totalCrossSection += PartialCrossSection(k,i,particle);
569 return totalCrossSection;
582 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
583 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
599 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
600 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
616 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
617 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
632 G4double tElectron = 0.511/3728. * t;
636 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (_slaterEffectiveCharge/shellNumber);
G4DLLIMPORT std::ostream G4cout
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
static G4DNAGenericIonsManager * Instance(void)
G4ParticleDefinition * GetIon(const G4String &name)
G4DNAMillerGreenExcitationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAMillerGreenExcitationModel")
G4ParticleChangeForGamma * fParticleChangeForGamma
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 void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual ~G4DNAMillerGreenExcitationModel()
static G4DNAMolecularMaterial * Instance()
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
const G4Track * GetCurrentTrack() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4int GetLeptonNumber() const
G4double GetPDGCharge() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
void SetHighEnergyLimit(G4double)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void SetLowEnergyLimit(G4double)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)