43 fpMolWaterDensity = 0;
45 numberOfPartialCrossSections[0] = 0;
46 numberOfPartialCrossSections[1] = 0;
58 G4cout <<
"Dingfelder charge increase model is constructed " <<
G4endl;
80 G4cout <<
"Calling G4DNADingfelderChargeIncreaseModel::Initialise()"
99 lowEnergyLimit[hydrogen] = 100. * eV;
100 highEnergyLimit[hydrogen] = 100. * MeV;
103 lowEnergyLimit[alphaPlus] = 1. * keV;
104 highEnergyLimit[alphaPlus] = 400. * MeV;
107 lowEnergyLimit[helium] = 1. * keV;
108 highEnergyLimit[helium] = 400. * MeV;
112 if (particle==hydrogenDef)
118 if (particle==alphaPlusDef)
124 if (particle==heliumDef)
145 numberOfPartialCrossSections[0]=1;
171 numberOfPartialCrossSections[1]=2;
177 G4cout <<
"Dingfelder charge increase model is initialized " <<
G4endl
188 if (isInitialised)
return;
191 isInitialised =
true;
202 if (verboseLevel > 3)
205 <<
"Calling CrossSectionPerVolume() of G4DNADingfelderChargeIncreaseModel"
215 particleDefinition != instance->
GetIon(
"hydrogen")
217 particleDefinition != instance->
GetIon(
"alpha+")
219 particleDefinition != instance->
GetIon(
"helium")
232 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
233 pos1 = lowEnergyLimit.find(particleName);
235 if (pos1 != lowEnergyLimit.end())
237 lowLim = pos1->second;
240 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
241 pos2 = highEnergyLimit.find(particleName);
243 if (pos2 != highEnergyLimit.end())
245 highLim = pos2->second;
248 if (k >= lowLim && k <= highLim)
251 if (particleDefinition == instance->
GetIon(
"hydrogen"))
260 G4double t = k / (proton_mass_c2/electron_mass_c2);
262 G4double temp = 4.0 * pi * Bohr_radius/nm * Bohr_radius/nm * fac;
263 G4double sigmal = temp * cc * (std::pow(x,dd));
264 G4double sigmah = temp * (aa * std::log(1.0 + x) + bb) / x;
265 totalCrossSection = 1.0/(1.0/sigmal + 1.0/sigmah) *m*m;
269 totalCrossSection = Sum(k,particleDefinition);
273 if (verboseLevel > 2)
275 G4cout <<
"__________________________________" <<
G4endl;
276 G4cout <<
"G4DNADingfelderChargeIncreaseModel - XS INFO START" <<
G4endl;
277 G4cout <<
"Kinetic energy(eV)=" << k/eV <<
" particle : " << particleName <<
G4endl;
278 G4cout <<
"Cross section per water molecule (cm^2)=" << totalCrossSection/cm/cm <<
G4endl;
279 G4cout <<
"Cross section per water molecule (cm^-1)=" << totalCrossSection*waterDensity/(1./cm) <<
G4endl;
282 G4cout <<
"G4DNADingfelderChargeIncreaseModel - XS INFO END" <<
G4endl;
285 return totalCrossSection*waterDensity;
299 if (verboseLevel > 3)
302 <<
"Calling SampleSecondaries() of G4DNADingfelderChargeIncreaseModel"
314 G4int finalStateIndex = RandomSelect(inK,definition);
316 G4int n = NumberOfFinalStates(definition,finalStateIndex);
320 if (!statCode) outK = inK - IncomingParticleBindingEnergyConstant(definition,finalStateIndex);
326 ProposeLocalEnergyDeposit(IncomingParticleBindingEnergyConstant(definition,finalStateIndex));
334 if (definition == instance->
GetIon(
"hydrogen")) electronK = inK*electron_mass_c2/proton_mass_c2;
335 else electronK = inK*electron_mass_c2/(particleMass);
339 G4Exception(
"G4DNADingfelderChargeIncreaseModel::SampleSecondaries",
"em0004",
347 fvect->push_back(dp);
362 G4int finalStateIndex)
368 if (particleDefinition == instance->
GetIon(
"hydrogen"))
371 if (particleDefinition == instance->
GetIon(
"alpha+"))
374 if (particleDefinition == instance->
GetIon(
"helium"))
376 if (finalStateIndex == 0)
387 G4int finalStateIndex)
391 if (particleDefinition == instance->
GetIon(
"hydrogen"))
394 if (particleDefinition == instance->
GetIon(
"alpha+"))
395 return instance->
GetIon(
"alpha++");
397 if (particleDefinition == instance->
GetIon(
"helium"))
399 if (finalStateIndex == 0)
400 return instance->
GetIon(
"alpha+");
401 return instance->
GetIon(
"alpha++");
410 G4int finalStateIndex)
414 if (particleDefinition == instance->
GetIon(
"hydrogen"))
417 if (particleDefinition == instance->
GetIon(
"alpha+"))
424 if (particleDefinition == instance->
GetIon(
"helium"))
429 if (finalStateIndex == 0)
431 return (54.509 + 24.587) * eV;
443 G4int particleTypeIndex = 0;
447 if (particleDefinition == instance->
GetIon(
"alpha+"))
448 particleTypeIndex = 0;
450 if (particleDefinition == instance->
GetIon(
"helium"))
451 particleTypeIndex = 1;
472 if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex])
483 x1[index][particleTypeIndex] = x0[index][particleTypeIndex]
484 + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
485 / (c0[index][particleTypeIndex]
486 * d0[index][particleTypeIndex]),
487 1. / (d0[index][particleTypeIndex] - 1.));
488 b1[index][particleTypeIndex] = (a0[index][particleTypeIndex]
489 - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
490 + b0[index][particleTypeIndex]
491 - c0[index][particleTypeIndex]
492 * std::pow(x1[index][particleTypeIndex]
493 - x0[index][particleTypeIndex],
494 d0[index][particleTypeIndex]);
500 if (x < x0[index][particleTypeIndex])
501 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
502 else if (x < x1[index][particleTypeIndex])
503 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]
504 - c0[index][particleTypeIndex]
505 * std::pow(x - x0[index][particleTypeIndex],
506 d0[index][particleTypeIndex]);
508 y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
510 return f0[index][particleTypeIndex] * std::pow(10., y) * m * m;
516G4int G4DNADingfelderChargeIncreaseModel::RandomSelect(
G4double k,
519 G4int particleTypeIndex = 0;
523 if (particleDefinition == instance->
GetIon(
"hydrogen"))
526 if (particleDefinition == instance->
GetIon(
"alpha+"))
527 particleTypeIndex = 0;
529 if (particleDefinition == instance->
GetIon(
"helium"))
530 particleTypeIndex = 1;
532 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
540 values[i] = PartialCrossSection(k, i, particleDefinition);
551 if (values[i] > value)
567 G4int particleTypeIndex = 0;
571 if (particleDefinition == instance->
GetIon(
"alpha+"))
572 particleTypeIndex = 0;
574 if (particleDefinition == instance->
GetIon(
"helium"))
575 particleTypeIndex = 1;
579 for (
G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++)
581 totalCrossSection += PartialCrossSection(k, i, particleDefinition);
584 return totalCrossSection;
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4GLOB_DLL std::ostream G4cout
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual ~G4DNADingfelderChargeIncreaseModel()
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4DNADingfelderChargeIncreaseModel(const G4ParticleDefinition *p=0, const G4String &nam="DNADingfelderChargeIncreaseModel")
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)
static G4DNAGenericIonsManager * Instance(void)
G4ParticleDefinition * GetIon(const G4String &name)
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
static G4Electron * Electron()
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
G4double GetPDGMass() const
const G4String & GetParticleName() const
static G4Proton * Proton()
void SetHighEnergyLimit(G4double)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void SetLowEnergyLimit(G4double)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)