44 fpMolWaterDensity = 0;
45 numberOfPartialCrossSections[0] = 0;
46 numberOfPartialCrossSections[1] = 0;
47 numberOfPartialCrossSections[2] = 0;
59 G4cout <<
"Dingfelder charge decrease model is constructed " <<
G4endl;
82 G4cout <<
"Calling G4DNADingfelderChargeDecreaseModel::Initialise()"
101 lowEnergyLimit[proton] = 100. * eV;
102 highEnergyLimit[proton] = 100. * MeV;
105 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
106 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
109 lowEnergyLimit[alphaPlus] = 1. * keV;
110 highEnergyLimit[alphaPlus] = 400. * MeV;
114 if (particle==protonDef)
120 if (particle==alphaPlusPlusDef)
126 if (particle==alphaPlusDef)
145 numberOfPartialCrossSections[0] = 1;
148 f0[0][1]=1.; a0[0][1]=0.95;
169 numberOfPartialCrossSections[1] = 2;
183 numberOfPartialCrossSections[2] = 1;
189 G4cout <<
"Dingfelder charge decrease model is initialized " <<
G4endl
203 isInitialised =
true;
215 if (verboseLevel > 3)
218 <<
"Calling CrossSectionPerVolume() of G4DNADingfelderChargeDecreaseModel"
230 particleDefinition != instance->
GetIon(
"alpha++")
232 particleDefinition != instance->
GetIon(
"alpha+")
245 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
246 pos1 = lowEnergyLimit.find(particleName);
248 if (pos1 != lowEnergyLimit.end())
250 lowLim = pos1->second;
253 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
254 pos2 = highEnergyLimit.find(particleName);
256 if (pos2 != highEnergyLimit.end())
258 highLim = pos2->second;
261 if (k >= lowLim && k <= highLim)
263 crossSection = Sum(k,particleDefinition);
266 if (verboseLevel > 2)
268 G4cout <<
"_______________________________________" <<
G4endl;
269 G4cout <<
"G4DNADingfelderChargeDecreaeModel" <<
G4endl;
270 G4cout <<
"Kinetic energy(eV)=" << k/eV <<
"particle :" << particleName <<
G4endl;
271 G4cout <<
"Cross section per water molecule (cm^2)=" << crossSection/cm/cm <<
G4endl;
272 G4cout <<
"Cross section per water molecule (cm^-1)=" << crossSection*
273 waterDensity/(1./cm) <<
G4endl;
277 return crossSection*waterDensity;
291 if (verboseLevel > 3)
294 <<
"Calling SampleSecondaries() of G4DNADingfelderChargeDecreaseModel"
304 G4int finalStateIndex = RandomSelect(inK,definition);
306 G4int n = NumberOfFinalStates(definition, finalStateIndex);
307 G4double waterBindingEnergy = WaterBindingEnergyConstant(definition, finalStateIndex);
308 G4double outgoingParticleBindingEnergy = OutgoingParticleBindingEnergyConstant(definition, finalStateIndex);
314 if (!statCode) outK = inK - n*(inK*electron_mass_c2/proton_mass_c2)
315 - waterBindingEnergy + outgoingParticleBindingEnergy;
321 if (!statCode) outK = inK - n*(inK*electron_mass_c2/particleMass)
322 - waterBindingEnergy + outgoingParticleBindingEnergy;
328 G4Exception(
"G4DNADingfelderChargeDecreaseModel::SampleSecondaries",
"em0004",
341 + waterBindingEnergy - outgoingParticleBindingEnergy);
344 + waterBindingEnergy - outgoingParticleBindingEnergy);
350 fvect->push_back(dp);
362 G4int finalStateIndex)
371 if (particleDefinition == instance->
GetIon(
"alpha++"))
373 if (finalStateIndex == 0)
378 if (particleDefinition == instance->
GetIon(
"alpha+"))
387 G4int finalStateIndex)
392 return instance->
GetIon(
"hydrogen");
394 if (particleDefinition == instance->
GetIon(
"alpha++"))
396 if (finalStateIndex == 0)
397 return instance->
GetIon(
"alpha+");
398 return instance->
GetIon(
"helium");
401 if (particleDefinition == instance->
GetIon(
"alpha+"))
402 return instance->
GetIon(
"helium");
410 G4int finalStateIndex)
421 if (particleDefinition == instance->
GetIon(
"alpha++"))
429 if (finalStateIndex == 0)
432 return 10.79 * 2 * eV;
435 if (particleDefinition == instance->
GetIon(
"alpha+"))
452 G4int finalStateIndex)
459 if (particleDefinition == instance->
GetIon(
"alpha++"))
464 if (finalStateIndex == 0)
467 return (54.509 + 24.587) * eV;
470 if (particleDefinition == instance->
GetIon(
"alpha+"))
487 G4int particleTypeIndex = 0;
492 particleTypeIndex = 0;
494 if (particleDefinition == instance->
GetIon(
"alpha++"))
495 particleTypeIndex = 1;
497 if (particleDefinition == instance->
GetIon(
"alpha+"))
498 particleTypeIndex = 2;
518 if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex])
528 x1[index][particleTypeIndex] = x0[index][particleTypeIndex]
529 + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
530 / (c0[index][particleTypeIndex]
531 * d0[index][particleTypeIndex]),
532 1. / (d0[index][particleTypeIndex] - 1.));
533 b1[index][particleTypeIndex] = (a0[index][particleTypeIndex]
534 - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
535 + b0[index][particleTypeIndex]
536 - c0[index][particleTypeIndex]
537 * std::pow(x1[index][particleTypeIndex]
538 - x0[index][particleTypeIndex],
539 d0[index][particleTypeIndex]);
545 if (x < x0[index][particleTypeIndex])
546 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
547 else if (x < x1[index][particleTypeIndex])
548 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]
549 - c0[index][particleTypeIndex]
550 * std::pow(x - x0[index][particleTypeIndex],
551 d0[index][particleTypeIndex]);
553 y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
555 return f0[index][particleTypeIndex] * std::pow(10., y) * m * m;
559G4int G4DNADingfelderChargeDecreaseModel::RandomSelect(
G4double k,
562 G4int particleTypeIndex = 0;
567 particleTypeIndex = 0;
569 if (particleDefinition == instance->
GetIon(
"alpha++"))
570 particleTypeIndex = 1;
572 if (particleDefinition == instance->
GetIon(
"alpha+"))
573 particleTypeIndex = 2;
575 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
583 values[i] = PartialCrossSection(k, i, particleDefinition);
594 if (values[i] > value)
610 G4int particleTypeIndex = 0;
615 particleTypeIndex = 0;
617 if (particleDefinition == instance->
GetIon(
"alpha++"))
618 particleTypeIndex = 1;
620 if (particleDefinition == instance->
GetIon(
"alpha+"))
621 particleTypeIndex = 2;
625 for (
G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++)
627 totalCrossSection += PartialCrossSection(k, i, particleDefinition);
630 return totalCrossSection;
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4GLOB_DLL std::ostream G4cout
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4DNADingfelderChargeDecreaseModel(const G4ParticleDefinition *p=0, const G4String &nam="DNADingfelderChargeDecreaseModel")
virtual ~G4DNADingfelderChargeDecreaseModel()
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 G4Material * GetMaterial(const G4String &name, G4bool warning=true)
const G4Track * GetCurrentTrack() const
G4double GetPDGMass() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
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)