45 fpMolWaterDensity = 0;
47 numberOfPartialCrossSections[0]=0;
48 numberOfPartialCrossSections[1]=0;
60 G4cout <<
"Dingfelder charge increase model is constructed " <<
G4endl;
77 G4cout <<
"Calling G4DNADingfelderChargeIncreaseModel::Initialise()" <<
G4endl;
94 lowEnergyLimit[hydrogen] = 100. * eV;
95 highEnergyLimit[hydrogen] = 100. * MeV;
98 lowEnergyLimit[alphaPlus] = 1. * keV;
99 highEnergyLimit[alphaPlus] = 400. * MeV;
102 lowEnergyLimit[helium] = 1. * keV;
103 highEnergyLimit[helium] = 400. * MeV;
107 if (particle==hydrogenDef)
113 if (particle==alphaPlusDef)
119 if (particle==heliumDef)
140 numberOfPartialCrossSections[0]=1;
166 numberOfPartialCrossSections[1]=2;
172 G4cout <<
"Dingfelder charge increase model is initialized " <<
G4endl
183 if (isInitialised) {
return; }
185 isInitialised =
true;
196 if (verboseLevel > 3)
197 G4cout <<
"Calling CrossSectionPerVolume() of G4DNADingfelderChargeIncreaseModel" <<
G4endl;
205 particleDefinition != instance->
GetIon(
"hydrogen")
207 particleDefinition != instance->
GetIon(
"alpha+")
209 particleDefinition != instance->
GetIon(
"helium")
220 if(waterDensity!= 0.0)
225 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
226 pos1 = lowEnergyLimit.find(particleName);
228 if (pos1 != lowEnergyLimit.end())
230 lowLim = pos1->second;
233 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
234 pos2 = highEnergyLimit.find(particleName);
236 if (pos2 != highEnergyLimit.end())
238 highLim = pos2->second;
241 if (k >= lowLim && k < highLim)
244 if (particleDefinition == instance->
GetIon(
"hydrogen"))
253 G4double t = k / (proton_mass_c2/electron_mass_c2);
255 G4double temp = 4.0 * pi * Bohr_radius/nm * Bohr_radius/nm * fac;
256 G4double sigmal = temp * cc * (std::pow(x,dd));
257 G4double sigmah = temp * (aa * std::log(1.0 + x) + bb) / x;
258 totalCrossSection = 1.0/(1.0/sigmal + 1.0/sigmah) *m*m;
262 totalCrossSection = Sum(k,particleDefinition);
266 if (verboseLevel > 2)
268 G4cout <<
"__________________________________" <<
G4endl;
269 G4cout <<
"°°° G4DNADingfelderChargeIncreaseModel - XS INFO START" <<
G4endl;
270 G4cout <<
"°°° Kinetic energy(eV)=" << k/eV <<
" particle : " << particleName <<
G4endl;
271 G4cout <<
"°°° Cross section per water molecule (cm^2)=" << totalCrossSection/cm/cm <<
G4endl;
272 G4cout <<
"°°° Cross section per water molecule (cm^-1)=" << totalCrossSection*waterDensity/(1./cm) <<
G4endl;
274 G4cout <<
"°°° G4DNADingfelderChargeIncreaseModel - XS INFO END" <<
G4endl;
279 return totalCrossSection*waterDensity;
292 if (verboseLevel > 3)
293 G4cout <<
"Calling SampleSecondaries() of G4DNADingfelderChargeIncreaseModel" <<
G4endl;
304 G4int finalStateIndex = RandomSelect(inK,definition);
306 G4int n = NumberOfFinalStates(definition,finalStateIndex);
308 G4double outK = inK - IncomingParticleBindingEnergyConstant(definition,finalStateIndex);
314 if (definition == instance->
GetIon(
"hydrogen")) electronK = inK*electron_mass_c2/proton_mass_c2;
315 else electronK = inK*electron_mass_c2/(particleMass);
319 G4Exception(
"G4DNADingfelderChargeIncreaseModel::SampleSecondaries",
"em0004",
325 (OutgoingParticleDefinition(definition,finalStateIndex),
329 fvect->push_back(dp);
345 G4int finalStateIndex )
350 if (particleDefinition == instance->
GetIon(
"hydrogen"))
return 2;
351 if (particleDefinition == instance->
GetIon(
"alpha+"))
return 2;
353 if (particleDefinition == instance->
GetIon(
"helium"))
354 {
if (finalStateIndex==0)
return 2;
363 G4int finalStateIndex)
367 if (particleDefinition == instance->
GetIon(
"alpha+"))
return instance->
GetIon(
"alpha++");
369 if (particleDefinition == instance->
GetIon(
"helium"))
371 if (finalStateIndex==0)
return instance->
GetIon(
"alpha+");
372 return instance->
GetIon(
"alpha++");
380 G4int finalStateIndex )
383 if(particleDefinition == instance->
GetIon(
"hydrogen"))
return 13.6*eV;
385 if(particleDefinition == instance->
GetIon(
"alpha+"))
392 if(particleDefinition == instance->
GetIon(
"helium"))
397 if (finalStateIndex==0)
return 24.587*eV;
398 return (54.509 + 24.587)*eV;
412 G4int particleTypeIndex = 0;
416 if (particleDefinition == instance->
GetIon(
"alpha+")) particleTypeIndex=0;
417 if (particleDefinition == instance->
GetIon(
"helium")) particleTypeIndex=1;
437 if (x1[index][particleTypeIndex]<x0[index][particleTypeIndex])
447 x1[index][particleTypeIndex]=x0[index][particleTypeIndex] + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) / (c0[index][particleTypeIndex] * d0[index][particleTypeIndex]), 1. / (d0[index][particleTypeIndex] - 1.));
448 b1[index][particleTypeIndex]=(a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex] + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x1[index][particleTypeIndex] - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
454 if (x<x0[index][particleTypeIndex])
455 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
456 else if (x<x1[index][particleTypeIndex])
457 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
459 y=a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
461 return f0[index][particleTypeIndex] * std::pow(10., y)*m*m;
467G4int G4DNADingfelderChargeIncreaseModel::RandomSelect(
G4double k,
471 G4int particleTypeIndex = 0;
475 if (particleDefinition == instance->
GetIon(
"hydrogen"))
return 0;
476 if (particleDefinition == instance->
GetIon(
"alpha+")) particleTypeIndex=0;
477 if (particleDefinition == instance->
GetIon(
"helium")) particleTypeIndex=1;
479 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
487 values[i]=PartialCrossSection(k, i, particleDefinition);
513 G4int particleTypeIndex = 0;
517 if (particleDefinition == instance->
GetIon(
"alpha+")) particleTypeIndex=0;
518 if (particleDefinition == instance->
GetIon(
"helium")) particleTypeIndex=1;
522 for (
G4int i=0; i<numberOfPartialCrossSections[particleTypeIndex]; i++)
524 totalCrossSection += PartialCrossSection(k,i,particleDefinition);
526 return totalCrossSection;
G4DLLIMPORT 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)
static G4DNAMolecularMaterial * Instance()
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
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)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)