59 G4cout <<
"Born ionisation model is constructed " <<
G4endl;
110 G4cout <<
"Calling G4DNABornIonisationModel2::Initialise()" <<
G4endl;
116 description <<
"You are trying to initialized G4DNABornIonisationModel2 "
120 description <<
"G4DNABornIonisationModel2 was already initiliased "
122 G4Exception(
"G4DNABornIonisationModel2::Initialise",
"bornIonInit",
129 char *path = getenv(
"G4LEDATA");
134 std::ostringstream fullFileName;
135 fullFileName << path;
137 if(particleName ==
"e-")
145 fullFileName <<
"/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
149 fullFileName <<
"/dna/sigmadiff_ionisation_e_born.dat";
152 else if(particleName ==
"proton")
160 fullFileName <<
"/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
164 fullFileName <<
"/dna/sigmadiff_ionisation_p_born.dat";
169 G4double scaleFactor = (1.e-22 / 3.343) *
m*
m;
175 std::ifstream diffCrossSection(fullFileName.str().c_str());
177 if (!diffCrossSection)
180 description <<
"Missing data file:" <<
G4endl << fullFileName.str() <<
G4endl;
181 G4Exception(
"G4DNABornIonisationModel2::Initialise",
"em0003",
193 for (
int j=0; j<5; j++)
202 while(!diffCrossSection.eof())
206 diffCrossSection>>tDummy>>eDummy;
210 for (
int j=0; j<5; j++)
212 diffCrossSection>> tmp;
219 fProbaShellMap[j][tDummy].push_back(fDiffCrossSectionData[j][tDummy][eDummy]);
236 G4cout <<
"Born ionisation model is initialized " <<
G4endl
267 G4cout <<
"Calling CrossSectionPerVolume() of G4DNABornIonisationModel2"
279 if(waterDensity!= 0.0)
289 G4double A = 1.39241700556072800000E-009 ;
290 G4double B = -8.52610412942622630000E-002 ;
291 sigma = sigma * std::exp(A*(ekin/
eV)+B);
298 G4cout <<
"__________________________________" <<
G4endl;
299 G4cout <<
"G4DNABornIonisationModel2 - XS INFO START" <<
G4endl;
301 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
302 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
303 G4cout <<
"G4DNABornIonisationModel2 - XS INFO END" <<
G4endl;
307 return sigma*waterDensity;
321 G4cout <<
"Calling SampleSecondaries() of G4DNABornIonisationModel2"
331 G4double totalEnergy = k + particleMass;
332 G4double pSquare = k * (totalEnergy + particleMass);
333 G4double totalMomentum = std::sqrt(pSquare);
335 G4int ionizationShell = 0;
355 G4int secNumberInit = 0;
356 G4int secNumberFinal = 0;
362 if (k<bindingEnergy)
return;
370 if (ionizationShell <5 && ionizationShell >1)
374 else if (ionizationShell <2)
389 secNumberInit = fvect->size();
391 secNumberFinal = fvect->size();
414 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
416 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
417 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
418 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
419 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
420 finalPx /= finalMomentum;
421 finalPy /= finalMomentum;
422 finalPz /= finalMomentum;
425 direction.set(finalPx,finalPy,finalPz);
433 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
435 for (
G4int j=secNumberInit; j < secNumberFinal; j++)
437 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
452 if (secondaryKinetic>0)
455 fvect->push_back(dp);
476 G4double maximumEnergyTransfer = 0.;
478 maximumEnergyTransfer = k;
496 G4double maxEnergy = maximumEnergyTransfer;
497 G4int nEnergySteps = 50;
500 G4double stpEnergy(std::pow(maxEnergy / value,
501 1. / static_cast<G4double>(nEnergySteps - 1)));
502 G4int step(nEnergySteps);
511 if (differentialCrossSection >= crossSectionMaximum)
512 crossSectionMaximum = differentialCrossSection;
517 G4double secondaryElectronKineticEnergy = 0.;
525 return secondaryElectronKineticEnergy;
531 G4double maximumKineticEnergyTransfer = 4.
532 * (electron_mass_c2 / proton_mass_c2) * k;
543 if (differentialCrossSection >= crossSectionMaximum)
544 crossSectionMaximum = differentialCrossSection;
547 G4double secondaryElectronKineticEnergy = 0.;
550 secondaryElectronKineticEnergy =
G4UniformRand()* maximumKineticEnergyTransfer;
555 return secondaryElectronKineticEnergy;
609 G4int ionizationLevelIndex)
629 std::vector<double>::iterator t2 = std::upper_bound(
fTdummyVec.begin(),
633 std::vector<double>::iterator t1 = t2 - 1;
636 if (energyTransfer <=
fVecm[(*t1)].back()
637 && energyTransfer <=
fVecm[(*t2)].back())
639 std::vector<double>::iterator e12 = std::upper_bound(
fVecm[(*t1)].begin(),
642 std::vector<double>::iterator e11 = e12 - 1;
644 std::vector<double>::iterator e22 = std::upper_bound(
fVecm[(*t2)].begin(),
647 std::vector<double>::iterator e21 = e22 - 1;
663 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
696 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
699 G4double a = (std::log10(xs2) - std::log10(xs1))
700 / (std::log10(e2) - std::log10(e1));
701 G4double b = std::log10(xs2) - a * std::log10(e2);
702 G4double sigma = a * std::log10(e) + b;
703 value = (std::pow(10., sigma));
717 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 &&
fasterCode)
721 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
727 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) &&
fasterCode)
731 value = (d1 + (d2 -
d1) * (e - e1) / (e2 -
e1));
799 value += valuesBuffer[i];
810 if (valuesBuffer[i] > value)
812 delete[] valuesBuffer;
815 value -= valuesBuffer[i];
819 delete[] valuesBuffer;
832 G4double secondaryElectronKineticEnergy = 0.;
844 if (secondaryElectronKineticEnergy < 0.)
848 return secondaryElectronKineticEnergy;
855 G4int ionizationLevelIndex,
874 std::vector<double>::iterator k2 = std::upper_bound(
fTdummyVec.begin(),
877 std::vector<double>::iterator k1 = k2 - 1;
894 std::vector<double>::iterator prob12 =
895 std::upper_bound(
fProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
899 std::vector<double>::iterator prob11 = prob12 - 1;
901 std::vector<double>::iterator prob22 =
902 std::upper_bound(
fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
906 std::vector<double>::iterator prob21 = prob22 - 1;
910 valuePROB21 = *prob21;
911 valuePROB22 = *prob22;
912 valuePROB12 = *prob12;
913 valuePROB11 = *prob11;
920 nrjTransf11 =
fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
921 nrjTransf12 =
fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
922 nrjTransf21 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
923 nrjTransf22 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
937 std::vector<double>::iterator prob22 =
938 std::upper_bound(
fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
942 std::vector<double>::iterator prob21 = prob22 - 1;
946 valuePROB21 = *prob21;
947 valuePROB22 = *prob22;
951 nrjTransf21 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
952 nrjTransf22 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
979 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
983 if (nrjTransfProduct != 0.)
static G4Electron * ElectronDefinition()
G4double LowEnergyLimit() const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4LossTableManager * Instance()
G4int RandomSelect(G4double energy)
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
G4double Interpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
static G4Proton * ProtonDefinition()
TriDimensionMap fDiffCrossSectionData[6]
G4VEmAngularDistribution * GetAngularDistribution()
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
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)
double B(double temperature)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4DNABornIonisationModel2(const G4ParticleDefinition *p=0, const G4String &nam="DNABornIonisationModel")
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
const G4ThreeVector & GetMomentumDirection() const
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
virtual G4double FindValue(G4double e, G4int componentId=0) const
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual ~G4DNABornIonisationModel2()
const std::vector< G4double > * fpMolWaterDensity
virtual size_t NumberOfComponents(void) const
G4double IonisationEnergy(G4int level)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &=*(new G4DataVector()))
G4double fHighEnergyLimit
G4double TransferedEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, G4double random)
static G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
G4double GetPDGMass() const
virtual G4double GetPartialCrossSection(const G4Material *, G4int, const G4ParticleDefinition *, G4double)
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Method used by DNA physics model to create a water molecule.
G4DNAWaterIonisationStructure waterStructure
void SetAngularDistribution(G4VEmAngularDistribution *)
const G4Track * GetCurrentTrack() const
G4DNACrossSectionDataSet * fTableData
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
const G4ParticleDefinition * fParticleDef
G4VAtomDeexcitation * fAtomDeexcitation
std::vector< double > fTdummyVec
G4VAtomDeexcitation * AtomDeexcitation()
void SetLowEnergyLimit(G4double)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
G4double RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
G4ThreeVector G4ParticleMomentum
G4double bindingEnergy(G4int A, G4int Z)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
const G4Material * GetMaterial() const
TriDimensionMap fNrjTransfData[6]