45 fpMolWaterDensity = 0;
47 numberOfPartialCrossSections[0] = 0;
48 numberOfPartialCrossSections[1] = 0;
60 G4cout <<
"Dingfelder charge increase model is constructed " <<
G4endl;
82 G4cout <<
"Calling G4DNADingfelderChargeIncreaseModel::Initialise()"
101 lowEnergyLimit[hydrogen] = 100. *
eV;
102 highEnergyLimit[hydrogen] = 100. *
MeV;
105 lowEnergyLimit[alphaPlus] = 1. *
keV;
106 highEnergyLimit[alphaPlus] = 400. *
MeV;
109 lowEnergyLimit[helium] = 1. *
keV;
110 highEnergyLimit[helium] = 400. *
MeV;
114 if (particle==hydrogenDef)
120 if (particle==alphaPlusDef)
126 if (particle==heliumDef)
147 numberOfPartialCrossSections[0]=1;
173 numberOfPartialCrossSections[1]=2;
179 G4cout <<
"Dingfelder charge increase model is initialized " <<
G4endl
193 isInitialised =
true;
204 if (verboseLevel > 3)
207 <<
"Calling CrossSectionPerVolume() of G4DNADingfelderChargeIncreaseModel"
217 particleDefinition != instance->
GetIon(
"hydrogen")
219 particleDefinition != instance->
GetIon(
"alpha+")
221 particleDefinition != instance->
GetIon(
"helium")
232 if(waterDensity!= 0.0)
237 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
238 pos1 = lowEnergyLimit.find(particleName);
240 if (pos1 != lowEnergyLimit.end())
242 lowLim = pos1->second;
245 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
246 pos2 = highEnergyLimit.find(particleName);
248 if (pos2 != highEnergyLimit.end())
250 highLim = pos2->second;
253 if (k >= lowLim && k < highLim)
256 if (particleDefinition == instance->
GetIon(
"hydrogen"))
268 G4double sigmal = temp * cc * (std::pow(x,dd));
269 G4double sigmah = temp * (aa * std::log(1.0 + x) +
bb) / x;
270 totalCrossSection = 1.0/(1.0/sigmal + 1.0/sigmah) *
m*
m;
274 totalCrossSection = Sum(k,particleDefinition);
278 if (verboseLevel > 2)
280 G4cout <<
"__________________________________" <<
G4endl;
281 G4cout <<
"G4DNADingfelderChargeIncreaseModel - XS INFO START" <<
G4endl;
282 G4cout <<
"Kinetic energy(eV)=" << k/
eV <<
" particle : " << particleName <<
G4endl;
283 G4cout <<
"Cross section per water molecule (cm^2)=" << totalCrossSection/
cm/
cm <<
G4endl;
284 G4cout <<
"Cross section per water molecule (cm^-1)=" << totalCrossSection*waterDensity/(1./
cm) << G4endl;
286 G4cout <<
"G4DNADingfelderChargeIncreaseModel - XS INFO END" <<
G4endl;
291 return totalCrossSection*waterDensity;
305 if (verboseLevel > 3)
308 <<
"Calling SampleSecondaries() of G4DNADingfelderChargeIncreaseModel"
320 G4int finalStateIndex = RandomSelect(inK,definition);
322 G4int n = NumberOfFinalStates(definition,finalStateIndex);
326 if (!statCode) outK = inK - IncomingParticleBindingEnergyConstant(definition,finalStateIndex);
344 G4Exception(
"G4DNADingfelderChargeIncreaseModel::SampleSecondaries",
"em0004",
352 fvect->push_back(dp);
367 G4int finalStateIndex)
372 if (particleDefinition == instance->
GetIon(
"hydrogen"))
374 if (particleDefinition == instance->
GetIon(
"alpha+"))
377 if (particleDefinition == instance->
GetIon(
"helium"))
379 if (finalStateIndex == 0)
389 G4int finalStateIndex)
392 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++");
409 G4int finalStateIndex)
412 if (particleDefinition == instance->
GetIon(
"hydrogen"))
415 if (particleDefinition == instance->
GetIon(
"alpha+"))
422 if (particleDefinition == instance->
GetIon(
"helium"))
427 if (finalStateIndex == 0)
429 return (54.509 + 24.587) *
eV;
439 G4double G4DNADingfelderChargeIncreaseModel::PartialCrossSection(
G4double k,
443 G4int particleTypeIndex = 0;
447 if (particleDefinition == instance->
GetIon(
"alpha+"))
448 particleTypeIndex = 0;
449 if (particleDefinition == instance->
GetIon(
"helium"))
450 particleTypeIndex = 1;
470 if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex])
480 x1[index][particleTypeIndex] = x0[index][particleTypeIndex]
481 + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
482 / (c0[index][particleTypeIndex]
483 * d0[index][particleTypeIndex]),
484 1. / (d0[index][particleTypeIndex] - 1.));
485 b1[index][particleTypeIndex] = (a0[index][particleTypeIndex]
486 - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
487 + b0[index][particleTypeIndex]
488 - c0[index][particleTypeIndex]
489 * std::pow(x1[index][particleTypeIndex]
490 - x0[index][particleTypeIndex],
491 d0[index][particleTypeIndex]);
497 if (x < x0[index][particleTypeIndex])
498 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
499 else if (x < x1[index][particleTypeIndex])
500 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]
501 - c0[index][particleTypeIndex]
502 * std::pow(x - x0[index][particleTypeIndex],
503 d0[index][particleTypeIndex]);
505 y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
507 return f0[index][particleTypeIndex] * std::pow(10., y) *
m *
m;
513 G4int G4DNADingfelderChargeIncreaseModel::RandomSelect(
G4double k,
516 G4int particleTypeIndex = 0;
520 if (particleDefinition == instance->
GetIon(
"hydrogen"))
522 if (particleDefinition == instance->
GetIon(
"alpha+"))
523 particleTypeIndex = 0;
524 if (particleDefinition == instance->
GetIon(
"helium"))
525 particleTypeIndex = 1;
527 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
535 values[i] = PartialCrossSection(k, i, particleDefinition);
546 if (values[i] > value)
562 G4int particleTypeIndex = 0;
566 if (particleDefinition == instance->
GetIon(
"alpha+"))
567 particleTypeIndex = 0;
568 if (particleDefinition == instance->
GetIon(
"helium"))
569 particleTypeIndex = 1;
573 for (
G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++)
575 totalCrossSection += PartialCrossSection(k, i, particleDefinition);
577 return totalCrossSection;
G4double LowEnergyLimit() const
static constexpr double proton_mass_c2
G4double GetKineticEnergy() const
static constexpr double Bohr_radius
G4double HighEnergyLimit() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
virtual ~G4DNADingfelderChargeIncreaseModel()
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4ParticleDefinition * GetDefinition() const
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double electron_mass_c2
void SetHighEnergyLimit(G4double)
G4GLOB_DLL std::ostream G4cout
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
static constexpr double m
const XML_Char int const XML_Char * value
const G4ThreeVector & GetMomentumDirection() const
static constexpr double cm
static G4Proton * Proton()
static constexpr double eV
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
static G4DNAGenericIonsManager * Instance(void)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static G4DNAMolecularMaterial * Instance()
G4double GetPDGMass() const
static constexpr double nm
G4ParticleChangeForGamma * fParticleChangeForGamma
static G4Electron * Electron()
static const G4double fac
static constexpr double MeV
static constexpr double pi
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
G4DNADingfelderChargeIncreaseModel(const G4ParticleDefinition *p=0, const G4String &nam="DNADingfelderChargeIncreaseModel")
static constexpr double keV
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4ParticleDefinition * GetIon(const G4String &name)