68 G4cout <<
"Ion elastic model is constructed " <<
G4endl<<
"Energy range: "
101 G4cout <<
"Calling G4DNAIonElasticModel::Initialise()" <<
G4endl;
108 G4cout <<
"G4DNAIonElasticModel: low energy limit increased from " <<
115 G4cout <<
"G4DNAIonElasticModel: high energy limit decreased from " <<
124 char *path = getenv(
"G4LEDATA");
128 G4Exception(
"G4IonElasticModel::Initialise",
"em0006",
134 std::ostringstream fullFileName;
147 (particleDefinition == protonDef && protonDef != 0)
149 (particleDefinition == hydrogenDef && hydrogenDef != 0)
154 totalXSFile =
"dna/sigma_elastic_proton_HTran";
157 fullFileName << path <<
"/dna/sigmadiff_cumulated_elastic_proton_HTran.dat";
161 (particleDefinition == instance->
GetIon(
"helium") && heliumDef)
163 (particleDefinition == instance->
GetIon(
"alpha+") && alphaplusDef)
165 (particleDefinition == instance->
GetIon(
"alpha++") && alphaplusplusDef)
170 totalXSFile =
"dna/sigma_elastic_alpha_HTran";
173 fullFileName << path <<
"/dna/sigmadiff_cumulated_elastic_alpha_HTran.dat";
178 std::ifstream diffCrossSection(fullFileName.str().c_str());
180 if (!diffCrossSection)
183 description <<
"Missing data file:"
184 <<fullFileName.str().c_str()<<
G4endl;
185 G4Exception(
"G4IonElasticModel::Initialise",
"em0003",
200 while(!diffCrossSection.eof())
204 diffCrossSection>>tDummy>>eDummy;
211 eVecm[tDummy].push_back(0.);
216 if (eDummy !=
eVecm[tDummy].back())
eVecm[tDummy].push_back(eDummy);
225 G4cout <<
"Loaded cross section files for ion elastic model" <<
G4endl;
254 G4cout <<
"Calling CrossSectionPerVolume() of G4DNAIonElasticModel"
264 if(waterDensity!= 0.0)
280 G4Exception(
"G4DNAIonElasticModel::ComputeCrossSectionPerVolume",
"em0002",
287 G4cout <<
"__________________________________" <<
G4endl;
288 G4cout <<
"G4DNAIonElasticModel - XS INFO START" <<
G4endl;
289 G4cout <<
"Kinetic energy(eV)=" << ekin/
eV <<
" particle : " << particleName <<
G4endl;
290 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
291 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
292 G4cout <<
"G4DNAIonElasticModel - XS INFO END" <<
G4endl;
297 return sigma*waterDensity;
304 std::vector<G4DynamicParticle*>* ,
311 G4cout <<
"Calling SampleSecondaries() of G4DNAIonElasticModel" <<
G4endl;
345 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
347 xDir *= std::cos(phi);
348 yDir *= std::sin(phi);
350 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
357 depositEnergyCM = 4. * particleEnergy0 *
fParticle_Mass * water_mass *
362 if (!
statCode && (particleEnergy0 >= depositEnergyCM) )
390 std::vector<G4double>::iterator
t2 = std::upper_bound(
eTdummyVec.begin(),
392 std::vector<G4double>::iterator
t1 = t2 - 1;
394 std::vector<G4double>::iterator e12 = std::upper_bound(
eVecm[(*t1)].begin(),
397 std::vector<G4double>::iterator e11 = e12 - 1;
399 std::vector<G4double>::iterator e22 = std::upper_bound(
eVecm[(*t2)].begin(),
402 std::vector<G4double>::iterator e21 = e22 - 1;
416 if(xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0)
return (0.);
418 theta =
QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12,
419 xs21, xs22, valueT1, valueT2, k, integrDiff);
454 G4double a = (std::log10(xs2) - std::log10(xs1))
455 / (std::log10(e2) - std::log10(e1));
456 G4double b = std::log10(xs2) - a * std::log10(e2);
457 G4double sigma = a * std::log10(e) + b;
502 return Theta(particleDefinition, k /
eV, integrdiff);
514 G4cout <<
"*** WARNING : the G4DNAIonElasticModel class is not "
515 "activated below 100 eV !"
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
virtual G4bool LoadData(const G4String &argFileName)
std::ostringstream G4ExceptionDescription
std::vector< ExP01TrackerHit * > a
static G4ParticleTable * GetParticleTable()
static constexpr double MeV
void SetHighEnergyLimit(G4double)
G4DNAIonElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAIonElasticModel")
const G4ThreeVector & GetMomentumDirection() const
const G4String & GetParticleName() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
const std::vector< G4double > * fpMolWaterDensity
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
G4ParticleChangeForGamma * fParticleChangeForGamma
G4double RandomizeThetaCM(G4double k, G4ParticleDefinition *aParticleDefinition)
G4double LowEnergyLimit() const
G4double Theta(G4ParticleDefinition *aParticleDefinition, G4double k, G4double integrDiff)
G4ParticleDefinition * GetDefinition() const
TriDimensionMap fDiffCrossSectionData
const XML_Char int const XML_Char * value
static G4DNAMolecularMaterial * Instance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
std::vector< G4double > eTdummyVec
static constexpr double eV
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
void SetKillBelowThreshold(G4double threshold)
Hep3Vector cross(const Hep3Vector &) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double LinLinInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4DNACrossSectionDataSet * fpTableData
G4ParticleDefinition * GetIon(const G4String &name)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
void SetLowEnergyLimit(G4double)
static MCTruthManager * instance
static G4DNAGenericIonsManager * Instance(void)
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
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)
G4double GetKineticEnergy() const
virtual G4double FindValue(G4double e, G4int componentId=0) const
static constexpr double cm
G4GLOB_DLL std::ostream G4cout
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Hep3Vector orthogonal() const
G4ParticleChangeForGamma * GetParticleChangeForGamma()
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double HighEnergyLimit() const
void ProposeTrackStatus(G4TrackStatus status)
virtual ~G4DNAIonElasticModel()
virtual void Initialise(const G4ParticleDefinition *particuleDefinition, const G4DataVector &)
static constexpr double pi