69 G4cout <<
"CPA100 ionisation model is constructed " <<
G4endl;
102 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
122 G4cout <<
"Calling G4DNACPA100IonisationModel::Initialise()" <<
G4endl;
128 G4String fileElectron(
"dna/sigma_ionisation_e_cpa100_form_rel");
136 char *path = getenv(
"G4LEDATA");
170 std::ostringstream eFullFileName;
172 if (
fasterCode) eFullFileName << path <<
"/dna/sigmadiff_cumulated_ionisation_e_cpa100_rel.dat";
174 if (!
fasterCode) eFullFileName << path <<
"/dna/sigmadiff_ionisation_e_cpa100_rel.dat";
176 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
178 if (!eDiffCrossSection)
181 FatalException,
"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_cpa100_rel.dat");
184 FatalException,
"Missing data file:/dna/sigmadiff_ionisation_e_cpa100_rel.dat");
199 while(!eDiffCrossSection.eof())
203 eDiffCrossSection>>tDummy>>eDummy;
205 for (
G4int j=0; j<5; j++)
211 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
212 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
216 if (!eDiffCrossSection.eof() && !
fasterCode) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
231 if (particle==electronDef)
239 G4cout <<
"CPA100 ionisation model is initialized " << G4endl
268 G4cout <<
"Calling CrossSectionPerVolume() of G4DNACPA100IonisationModel" <<
G4endl;
280 if(waterDensity!= 0.0)
285 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
289 lowLim = pos1->second;
292 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
296 highLim = pos2->second;
299 if (ekin > lowLim && ekin < highLim)
301 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
314 G4Exception(
"G4DNACPA100IonisationModel::CrossSectionPerVolume",
"em0002",
321 G4cout <<
"__________________________________" <<
G4endl;
322 G4cout <<
"G4DNACPA100IonisationModel - XS INFO START" <<
G4endl;
323 G4cout <<
"Kinetic energy(eV)=" << ekin/
eV <<
" particle : " << particleName <<
G4endl;
324 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
325 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
326 G4cout <<
"G4DNACPA100IonisationModel - XS INFO END" <<
G4endl;
331 return sigma*waterDensity;
343 G4cout <<
"Calling SampleSecondaries() of G4DNACPA100IonisationModel" <<
G4endl;
352 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
357 lowLim = pos1->second;
360 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
365 highLim = pos2->second;
368 if (k >= lowLim && k < highLim)
372 G4double totalEnergy = k + particleMass;
373 G4double pSquare = k * (totalEnergy + particleMass);
374 G4double totalMomentum = std::sqrt(pSquare);
376 G4int ionizationShell = -1;
409 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
410 G4double dirX = sinTheta*std::cos(phi);
411 G4double dirY = sinTheta*std::sin(phi);
414 deltaDirection.
rotateUz(primaryDirection);
417 if (secondaryKinetic>0)
420 fvect->push_back(dp);
428 G4double finalPx = totalMomentum*primaryDirection.
x() - deltaTotalMomentum*deltaDirection.
x();
429 G4double finalPy = totalMomentum*primaryDirection.
y() - deltaTotalMomentum*deltaDirection.
y();
430 G4double finalPz = totalMomentum*primaryDirection.
z() - deltaTotalMomentum*deltaDirection.
z();
431 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
432 finalPx /= finalMomentum;
433 finalPy /= finalMomentum;
434 finalPz /= finalMomentum;
437 direction.
set(finalPx,finalPy,finalPz);
450 G4int secNumberInit = 0;
451 G4int secNumberFinal = 0;
458 if (ionizationShell <5 && ionizationShell >1)
462 else if (ionizationShell <2)
477 secNumberInit = fvect->size();
479 secNumberFinal = fvect->size();
483 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
485 for (
G4int j=secNumberInit; j < secNumberFinal; j++) {
486 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
543 G4double maxEnergy = maximumEnergyTransfer;
546 G4int nEnergySteps = 50;
567 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
568 G4int step(nEnergySteps);
569 G4double differentialCrossSection = 0.;
574 if(differentialCrossSection >0)
576 crossSectionMaximum=differentialCrossSection;
584 G4double secondaryElectronKineticEnergy=0.;
592 return secondaryElectronKineticEnergy;
610 cosTheta = std::sqrt(1.-sin2O);
619 G4int ionizationLevelIndex)
643 std::vector<G4double>::iterator
t1 = t2-1;
647 if (energyTransfer <=
eVecm[(*t1)].back() && energyTransfer <=
eVecm[(*t2)].back() )
649 std::vector<G4double>::iterator e12 = std::upper_bound(
eVecm[(*t1)].begin(),
eVecm[(*t1)].end(), energyTransfer);
650 std::vector<G4double>::iterator e11 = e12-1;
652 std::vector<G4double>::iterator e22 = std::upper_bound(
eVecm[(*t2)].begin(),
eVecm[(*t2)].end(), energyTransfer);
653 std::vector<G4double>::iterator e21 = e22-1;
671 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
699 if (e1!=0 && e2!=0 && (std::log10(e2)-std::log10(e1)) !=0 && !
fasterCode &&
useDcs)
701 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
702 G4double b = std::log10(xs2) - a*std::log10(e2);
703 G4double sigma = a*std::log10(e) + b;
704 value = (std::pow(10.,sigma));
723 value = std::pow(10.,(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)) );
733 value = (d1 + (d2 -
d1)*(e - e1)/ (e2 - e1));
772 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
854 value += valuesBuffer[i];
864 if (valuesBuffer[i] > value)
866 delete[] valuesBuffer;
870 value -= valuesBuffer[i];
873 if (valuesBuffer)
delete[] valuesBuffer;
879 G4Exception(
"G4DNACPA100IonisationModel::RandomSelect",
"em0002",
893 G4double secondaryElectronKineticEnergy = 0.;
895 secondaryElectronKineticEnergy=
896 RandomTransferedEnergy(particleDefinition, k/
eV, shell)*
eV-waterStructure.IonisationEnergy(shell);
899 if (secondaryElectronKineticEnergy<0.)
return 0.;
902 return secondaryElectronKineticEnergy;
932 std::vector<G4double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
934 std::vector<G4double>::iterator k1 = k2-1;
949 if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
950 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
954 std::vector<G4double>::iterator prob12 =
955 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
956 eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
958 std::vector<G4double>::iterator prob11 = prob12-1;
961 std::vector<G4double>::iterator prob22 =
962 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
963 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
965 std::vector<G4double>::iterator prob21 = prob22-1;
969 valuePROB21 =*prob21;
970 valuePROB22 =*prob22;
971 valuePROB12 =*prob12;
972 valuePROB11 =*prob11;
980 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
981 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
982 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
983 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
998 if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
1002 std::vector<G4double>::iterator prob22 =
1004 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1005 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
1007 std::vector<G4double>::iterator prob21 = prob22-1;
1011 valuePROB21 =*prob21;
1012 valuePROB22 =*prob22;
1016 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1017 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1019 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
1023 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1044 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
1048 if (nrjTransfProduct != 0.)
1050 nrj = QuadInterpolator( valuePROB11, valuePROB12,
1051 valuePROB21, valuePROB22,
1052 nrjTransf11, nrjTransf12,
1053 nrjTransf21, nrjTransf22,
1215 G4double bb = waterStructure.IonisationEnergy(shell);
1216 G4double uu = waterStructure.UEnergy(shell);
1218 if (tt<=bb)
return 0.;
1228 G4double a1 = t * tm1 / tu1 / tp12;
1229 G4double a2 = tm1 / tu1 / t / tp1 / deux;
1230 G4double a3 = dlt * (tp12 - deux * deux ) / tu1 / tp12;
1255 else if ((r1>A1) && (r1< A2))
1274 gx=deux*(un-(t-wx)/tp1);
1280 fx=un-r2*(tp12-deux*deux)/tp12;
1283 gx=(un+gg*gg*gg)/deux;
void set(double x, double y, double z)
TriDimensionMap eDiffCrossSectionData[6]
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &=*(new G4DataVector()))
virtual G4bool LoadData(const G4String &argFileName)
std::vector< ExP01TrackerHit * > a
static const G4double pos
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4VAtomDeexcitation * AtomDeexcitation()
void SetHighEnergyLimit(G4double)
virtual size_t NumberOfComponents(void) const
static constexpr double keV
G4ParticleChangeForGamma * fParticleChangeForGamma
const G4ThreeVector & GetMomentumDirection() const
const G4String & GetParticleName() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Hep3Vector & rotateUz(const Hep3Vector &)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
G4double GetPDGMass() const
G4double LowEnergyLimit() const
G4double IonisationEnergy(G4int level)
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
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)
G4DNACPA100WaterIonisationStructure waterStructure
G4VAtomDeexcitation * fAtomDeexcitation
static constexpr double m
void SetDeexcitationFlag(G4bool val)
G4ParticleDefinition * GetDefinition() const
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
const XML_Char int const XML_Char * value
static constexpr double electron_mass_c2
TriDimensionMap eNrjTransfData[6]
static G4DNAMolecularMaterial * Instance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
virtual ~G4DNACPA100IonisationModel()
static constexpr double twopi
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
static constexpr double eV
const G4Track * GetCurrentTrack() const
static G4Electron * Electron()
static G4Electron * ElectronDefinition()
std::vector< G4double > eTdummyVec
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
const std::vector< G4double > * fpMolWaterDensity
void SetLowEnergyLimit(G4double)
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
void RandomizeEjectedElectronDirection(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4double outgoingParticleEnergy, G4double &cosTheta, G4double &phi)
G4double GetKineticEnergy() const
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
virtual G4double FindValue(G4double e, G4int componentId=0) const
static constexpr double cm
G4GLOB_DLL std::ostream G4cout
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
G4double RandomizeEjectedElectronEnergyFromCompositionSampling(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
static G4LossTableManager * Instance()
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double Interpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4int RandomSelect(G4double energy, const G4String &particle)
G4double RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
G4double bindingEnergy(G4int A, G4int Z)
G4double HighEnergyLimit() const
static G4DNAChemistryManager * Instance()
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4double RandomTransferedEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
G4DNACPA100IonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNACPA100IonisationModel")