59 G4cout <<
"Born ionisation model is constructed " <<
G4endl;
82 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator
pos;
102 G4cout <<
"Calling G4DNABornIonisationModel::Initialise()" <<
G4endl;
106 G4String fileElectron(
"dna/sigma_ionisation_e_born");
107 G4String fileProton(
"dna/sigma_ionisation_p_born");
115 G4double scaleFactor = (1.e-22 / 3.343) *
m*
m;
117 char *path = getenv(
"G4LEDATA");
137 std::ostringstream eFullFileName;
139 if (
fasterCode) eFullFileName << path <<
"/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
140 if (!
fasterCode) eFullFileName << path <<
"/dna/sigmadiff_ionisation_e_born.dat";
142 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
144 if (!eDiffCrossSection)
147 FatalException,
"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat");
150 FatalException,
"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
176 while(!eDiffCrossSection.eof())
180 eDiffCrossSection>>tDummy>>eDummy;
182 for (
int j=0; j<5; j++)
188 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
189 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
193 if (!eDiffCrossSection.eof() && !
fasterCode) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
218 std::ostringstream pFullFileName;
220 if (
fasterCode) pFullFileName << path <<
"/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
222 if (!
fasterCode) pFullFileName << path <<
"/dna/sigmadiff_ionisation_p_born.dat";
224 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
226 if (!pDiffCrossSection)
229 FatalException,
"Missing data file:/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat");
232 FatalException,
"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
236 while(!pDiffCrossSection.eof())
240 pDiffCrossSection>>tDummy>>eDummy;
242 for (
int j=0; j<5; j++)
248 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
249 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
253 if (!pDiffCrossSection.eof() && !
fasterCode) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
261 if (particle==electronDef)
267 if (particle==protonDef)
275 G4cout <<
"Born ionisation model is initialized " << G4endl
305 G4cout <<
"Calling CrossSectionPerVolume() of G4DNABornIonisationModel"
324 if(waterDensity!= 0.0)
329 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
333 lowLim = pos1->second;
336 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
340 highLim = pos2->second;
343 if (ekin >= lowLim && ekin < highLim)
345 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
358 G4Exception(
"G4DNABornIonisationModel::CrossSectionPerVolume",
"em0002",
365 G4cout <<
"__________________________________" <<
G4endl;
366 G4cout <<
"G4DNABornIonisationModel - XS INFO START" <<
G4endl;
367 G4cout <<
"Kinetic energy(eV)=" << ekin/
eV <<
" particle : " << particleName <<
G4endl;
368 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
369 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
370 G4cout <<
"G4DNABornIonisationModel - XS INFO END" <<
G4endl;
374 return sigma*waterDensity;
387 G4cout <<
"Calling SampleSecondaries() of G4DNABornIonisationModel" <<
G4endl;
396 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
401 lowLim = pos1->second;
404 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
409 highLim = pos2->second;
412 if (k >= lowLim && k < highLim)
416 G4double totalEnergy = k + particleMass;
417 G4double pSquare = k * (totalEnergy + particleMass);
418 G4double totalMomentum = std::sqrt(pSquare);
420 G4int ionizationShell = 0;
440 G4int secNumberInit = 0;
441 G4int secNumberFinal = 0;
451 if (ionizationShell <5 && ionizationShell >1)
455 else if (ionizationShell <2)
470 secNumberInit = fvect->size();
472 secNumberFinal = fvect->size();
495 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
497 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
498 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
499 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
500 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
501 finalPx /= finalMomentum;
502 finalPy /= finalMomentum;
503 finalPz /= finalMomentum;
506 direction.set(finalPx,finalPy,finalPz);
514 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
516 for (
G4int j=secNumberInit; j < secNumberFinal; j++)
518 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
525 if (secondaryKinetic>0)
528 fvect->push_back(dp);
549 G4double maximumEnergyTransfer = 0.;
569 G4double maxEnergy = maximumEnergyTransfer;
570 G4int nEnergySteps = 50;
574 std::pow(maxEnergy / value,
575 1. / static_cast<G4double>(nEnergySteps - 1)));
576 G4int step(nEnergySteps);
581 particleDefinition, k /
eV, value /
eV, shell);
582 if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum =
583 differentialCrossSection;
588 G4double secondaryElectronKineticEnergy = 0.;
596 return secondaryElectronKineticEnergy;
602 G4double maximumKineticEnergyTransfer = 4.
603 * (electron_mass_c2 / proton_mass_c2) * k;
610 particleDefinition, k /
eV, value /
eV, shell);
611 if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum =
612 differentialCrossSection;
615 G4double secondaryElectronKineticEnergy = 0.;
618 secondaryElectronKineticEnergy =
G4UniformRand()* maximumKineticEnergyTransfer;
623 return secondaryElectronKineticEnergy;
678 G4int ionizationLevelIndex)
700 std::vector<double>::iterator t2 = std::upper_bound(
eTdummyVec.begin(),
703 std::vector<double>::iterator t1 = t2 - 1;
706 if (energyTransfer <=
eVecm[(*t1)].back() && energyTransfer
707 <=
eVecm[(*t2)].back())
709 std::vector<double>::iterator e12 = std::upper_bound(
710 eVecm[(*t1)].begin(),
eVecm[(*t1)].end(), energyTransfer);
711 std::vector<double>::iterator e11 = e12 - 1;
713 std::vector<double>::iterator e22 = std::upper_bound(
714 eVecm[(*t2)].begin(),
eVecm[(*t2)].end(), energyTransfer);
715 std::vector<double>::iterator e21 = e22 - 1;
736 std::vector<double>::iterator t2 = std::upper_bound(
pTdummyVec.begin(),
738 std::vector<double>::iterator t1 = t2 - 1;
740 std::vector<double>::iterator e12 = std::upper_bound(
pVecm[(*t1)].begin(),
743 std::vector<double>::iterator e11 = e12 - 1;
745 std::vector<double>::iterator e22 = std::upper_bound(
pVecm[(*t2)].begin(),
748 std::vector<double>::iterator e21 = e22 - 1;
764 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
768 xs12, xs21, xs22, valueT1, valueT2, k,
788 if (e1!=0 && e2!=0 && (std::log10(e2)-std::log10(e1)) !=0 && !
fasterCode)
790 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
791 G4double b = std::log10(xs2) - a*std::log10(e2);
792 G4double sigma = a*std::log10(e) + b;
793 value = (std::pow(10.,sigma));
808 if ((e2-e1)!=0 && xs1 !=0 && xs2 !=0 &&
fasterCode )
812 value = std::pow(10.,(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)) );
818 if ((e2-e1)!=0 && (xs1 ==0 || xs2 ==0) &&
fasterCode )
822 value = (d1 + (d2 -
d1)*(e - e1)/ (e2 -
e1));
861 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
879 value += valuesBuffer[i];
890 if (valuesBuffer[i] > value)
892 delete[] valuesBuffer;
895 value -= valuesBuffer[i];
898 if (valuesBuffer)
delete[] valuesBuffer;
904 G4Exception(
"G4DNABornIonisationModel::RandomSelect",
"em0002",
918 G4double secondaryElectronKineticEnergy = 0.;
920 secondaryElectronKineticEnergy=
921 RandomTransferedEnergy(particleDefinition, k/
eV, shell)*
eV-waterStructure.IonisationEnergy(shell);
925 if (secondaryElectronKineticEnergy<0.)
return 0.;
928 return secondaryElectronKineticEnergy;
955 std::vector<double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),
958 std::vector<double>::iterator k1 = k2-1;
972 if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
973 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
975 std::vector<double>::iterator prob12 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
976 eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
978 std::vector<double>::iterator prob11 = prob12-1;
980 std::vector<double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
981 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
983 std::vector<double>::iterator prob21 = prob22-1;
987 valuePROB21 =*prob21;
988 valuePROB22 =*prob22;
989 valuePROB12 =*prob12;
990 valuePROB11 =*prob11;
997 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
998 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
999 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1000 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1012 if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
1014 std::vector<double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1015 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
1017 std::vector<double>::iterator prob21 = prob22-1;
1021 valuePROB21 =*prob21;
1022 valuePROB22 =*prob22;
1026 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1027 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1029 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
1033 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1053 std::vector<double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
1055 std::vector<double>::iterator k1 = k2-1;
1070 if ( random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1071 && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back() )
1073 std::vector<double>::iterator prob12 = std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1074 pProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
1076 std::vector<double>::iterator prob11 = prob12-1;
1078 std::vector<double>::iterator prob22 = std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1079 pProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
1081 std::vector<double>::iterator prob21 = prob22-1;
1085 valuePROB21 =*prob21;
1086 valuePROB22 =*prob22;
1087 valuePROB12 =*prob12;
1088 valuePROB11 =*prob11;
1095 nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1096 nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1097 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1098 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1111 if ( random > pProbaShellMap[ionizationLevelIndex][(*k1)].back() )
1113 std::vector<double>::iterator prob22 = std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1114 pProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
1116 std::vector<double>::iterator prob21 = prob22-1;
1120 valuePROB21 =*prob21;
1121 valuePROB22 =*prob22;
1125 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1126 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1128 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
1132 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1149 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
1152 if (nrjTransfProduct != 0.)
1154 nrj = QuadInterpolator( valuePROB11, valuePROB12,
1155 valuePROB21, valuePROB22,
1156 nrjTransf11, nrjTransf12,
1157 nrjTransf21, nrjTransf22,
static G4Electron * ElectronDefinition()
G4double RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
G4double LowEnergyLimit() const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4LossTableManager * Instance()
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
const std::vector< G4double > * fpMolWaterDensity
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static G4Proton * ProtonDefinition()
G4VEmAngularDistribution * GetAngularDistribution()
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
std::vector< double > eTdummyVec
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
std::vector< double > pTdummyVec
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)
G4VAtomDeexcitation * fAtomDeexcitation
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
virtual ~G4DNABornIonisationModel()
void SetHighEnergyLimit(G4double)
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
TriDimensionMap pDiffCrossSectionData[6]
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
G4DNABornIonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNABornIonisationModel")
virtual G4double FindValue(G4double e, G4int componentId=0) const
G4int RandomSelect(G4double energy, const G4String &particle)
virtual size_t NumberOfComponents(void) const
G4double IonisationEnergy(G4int level)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4DNAWaterIonisationStructure waterStructure
TriDimensionMap eNrjTransfData[6]
static G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
G4double RandomTransferedEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
G4double GetPDGMass() const
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Method used by DNA physics model to create a water molecule.
void SetAngularDistribution(G4VEmAngularDistribution *)
TriDimensionMap pNrjTransfData[6]
G4ParticleChangeForGamma * fParticleChangeForGamma
G4double Interpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
const G4Track * GetCurrentTrack() const
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
TriDimensionMap eDiffCrossSectionData[6]
G4VAtomDeexcitation * AtomDeexcitation()
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
void SetLowEnergyLimit(G4double)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
G4ThreeVector G4ParticleMomentum
G4double bindingEnergy(G4int A, G4int Z)
static const G4double pos
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &=*(new G4DataVector()))
G4ParticleChangeForGamma * GetParticleChangeForGamma()
const G4Material * GetMaterial() const