67 G4cout <<
"Ion elastic model is constructed " <<
G4endl<<
"Energy range: "
100 G4cout <<
"Calling G4DNAIonElasticModel::Initialise()" <<
G4endl;
107 G4cout <<
"G4DNAIonElasticModel: low energy limit increased from " <<
114 G4cout <<
"G4DNAIonElasticModel: high energy limit decreased from " <<
123 char *path = getenv(
"G4LEDATA");
127 G4Exception(
"G4IonElasticModel::Initialise",
"em0006",
133 std::ostringstream fullFileName;
146 (particleDefinition == protonDef && protonDef != 0)
148 (particleDefinition == hydrogenDef && hydrogenDef != 0)
153 totalXSFile =
"dna/sigma_elastic_proton_HTran";
156 fullFileName << path <<
"/dna/sigmadiff_cumulated_elastic_proton_HTran.dat";
160 (particleDefinition == instance->
GetIon(
"helium") && heliumDef)
162 (particleDefinition == instance->
GetIon(
"alpha+") && alphaplusDef)
164 (particleDefinition == instance->
GetIon(
"alpha++") && alphaplusplusDef)
169 totalXSFile =
"dna/sigma_elastic_alpha_HTran";
172 fullFileName << path <<
"/dna/sigmadiff_cumulated_elastic_alpha_HTran.dat";
177 std::ifstream diffCrossSection(fullFileName.str().c_str());
179 if (!diffCrossSection)
182 description <<
"Missing data file:"
183 <<fullFileName.str().c_str()<<
G4endl;
184 G4Exception(
"G4IonElasticModel::Initialise",
"em0003",
199 while(!diffCrossSection.eof())
203 diffCrossSection>>tDummy>>eDummy;
210 eVecm[tDummy].push_back(0.);
215 if (eDummy !=
eVecm[tDummy].back())
eVecm[tDummy].push_back(eDummy);
224 G4cout <<
"Loaded cross section files for ion elastic model" <<
G4endl;
253 G4cout <<
"Calling CrossSectionPerVolume() of G4DNAIonElasticModel"
263 if(waterDensity!= 0.0)
279 G4Exception(
"G4DNAIonElasticModel::ComputeCrossSectionPerVolume",
"em0002",
286 G4cout <<
"__________________________________" <<
G4endl;
287 G4cout <<
"G4DNAIonElasticModel - XS INFO START" <<
G4endl;
288 G4cout <<
"Kinetic energy(eV)=" << ekin/
eV <<
" particle : " << particleName <<
G4endl;
289 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
290 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
291 G4cout <<
"G4DNAIonElasticModel - XS INFO END" <<
G4endl;
296 return sigma*waterDensity;
303 std::vector<G4DynamicParticle*>* ,
310 G4cout <<
"Calling SampleSecondaries() of G4DNAIonElasticModel" <<
G4endl;
344 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
346 xDir *= std::cos(phi);
347 yDir *= std::sin(phi);
349 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
356 depositEnergyCM = 4. * particleEnergy0 *
fParticle_Mass * water_mass *
361 if (!
statCode && (particleEnergy0 >= depositEnergyCM) )
389 std::vector<double>::iterator t2 = std::upper_bound(
eTdummyVec.begin(),
391 std::vector<double>::iterator t1 = t2 - 1;
393 std::vector<double>::iterator e12 = std::upper_bound(
eVecm[(*t1)].begin(),
396 std::vector<double>::iterator e11 = e12 - 1;
398 std::vector<double>::iterator e22 = std::upper_bound(
eVecm[(*t2)].begin(),
401 std::vector<double>::iterator e21 = e22 - 1;
415 if(xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0)
return (0.);
417 theta =
QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12,
418 xs21, xs22, valueT1, valueT2, k, integrDiff);
431 G4double value = (d1 + (d2 -
d1) * (e - e1) / (e2 -
e1));
443 G4double value = std::exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1));
453 G4double a = (std::log10(xs2) - std::log10(xs1))
454 / (std::log10(e2) - std::log10(e1));
455 G4double b = std::log10(xs2) - a * std::log10(e2);
456 G4double sigma = a * std::log10(e) + b;
457 G4double value = (std::pow(10., sigma));
501 return Theta(particleDefinition, k /
eV, integrdiff);
513 G4cout <<
"*** WARNING : the G4DNAIonElasticModel class is not "
514 "activated below 100 eV !"
G4double LowEnergyLimit() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
TriDimensionMap fDiffCrossSectionData
std::vector< double > eTdummyVec
virtual G4bool LoadData(const G4String &argFileName)
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4ParticleDefinition * GetDefinition() const
G4double RandomizeThetaCM(G4double k, G4ParticleDefinition *aParticleDefinition)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void SetKillBelowThreshold(G4double threshold)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
G4double LinLinInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4GLOB_DLL std::ostream G4cout
G4DNACrossSectionDataSet * fpTableData
G4ParticleChangeForGamma * fParticleChangeForGamma
const G4ThreeVector & GetMomentumDirection() const
virtual void Initialise(const G4ParticleDefinition *particuleDefinition, const G4DataVector &)
const std::vector< G4double > * fpMolWaterDensity
static G4DNAGenericIonsManager * Instance(void)
virtual G4double FindValue(G4double e, G4int componentId=0) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
virtual ~G4DNAIonElasticModel()
static G4DNAMolecularMaterial * Instance()
static G4ParticleTable * GetParticleTable()
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
static MCTruthManager * instance
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double Theta(G4ParticleDefinition *aParticleDefinition, G4double k, G4double integrDiff)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4DNAIonElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAIonElasticModel")
G4ParticleDefinition * GetIon(const G4String &name)