45 fpMolWaterDensity = 0;
46 numberOfPartialCrossSections[0] = 0;
47 numberOfPartialCrossSections[1] = 0;
48 numberOfPartialCrossSections[2] = 0;
60 G4cout <<
"Dingfelder charge decrease model is constructed " <<
G4endl;
77 G4cout <<
"Calling G4DNADingfelderChargeDecreaseModel::Initialise()" <<
G4endl;
94 lowEnergyLimit[proton] = 100. * eV;
95 highEnergyLimit[proton] = 100. * MeV;
98 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
99 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
102 lowEnergyLimit[alphaPlus] = 1. * keV;
103 highEnergyLimit[alphaPlus] = 400. * MeV;
107 if (particle==protonDef)
113 if (particle==alphaPlusPlusDef)
119 if (particle==alphaPlusDef)
138 numberOfPartialCrossSections[0] = 1;
141 f0[0][1]=1.; a0[0][1]=0.95;
162 numberOfPartialCrossSections[1] = 2;
176 numberOfPartialCrossSections[2] = 1;
182 G4cout <<
"Dingfelder charge decrease model is initialized " <<
G4endl
193 if (isInitialised) {
return; }
195 isInitialised =
true;
207 if (verboseLevel > 3)
208 G4cout <<
"Calling CrossSectionPerVolume() of G4DNADingfelderChargeDecreaseModel" <<
G4endl;
218 particleDefinition != instance->
GetIon(
"alpha++")
220 particleDefinition != instance->
GetIon(
"alpha+")
231 if(waterDensity!= 0.0)
236 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
237 pos1 = lowEnergyLimit.find(particleName);
239 if (pos1 != lowEnergyLimit.end())
241 lowLim = pos1->second;
244 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
245 pos2 = highEnergyLimit.find(particleName);
247 if (pos2 != highEnergyLimit.end())
249 highLim = pos2->second;
252 if (k >= lowLim && k < highLim)
254 crossSection = Sum(k,particleDefinition);
257 if (verboseLevel > 2)
259 G4cout <<
"_______________________________________" <<
G4endl;
260 G4cout <<
"°°° G4DNADingfelderChargeDecreaeModel" <<
G4endl;
261 G4cout <<
"---> Kinetic energy(eV)=" << k/eV <<
"particle :" << particleName <<
G4endl;
262 G4cout <<
" - Cross section per water molecule (cm^2)=" << crossSection/cm/cm <<
G4endl;
263 G4cout <<
" - Cross section per water molecule (cm^-1)=" << crossSection*
264 waterDensity/(1./cm) <<
G4endl;
269 return crossSection*waterDensity;
282 if (verboseLevel > 3)
283 G4cout <<
"Calling SampleSecondaries() of G4DNADingfelderChargeDecreaseModel" <<
G4endl;
291 G4int finalStateIndex = RandomSelect(inK,definition);
293 G4int n = NumberOfFinalStates(definition, finalStateIndex);
294 G4double waterBindingEnergy = WaterBindingEnergyConstant(definition, finalStateIndex);
295 G4double outgoingParticleBindingEnergy = OutgoingParticleBindingEnergyConstant(definition, finalStateIndex);
299 outK = inK - n*(inK*electron_mass_c2/proton_mass_c2) - waterBindingEnergy + outgoingParticleBindingEnergy;
301 outK = inK - n*(inK*electron_mass_c2/particleMass) - waterBindingEnergy + outgoingParticleBindingEnergy;
305 G4Exception(
"G4DNADingfelderChargeDecreaseModel::SampleSecondaries",
"em0004",
315 fvect->push_back(dp);
322 G4int finalStateIndex )
330 if (particleDefinition == instance->
GetIon(
"alpha++") )
332 if (finalStateIndex==0)
return 1;
336 if (particleDefinition == instance->
GetIon(
"alpha+") )
return 1;
344 G4int finalStateIndex)
350 if (particleDefinition == instance->
GetIon(
"alpha++") )
352 if (finalStateIndex == 0)
return instance->
GetIon(
"alpha+");
353 return instance->
GetIon(
"helium");
356 if (particleDefinition == instance->
GetIon(
"alpha+") )
return instance->
GetIon(
"helium");
364 G4int finalStateIndex )
374 if (particleDefinition == instance->
GetIon(
"alpha++") )
382 if (finalStateIndex == 0)
return 10.79*eV;
387 if (particleDefinition == instance->
GetIon(
"alpha+") )
404 G4int finalStateIndex)
410 if (particleDefinition == instance->
GetIon(
"alpha++") )
415 if (finalStateIndex==0)
return 54.509*eV;
417 return (54.509 + 24.587)*eV;
420 if (particleDefinition == instance->
GetIon(
"alpha+") )
436 G4int particleTypeIndex = 0;
441 if (particleDefinition == instance->
GetIon(
"alpha++")) particleTypeIndex=1;
442 if (particleDefinition == instance->
GetIon(
"alpha+")) particleTypeIndex=2;
462 if (x1[index][particleTypeIndex]<x0[index][particleTypeIndex])
472 x1[index][particleTypeIndex]=x0[index][particleTypeIndex] + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
473 / (c0[index][particleTypeIndex] * d0[index][particleTypeIndex]), 1. / (d0[index][particleTypeIndex] - 1.));
474 b1[index][particleTypeIndex]=(a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
475 + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x1[index][particleTypeIndex]
476 - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
482 if (x<x0[index][particleTypeIndex])
483 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
484 else if (x<x1[index][particleTypeIndex])
485 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] - c0[index][particleTypeIndex]
486 * std::pow(x - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
488 y=a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
490 return f0[index][particleTypeIndex] * std::pow(10., y)*m*m;
494G4int G4DNADingfelderChargeDecreaseModel::RandomSelect(
G4double k,
497 G4int particleTypeIndex = 0;
502 if (particleDefinition == instance->
GetIon(
"alpha++")) particleTypeIndex = 1;
503 if (particleDefinition == instance->
GetIon(
"alpha+")) particleTypeIndex = 2;
505 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
513 values[i]=PartialCrossSection(k, i, particleDefinition);
539 G4int particleTypeIndex = 0;
544 if (particleDefinition == instance->
GetIon(
"alpha++")) particleTypeIndex=1;
545 if (particleDefinition == instance->
GetIon(
"alpha+")) particleTypeIndex=2;
549 for (
G4int i=0; i<numberOfPartialCrossSections[particleTypeIndex]; i++)
551 totalCrossSection += PartialCrossSection(k,i,particleDefinition);
553 return totalCrossSection;
G4DLLIMPORT std::ostream G4cout
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)
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)
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)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)