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;
366 if (ionizationShell <5 && ionizationShell >1)
370 else if (ionizationShell <2)
385 secNumberInit = fvect->size();
387 secNumberFinal = fvect->size();
410 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
412 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
413 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
414 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
415 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
416 finalPx /= finalMomentum;
417 finalPy /= finalMomentum;
418 finalPz /= finalMomentum;
421 direction.set(finalPx,finalPy,finalPz);
429 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
431 for (
G4int j=secNumberInit; j < secNumberFinal; j++)
433 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
448 if (secondaryKinetic>0)
451 fvect->push_back(dp);
472 G4double maximumEnergyTransfer = 0.;
474 maximumEnergyTransfer = k;
492 G4double maxEnergy = maximumEnergyTransfer;
493 G4int nEnergySteps = 50;
496 G4double stpEnergy(std::pow(maxEnergy / value,
497 1. / static_cast<G4double>(nEnergySteps - 1)));
498 G4int step(nEnergySteps);
507 if (differentialCrossSection >= crossSectionMaximum)
508 crossSectionMaximum = differentialCrossSection;
513 G4double secondaryElectronKineticEnergy = 0.;
521 return secondaryElectronKineticEnergy;
527 G4double maximumKineticEnergyTransfer = 4.
528 * (electron_mass_c2 / proton_mass_c2) * k;
539 if (differentialCrossSection >= crossSectionMaximum)
540 crossSectionMaximum = differentialCrossSection;
543 G4double secondaryElectronKineticEnergy = 0.;
546 secondaryElectronKineticEnergy =
G4UniformRand()* maximumKineticEnergyTransfer;
551 return secondaryElectronKineticEnergy;
605 G4int ionizationLevelIndex)
625 std::vector<double>::iterator t2 = std::upper_bound(
fTdummyVec.begin(),
629 std::vector<double>::iterator t1 = t2 - 1;
632 if (energyTransfer <=
fVecm[(*t1)].back()
633 && energyTransfer <=
fVecm[(*t2)].back())
635 std::vector<double>::iterator e12 = std::upper_bound(
fVecm[(*t1)].begin(),
638 std::vector<double>::iterator e11 = e12 - 1;
640 std::vector<double>::iterator e22 = std::upper_bound(
fVecm[(*t2)].begin(),
643 std::vector<double>::iterator e21 = e22 - 1;
659 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
692 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
695 G4double a = (std::log10(xs2) - std::log10(xs1))
696 / (std::log10(e2) - std::log10(e1));
697 G4double b = std::log10(xs2) - a * std::log10(e2);
698 G4double sigma = a * std::log10(e) + b;
699 value = (std::pow(10., sigma));
713 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 &&
fasterCode)
717 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
723 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) &&
fasterCode)
727 value = (d1 + (d2 -
d1) * (e - e1) / (e2 -
e1));
795 value += valuesBuffer[i];
806 if (valuesBuffer[i] > value)
808 delete[] valuesBuffer;
811 value -= valuesBuffer[i];
815 delete[] valuesBuffer;
828 G4double secondaryElectronKineticEnergy = 0.;
840 if (secondaryElectronKineticEnergy < 0.)
844 return secondaryElectronKineticEnergy;
851 G4int ionizationLevelIndex,
870 std::vector<double>::iterator k2 = std::upper_bound(
fTdummyVec.begin(),
873 std::vector<double>::iterator k1 = k2 - 1;
890 std::vector<double>::iterator prob12 =
891 std::upper_bound(
fProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
895 std::vector<double>::iterator prob11 = prob12 - 1;
897 std::vector<double>::iterator prob22 =
898 std::upper_bound(
fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
902 std::vector<double>::iterator prob21 = prob22 - 1;
906 valuePROB21 = *prob21;
907 valuePROB22 = *prob22;
908 valuePROB12 = *prob12;
909 valuePROB11 = *prob11;
916 nrjTransf11 =
fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
917 nrjTransf12 =
fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
918 nrjTransf21 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
919 nrjTransf22 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
933 std::vector<double>::iterator prob22 =
934 std::upper_bound(
fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
938 std::vector<double>::iterator prob21 = prob22 - 1;
942 valuePROB21 = *prob21;
943 valuePROB22 = *prob22;
947 nrjTransf21 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
948 nrjTransf22 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
975 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
979 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]