57 G4cout <<
"Born ionisation model is constructed " <<
G4endl;
62 fAtomDeexcitation = 0;
64 fpMolWaterDensity = 0;
77 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
78 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
98 G4cout <<
"Calling G4DNABornIonisationModel::Initialise()" <<
G4endl;
102 G4String fileElectron(
"dna/sigma_ionisation_e_born");
103 G4String fileProton(
"dna/sigma_ionisation_p_born");
111 G4double scaleFactor = (1.e-22 / 3.343) *
m*
m;
113 char *path = getenv(
"G4LEDATA");
133 std::ostringstream eFullFileName;
135 if (fasterCode) eFullFileName << path <<
"/dna/sigmadiff_cumulated_ionisation_e_born.dat";
136 if (!fasterCode) eFullFileName << path <<
"/dna/sigmadiff_ionisation_e_born.dat";
138 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
140 if (!eDiffCrossSection)
142 if (fasterCode)
G4Exception(
"G4DNABornIonisationModel::Initialise",
"em0003",
143 FatalException,
"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_born.dat");
145 if (!fasterCode)
G4Exception(
"G4DNABornIonisationModel::Initialise",
"em0003",
146 FatalException,
"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
149 eTdummyVec.push_back(0.);
150 while(!eDiffCrossSection.eof())
154 eDiffCrossSection>>tDummy>>eDummy;
155 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
156 for (
int j=0; j<5; j++)
158 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
162 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
163 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
167 if (!eDiffCrossSection.eof() && !fasterCode) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
169 if (!fasterCode) eVecm[tDummy].push_back(eDummy);
178 tableFile[
proton] = fileProton;
188 tableData[
proton] = tableP;
192 std::ostringstream pFullFileName;
194 if (fasterCode) pFullFileName << path <<
"/dna/sigmadiff_cumulated_ionisation_p_born.dat";
196 if (!fasterCode) pFullFileName << path <<
"/dna/sigmadiff_ionisation_p_born.dat";
198 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
200 if (!pDiffCrossSection)
202 if (fasterCode)
G4Exception(
"G4DNABornIonisationModel::Initialise",
"em0003",
203 FatalException,
"Missing data file:/dna/sigmadiff_cumulated_ionisation_p_born.dat");
205 if (!fasterCode)
G4Exception(
"G4DNABornIonisationModel::Initialise",
"em0003",
206 FatalException,
"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
209 pTdummyVec.push_back(0.);
210 while(!pDiffCrossSection.eof())
214 pDiffCrossSection>>tDummy>>eDummy;
215 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
216 for (
int j=0; j<5; j++)
218 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
222 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
223 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
227 if (!pDiffCrossSection.eof() && !fasterCode) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
229 if (!fasterCode) pVecm[tDummy].push_back(eDummy);
235 if (particle==electronDef)
241 if (particle==protonDef)
249 G4cout <<
"Born ionisation model is initialized " << G4endl
263 if (isInitialised) {
return; }
265 isInitialised =
true;
276 if (verboseLevel > 3)
277 G4cout <<
"Calling CrossSectionPerVolume() of G4DNABornIonisationModel" <<
G4endl;
295 if(waterDensity!= 0.0)
300 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
301 pos1 = lowEnergyLimit.find(particleName);
302 if (pos1 != lowEnergyLimit.end())
304 lowLim = pos1->second;
307 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
308 pos2 = highEnergyLimit.find(particleName);
309 if (pos2 != highEnergyLimit.end())
311 highLim = pos2->second;
314 if (ekin >= lowLim && ekin < highLim)
316 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
317 pos = tableData.find(particleName);
319 if (pos != tableData.end())
329 G4Exception(
"G4DNABornIonisationModel::CrossSectionPerVolume",
"em0002",
334 if (verboseLevel > 2)
336 G4cout <<
"__________________________________" <<
G4endl;
337 G4cout <<
"G4DNABornIonisationModel - XS INFO START" <<
G4endl;
338 G4cout <<
"Kinetic energy(eV)=" << ekin/
eV <<
" particle : " << particleName <<
G4endl;
339 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
340 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
341 G4cout <<
"G4DNABornIonisationModel - XS INFO END" <<
G4endl;
346 return sigma*waterDensity;
358 if (verboseLevel > 3)
359 G4cout <<
"Calling SampleSecondaries() of G4DNABornIonisationModel" <<
G4endl;
368 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
369 pos1 = lowEnergyLimit.find(particleName);
371 if (pos1 != lowEnergyLimit.end())
373 lowLim = pos1->second;
376 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
377 pos2 = highEnergyLimit.find(particleName);
379 if (pos2 != highEnergyLimit.end())
381 highLim = pos2->second;
384 if (k >= lowLim && k < highLim)
388 G4double totalEnergy = k + particleMass;
389 G4double pSquare = k * (totalEnergy + particleMass);
390 G4double totalMomentum = std::sqrt(pSquare);
392 G4int ionizationShell = -1;
394 if (!fasterCode) ionizationShell = RandomSelect(k,particleName);
405 ionizationShell = RandomSelect(k,particleName);
413 G4int secNumberInit = 0;
414 G4int secNumberFinal = 0;
419 if(fAtomDeexcitation) {
423 if (ionizationShell <5 && ionizationShell >1)
427 else if (ionizationShell <2)
442 secNumberInit = fvect->size();
444 secNumberFinal = fvect->size();
449 if (!fasterCode) secondaryKinetic = RandomizeEjectedElectronEnergy(particle->
GetDefinition(),k,ionizationShell);
454 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->
GetDefinition(),k,ionizationShell);
457 while (secondaryKinetic<0) ;
461 RandomizeEjectedElectronDirection(particle->
GetDefinition(), k,secondaryKinetic, cosTheta, phi);
463 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
464 G4double dirX = sinTheta*std::cos(phi);
465 G4double dirY = sinTheta*std::sin(phi);
468 deltaDirection.
rotateUz(primaryDirection);
474 G4double finalPx = totalMomentum*primaryDirection.
x() - deltaTotalMomentum*deltaDirection.
x();
475 G4double finalPy = totalMomentum*primaryDirection.
y() - deltaTotalMomentum*deltaDirection.
y();
476 G4double finalPz = totalMomentum*primaryDirection.
z() - deltaTotalMomentum*deltaDirection.
z();
477 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
478 finalPx /= finalMomentum;
479 finalPy /= finalMomentum;
480 finalPz /= finalMomentum;
483 direction.
set(finalPx,finalPy,finalPz);
491 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
493 for (
G4int j=secNumberInit; j < secNumberFinal; j++) {
495 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
504 fvect->push_back(dp);
524 if ((k+waterStructure.
IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
542 G4double maxEnergy = maximumEnergyTransfer;
543 G4int nEnergySteps = 50;
546 G4double stpEnergy(std::pow(maxEnergy/
value, 1./static_cast<G4double>(nEnergySteps-1)));
547 G4int step(nEnergySteps);
552 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
557 G4double secondaryElectronKineticEnergy=0.;
565 return secondaryElectronKineticEnergy;
579 if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
582 G4double secondaryElectronKineticEnergy = 0.;
585 secondaryElectronKineticEnergy =
G4UniformRand() * maximumKineticEnergyTransfer;
590 return secondaryElectronKineticEnergy;
598 void G4DNABornIonisationModel::RandomizeEjectedElectronDirection(
G4ParticleDefinition* particleDefinition,
608 else if (secKinetic <= 200.*
eV)
616 cosTheta = std::sqrt(1.-sin2O);
629 if (secKinetic>100*
eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
640 G4int ionizationLevelIndex)
644 if (energyTransfer >= waterStructure.
IonisationEnergy(ionizationLevelIndex))
662 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
664 std::vector<double>::iterator
t1 = t2-1;
667 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
669 std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
670 std::vector<double>::iterator e11 = e12-1;
672 std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
673 std::vector<double>::iterator e21 = e22-1;
682 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
683 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
684 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
685 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
694 std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
695 std::vector<double>::iterator
t1 = t2-1;
697 std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
698 std::vector<double>::iterator e11 = e12-1;
700 std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
701 std::vector<double>::iterator e21 = e22-1;
710 xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
711 xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
712 xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
713 xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
717 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
720 sigma = QuadInterpolator( valueE11, valueE12,
746 if (e1!=0 && e2!=0 && (std::log10(e2)-std::log10(e1)) !=0 && !fasterCode)
748 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
749 G4double b = std::log10(xs2) - a*std::log10(e2);
751 value = (std::pow(10.,sigma));
766 if ((e2-e1)!=0 && xs1 !=0 && xs2 !=0 && fasterCode )
770 value = std::pow(10.,(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)) );
776 if ((e2-e1)!=0 && (xs1 ==0 || xs2 ==0) && fasterCode )
780 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
806 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
807 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
808 G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
819 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
820 pos = tableData.find(particle);
822 if (pos != tableData.end())
837 value += valuesBuffer[i];
848 if (valuesBuffer[i] > value)
850 delete[] valuesBuffer;
853 value -= valuesBuffer[i];
856 if (valuesBuffer)
delete[] valuesBuffer;
862 G4Exception(
"G4DNABornIonisationModel::RandomSelect",
"em0002",
871 G4double G4DNABornIonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs
876 G4double secondaryElectronKineticEnergy = 0.;
878 secondaryElectronKineticEnergy=
879 RandomTransferedEnergy(particleDefinition, k/
eV, shell)*
eV-waterStructure.
IonisationEnergy(shell);
881 return secondaryElectronKineticEnergy;
886 G4double G4DNABornIonisationModel::RandomTransferedEnergy
910 std::vector<double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
912 std::vector<double>::iterator k1 = k2-1;
927 if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
928 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
931 std::vector<double>::iterator prob12 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
932 eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
934 std::vector<double>::iterator prob11 = prob12-1;
937 std::vector<double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
938 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
940 std::vector<double>::iterator prob21 = prob22-1;
944 valuePROB21 =*prob21;
945 valuePROB22 =*prob22;
946 valuePROB12 =*prob12;
947 valuePROB11 =*prob11;
955 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
956 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
957 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
958 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
973 if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
976 std::vector<double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
977 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
979 std::vector<double>::iterator prob21 = prob22-1;
983 valuePROB21 =*prob21;
984 valuePROB22 =*prob22;
988 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
989 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
991 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
995 G4double value = Interpolate(0., valueK2, k, 0., interpolatedvalue2);
1019 std::vector<double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
1021 std::vector<double>::iterator k1 = k2-1;
1037 if ( random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1038 && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back() )
1040 std::vector<double>::iterator prob12 = std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1041 pProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
1043 std::vector<double>::iterator prob11 = prob12-1;
1046 std::vector<double>::iterator prob22 = std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1047 pProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
1049 std::vector<double>::iterator prob21 = prob22-1;
1053 valuePROB21 =*prob21;
1054 valuePROB22 =*prob22;
1055 valuePROB12 =*prob12;
1056 valuePROB11 =*prob11;
1063 nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1064 nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1065 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1066 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1079 if ( random > pProbaShellMap[ionizationLevelIndex][(*k1)].back() )
1082 std::vector<double>::iterator prob22 = std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1083 pProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
1085 std::vector<double>::iterator prob21 = prob22-1;
1089 valuePROB21 =*prob21;
1090 valuePROB22 =*prob22;
1094 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1095 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1097 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
1101 G4double value = Interpolate(0., valueK2, k, 0., interpolatedvalue2);
1120 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
1124 if (nrjTransfProduct != 0.)
1126 nrj = QuadInterpolator( valuePROB11, valuePROB12,
1127 valuePROB21, valuePROB22,
1128 nrjTransf11, nrjTransf12,
1129 nrjTransf21, nrjTransf22,
void set(double x, double y, double z)
static G4Electron * ElectronDefinition()
G4double LowEnergyLimit() const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4LossTableManager * Instance()
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
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()
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
virtual ~G4DNABornIonisationModel()
void SetHighEnergyLimit(G4double)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4GLOB_DLL std::ostream G4cout
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
G4DNABornIonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNABornIonisationModel")
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual size_t NumberOfComponents(void) const
G4double IonisationEnergy(G4int level)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
G4double GetPDGMass() const
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
G4ParticleChangeForGamma * fParticleChangeForGamma
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
const G4Track * GetCurrentTrack() const
const XML_Char int const XML_Char * value
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4VAtomDeexcitation * AtomDeexcitation()
void SetLowEnergyLimit(G4double)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
G4double bindingEnergy(G4int A, G4int Z)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &=*(new G4DataVector()))
G4ParticleChangeForGamma * GetParticleChangeForGamma()