Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNABornIonisationModel2 Class Reference

#include <G4DNABornIonisationModel2.hh>

Inheritance diagram for G4DNABornIonisationModel2:

Public Member Functions

 G4DNABornIonisationModel2 (const G4ParticleDefinition *p=nullptr, const G4String &nam="DNABornIonisationModel")
 ~G4DNABornIonisationModel2 () override
G4DNABornIonisationModel2 & operator= (const G4DNABornIonisationModel2 &right)=delete
 G4DNABornIonisationModel2 (const G4DNABornIonisationModel2 &)=delete
void Initialise (const G4ParticleDefinition *, const G4DataVector &= *(new G4DataVector())) override
G4double CrossSectionPerVolume (const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax) override
void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4double GetPartialCrossSection (const G4Material *, G4int, const G4ParticleDefinition *, G4double) override
G4double DifferentialCrossSection (G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
G4double TransferedEnergy (G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, G4double random)
void SelectFasterComputation (G4bool input)
void SelectStationary (G4bool input)
void SelectSPScaling (G4bool input)
Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
virtual ~G4VEmModel ()
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual G4double ChargeSquareRatio (const G4Track &)
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual void StartTracking (G4Track *)
virtual void CorrectionsAlongStep (const G4Material *, const G4ParticleDefinition *, const G4double kinEnergy, const G4double cutEnergy, const G4double &length, G4double &eloss)
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual void DefineForRegion (const G4Region *)
virtual void FillNumberOfSecondaries (G4int &numberOfTriplets, G4int &numberOfRecoil)
virtual void ModelDescription (std::ostream &outFile) const
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
std::vector< G4EmElementSelector * > * GetElementSelectors ()
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element * SelectTargetAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element * SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element * GetCurrentElement (const G4Material *mat=nullptr) const
G4int SelectRandomAtomNumber (const G4Material *) const
const G4Isotope * GetCurrentIsotope (const G4Element *elm=nullptr) const
G4int SelectIsotopeNumber (const G4Element *) const
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
G4ElementData * GetElementData ()
G4PhysicsTable * GetCrossSectionTable ()
G4VEmFluctuationModel * GetModelOfFluctuations ()
G4VEmAngularDistribution * GetAngularDistribution ()
G4VEmModel * GetTripletModel ()
void SetTripletModel (G4VEmModel *)
void SetAngularDistribution (G4VEmAngularDistribution *)
G4double HighEnergyLimit () const
G4double LowEnergyLimit () const
G4double HighEnergyActivationLimit () const
G4double LowEnergyActivationLimit () const
G4double PolarAngleLimit () const
G4double SecondaryThreshold () const
G4bool DeexcitationFlag () const
G4bool ForceBuildTableFlag () const
G4bool UseAngularGeneratorFlag () const
void SetAngularGeneratorFlag (G4bool)
void SetHighEnergyLimit (G4double)
void SetLowEnergyLimit (G4double)
void SetActivationHighEnergyLimit (G4double)
void SetActivationLowEnergyLimit (G4double)
G4bool IsActive (G4double kinEnergy) const
void SetPolarAngleLimit (G4double)
void SetSecondaryThreshold (G4double)
void SetDeexcitationFlag (G4bool val)
void SetForceBuildTable (G4bool val)
void SetFluctuationFlag (G4bool val)
G4bool IsMaster () const
void SetUseBaseMaterials (G4bool val)
G4bool UseBaseMaterials () const
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
const G4String & GetName () const
void SetCurrentCouple (const G4MaterialCutsCouple *)
G4bool IsLocked () const
void SetLocked (G4bool)
void SetLPMFlag (G4bool)
void SetMasterThread (G4bool)
G4VEmModel & operator= (const G4VEmModel &right)=delete
 G4VEmModel (const G4VEmModel &)=delete

Protected Attributes

G4ParticleChangeForGamma * fParticleChangeForGamma
Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData = nullptr
G4VParticleChange * pParticleChange = nullptr
G4PhysicsTable * xSectionTable = nullptr
const G4Material * pBaseMaterial = nullptr
const std::vector< G4double > * theDensityFactor = nullptr
const std::vector< G4int > * theDensityIdx = nullptr
G4double inveplus
G4double pFactor = 1.0
std::size_t currentCoupleIndex = 0
std::size_t basedCoupleIndex = 0
G4bool lossFlucFlag = true

Additional Inherited Members

Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
const G4MaterialCutsCouple * CurrentCouple () const
void SetCurrentElement (const G4Element *)

Detailed Description

Definition at line 47 of file G4DNABornIonisationModel2.hh.

Constructor & Destructor Documentation

◆ G4DNABornIonisationModel2() [1/2]

G4DNABornIonisationModel2::G4DNABornIonisationModel2 ( const G4ParticleDefinition * p = nullptr,
const G4String & nam = "DNABornIonisationModel" )

Definition at line 46 of file G4DNABornIonisationModel2.cc.

47 :
48G4VEmModel(nam)
49{
50 verboseLevel = 0;
51 // Verbosity scale:
52 // 0 = nothing
53 // 1 = warning for energy non-conservation
54 // 2 = details of energy budget
55 // 3 = calculation of cross sections, file openings, sampling of atoms
56 // 4 = entering in methods
57
58 if (verboseLevel > 0)
59 {
60 G4cout << "Born ionisation model is constructed " << G4endl;
61 }
62
63 // Mark this model as "applicable" for atomic deexcitation
64
65 SetDeexcitationFlag(true);
66 fAtomDeexcitation = nullptr;
67 fParticleChangeForGamma = nullptr;
68 fpMolWaterDensity = nullptr;
69 fTableData = nullptr;
70 fLowEnergyLimit = 0;
71 fHighEnergyLimit = 0;
72 fParticleDef = nullptr;
73
74 // Define default angular generator
75
76 SetAngularDistribution(new G4DNABornAngle());
77
78 // Selection of computation method
79
80 fasterCode = false;
81
82 // Selection of stationary mode
83
84 statCode = false;
85
86 // Selection of SP scaling
87
88 spScaling = true;
89}
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4DNABornIonisationModel2::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition G4DNABornIonisationModel2.hh:91
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition G4VEmModel.hh:783
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition G4VEmModel.cc:67
G4VEmModel::SetAngularDistribution
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition G4VEmModel.hh:609

Referenced by G4DNABornIonisationModel2(), and operator=().

◆ ~G4DNABornIonisationModel2()

G4DNABornIonisationModel2::~G4DNABornIonisationModel2 ( )
override

Definition at line 93 of file G4DNABornIonisationModel2.cc.

94{
95 delete fTableData;
96 fVecm.clear();
97}

◆ G4DNABornIonisationModel2() [2/2]

G4DNABornIonisationModel2::G4DNABornIonisationModel2 ( const G4DNABornIonisationModel2 & )
delete

Member Function Documentation

◆ CrossSectionPerVolume()

G4double G4DNABornIonisationModel2::CrossSectionPerVolume ( const G4Material * material,
const G4ParticleDefinition * p,
G4double ekin,
G4double emin,
G4double emax )
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 259 of file G4DNABornIonisationModel2.cc.

264{
265 if (verboseLevel > 3)
266 {
267 G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel2"
268 << G4endl;
269 }
270
271 if (particleDefinition != fParticleDef) return 0;
272
273 // Calculate total cross section for model
274
275 G4double sigma=0;
276
277 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
278
279 if (ekin >= fLowEnergyLimit && ekin <= fHighEnergyLimit)
280 {
281 sigma = fTableData->FindValue(ekin);
282
283 // ICRU49 electronic SP scaling - ZF, SI
284
285 if (particleDefinition == G4Proton::ProtonDefinition() && ekin < 70*MeV && spScaling)
286 {
287 G4double A = 1.39241700556072800000E-009 ;
288 G4double B = -8.52610412942622630000E-002 ;
289 sigma = sigma * G4Exp(A*(ekin/eV)+B);
290 }
291 //
292 }
293
294 if (verboseLevel > 2)
295 {
296 G4cout << "__________________________________" << G4endl;
297 G4cout << "G4DNABornIonisationModel2 - XS INFO START" << G4endl;
298 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
299 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
300 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
301 G4cout << "G4DNABornIonisationModel2 - XS INFO END" << G4endl;
302 }
303
304 return sigma*waterDensity;
305}
B
G4double B(G4double temperature)
Definition G4DNAElectronHoleRecombination.cc:71
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:132
G4double
double G4double
Definition G4Types.hh:83
A
const G4double A[17]
Definition G4WaterStopping.cc:53
G4Material::GetIndex
std::size_t GetIndex() const
Definition G4Material.hh:239
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition G4Proton.cc:85

◆ DifferentialCrossSection()

G4double G4DNABornIonisationModel2::DifferentialCrossSection ( G4ParticleDefinition * aParticleDefinition,
G4double k,
G4double energyTransfer,
G4int shell )

Definition at line 611 of file G4DNABornIonisationModel2.cc.

615{
616 G4double sigma = 0.;
617
618 if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
619 {
620 G4double valueT1 = 0;
621 G4double valueT2 = 0;
622 G4double valueE21 = 0;
623 G4double valueE22 = 0;
624 G4double valueE12 = 0;
625 G4double valueE11 = 0;
626
627 G4double xs11 = 0;
628 G4double xs12 = 0;
629 G4double xs21 = 0;
630 G4double xs22 = 0;
631
632 // Protection against out of boundary access - proton case : 100 MeV
633 if (k==fTdummyVec.back()) k=k*(1.-1e-12);
634 //
635
636 // k should be in eV and energy transfer eV also
637
638 auto t2 = std::upper_bound(fTdummyVec.begin(),
639 fTdummyVec.end(),
640 k);
641
642 auto t1 = t2 - 1;
643
644 // SI : the following condition avoids situations where energyTransfer >last vector element
645
646 if (energyTransfer <= fVecm[(*t1)].back()
647 && energyTransfer <= fVecm[(*t2)].back())
648 {
649 auto e12 = std::upper_bound(fVecm[(*t1)].begin(),
650 fVecm[(*t1)].end(),
651 energyTransfer);
652 auto e11 = e12 - 1;
653
654 auto e22 = std::upper_bound(fVecm[(*t2)].begin(),
655 fVecm[(*t2)].end(),
656 energyTransfer);
657 auto e21 = e22 - 1;
658
659 valueT1 = *t1;
660 valueT2 = *t2;
661 valueE21 = *e21;
662 valueE22 = *e22;
663 valueE12 = *e12;
664 valueE11 = *e11;
665
666 xs11 = fDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
667 xs12 = fDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
668 xs21 = fDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
669 xs22 = fDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
670
671 }
672
673 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
674 if (xsProduct != 0.)
675 {
676 sigma = QuadInterpolator(valueE11,
677 valueE12,
678 valueE21,
679 valueE22,
680 xs11,
681 xs12,
682 xs21,
683 xs22,
684 valueT1,
685 valueT2,
686 k,
687 energyTransfer);
688 }
689 }
690
691 return sigma;
692}

◆ GetPartialCrossSection()

G4double G4DNABornIonisationModel2::GetPartialCrossSection ( const G4Material * ,
G4int level,
const G4ParticleDefinition * ,
G4double kineticEnergy )
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 786 of file G4DNABornIonisationModel2.cc.

790{
791 return fTableData->GetComponent(level)->FindValue(kineticEnergy);
792}

◆ Initialise()

void G4DNABornIonisationModel2::Initialise ( const G4ParticleDefinition * particle,
const G4DataVector & = *(new G4DataVector()) )
overridevirtual

Implements G4VEmModel.

Definition at line 101 of file G4DNABornIonisationModel2.cc.

103{
104
105 if (verboseLevel > 3)
106 {
107 G4cout << "Calling G4DNABornIonisationModel2::Initialise()" << G4endl;
108 }
109
110 if(fParticleDef != nullptr && particle != fParticleDef)
111 {
112 G4ExceptionDescription description;
113 description << "You are trying to initialized G4DNABornIonisationModel2 "
114 "for particle "
115 << particle->GetParticleName()
116 << G4endl;
117 description << "G4DNABornIonisationModel2 was already initialised "
118 "for particle:" << fParticleDef->GetParticleName() << G4endl;
119 G4Exception("G4DNABornIonisationModel2::Initialise","bornIonInit",
120 FatalException,description);
121 }
122
123 fParticleDef = particle;
124
125 // Energy limits
126 const char* path = G4FindDataDir("G4LEDATA");
127
128 // ***
129
130 G4String particleName = particle->GetParticleName();
131 std::ostringstream fullFileName;
132 fullFileName << path;
133
134 if(particleName == "e-")
135 {
136 fTableFile = "dna/sigma_ionisation_e_born";
137 fLowEnergyLimit = 11.*eV;
138 fHighEnergyLimit = 1.*MeV;
139
140 if (fasterCode)
141 {
142 fullFileName << "/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
143 }
144 else
145 {
146 fullFileName << "/dna/sigmadiff_ionisation_e_born.dat";
147 }
148 }
149 else if(particleName == "proton")
150 {
151 fTableFile = "dna/sigma_ionisation_p_born";
152 fLowEnergyLimit = 500. * keV;
153 fHighEnergyLimit = 100. * MeV;
154
155 if (fasterCode)
156 {
157 fullFileName << "/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
158 }
159 else
160 {
161 fullFileName << "/dna/sigmadiff_ionisation_p_born.dat";
162 }
163 }
164
165 // Cross section
166
167 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
168 fTableData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
169 fTableData->LoadData(fTableFile);
170
171 // Final state
172
173 std::ifstream diffCrossSection(fullFileName.str().c_str());
174
175 if (!diffCrossSection)
176 {
177 G4ExceptionDescription description;
178 description << "Missing data file:" << G4endl << fullFileName.str() << G4endl;
179 G4Exception("G4DNABornIonisationModel2::Initialise","em0003",
180 FatalException,description);
181 }
182
183 // Clear the arrays for re-initialization case (MT mode)
184 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
185
186 fTdummyVec.clear();
187 fVecm.clear();
188
189 for (int j=0; j<5; j++)
190 {
191 fProbaShellMap[j].clear();
192 fDiffCrossSectionData[j].clear();
193 fNrjTransfData[j].clear();
194 }
195
196 //
197
198 fTdummyVec.push_back(0.);
199 while(!diffCrossSection.eof())
200 {
201 G4double tDummy;
202 G4double eDummy;
203 diffCrossSection>>tDummy>>eDummy;
204 if (tDummy != fTdummyVec.back()) fTdummyVec.push_back(tDummy);
205
206 G4double tmp;
207 for (int j=0; j<5; j++)
208 {
209 diffCrossSection>> tmp;
210
211 fDiffCrossSectionData[j][tDummy][eDummy] = tmp;
212
213 if (fasterCode)
214 {
215 fNrjTransfData[j][tDummy][fDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
216 fProbaShellMap[j][tDummy].push_back(fDiffCrossSectionData[j][tDummy][eDummy]);
217 }
218
219 // SI - only if eof is not reached
220 if (!diffCrossSection.eof() && !fasterCode) fDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
221
222 if (!fasterCode) fVecm[tDummy].push_back(eDummy);
223
224 }
225 }
226
227 //
228 SetLowEnergyLimit(fLowEnergyLimit);
229 SetHighEnergyLimit(fHighEnergyLimit);
230
231 if( verboseLevel>0 )
232 {
233 G4cout << "Born ionisation model is initialized " << G4endl
234 << "Energy range: "
235 << LowEnergyLimit() / eV << " eV - "
236 << HighEnergyLimit() / keV << " keV for "
237 << particle->GetParticleName()
238 << G4endl;
239 }
240
241 // Initialize water density pointer
242
243 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
244 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
245
246 if (isInitialised)
247 { return;}
248 fParticleChangeForGamma = GetParticleChangeForGamma();
249 statCode = G4EmParameters::Instance()->DNAStationary();
250
251 if (!statCode)
252 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
253
254 isInitialised = true;
255}
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition G4FindDataDir.cc:94
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition G4Exception.hh:40
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition G4DNAMolecularMaterial.cc:80
G4EmParameters::Instance
static G4EmParameters * Instance()
Definition G4EmParameters.cc:71
G4EmParameters::DNAStationary
G4bool DNAStationary() const
Definition G4EmParameters.cc:451
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition G4LossTableManager.cc:90
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition G4LossTableManager.hh:396
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition G4Material.cc:663
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition G4ParticleDefinition.hh:96
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition G4VEmModel.hh:734
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition G4VEmModel.cc:121
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition G4VEmModel.hh:643
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition G4VEmModel.hh:636
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition G4VEmModel.hh:741

◆ operator=()

G4DNABornIonisationModel2 & G4DNABornIonisationModel2::operator= ( const G4DNABornIonisationModel2 & right)
delete

◆ SampleSecondaries()

void G4DNABornIonisationModel2::SampleSecondaries ( std::vector< G4DynamicParticle * > * fvect,
const G4MaterialCutsCouple * couple,
const G4DynamicParticle * particle,
G4double tmin,
G4double maxEnergy )
overridevirtual

Implements G4VEmModel.

Definition at line 309 of file G4DNABornIonisationModel2.cc.

314{
315
316 if (verboseLevel > 3)
317 {
318 G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel2"
319 << G4endl;
320 }
321
322 G4double k = particle->GetKineticEnergy();
323
324 if (k >= fLowEnergyLimit && k <= fHighEnergyLimit)
325 {
326 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
327 G4double particleMass = particle->GetDefinition()->GetPDGMass();
328 G4double totalEnergy = k + particleMass;
329 G4double pSquare = k * (totalEnergy + particleMass);
330 G4double totalMomentum = std::sqrt(pSquare);
331
332 G4int ionizationShell = 0;
333
334 if (!fasterCode) ionizationShell = RandomSelect(k);
335
336 // SI: The following protection is necessary to avoid infinite loops :
337 // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
338 // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
339 // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
340 // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
341
342 if (fasterCode)
343 do
344 {
345 ionizationShell = RandomSelect(k);
346 } while (k<19*eV && ionizationShell==2 && particle->GetDefinition()==G4Electron::ElectronDefinition());
347
348 G4double secondaryKinetic=-1000*eV;
349
350 if (!fasterCode)
351 {
352 secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
353 }
354 else
355 {
356 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
357 }
358
359 G4int Z = 8;
360
361 G4ThreeVector deltaDirection =
362 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
363 Z, ionizationShell,
364 couple->GetMaterial());
365
366 if (secondaryKinetic>0)
367 {
368 auto dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
369 fvect->push_back(dp);
370 }
371
372 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
373 {
374 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
375
376 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
377 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
378 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
379 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
380 finalPx /= finalMomentum;
381 finalPy /= finalMomentum;
382 finalPz /= finalMomentum;
383
384 G4ThreeVector direction;
385 direction.set(finalPx,finalPy,finalPz);
386
387 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
388 }
389
390 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
391
392 // AM: sample deexcitation
393 // here we assume that H_{2}O electronic levels are the same as Oxygen.
394 // this can be considered true with a rough 10% error in energy on K-shell,
395
396 std::size_t secNumberInit = 0;
397 std::size_t secNumberFinal = 0;
398
399 G4double bindingEnergy = 0;
400 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
401
402 // SI: additional protection if tcs interpolation method is modified
403 if (k<bindingEnergy) return;
404 //
405
406 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
407
408 // SI: only atomic deexcitation from K shell is considered
409 if((fAtomDeexcitation != nullptr) && ionizationShell == 4)
410 {
411 const G4AtomicShell* shell =
412 fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
413 secNumberInit = fvect->size();
414 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
415 secNumberFinal = fvect->size();
416
417 if(secNumberFinal > secNumberInit)
418 {
419 for (std::size_t i=secNumberInit; i<secNumberFinal; ++i)
420 {
421 //Check if there is enough residual energy
422 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
423 {
424 //Ok, this is a valid secondary: keep it
425 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
426 }
427 else
428 {
429 //Invalid secondary: not enough energy to create it!
430 //Keep its energy in the local deposit
431 delete (*fvect)[i];
432 (*fvect)[i]=nullptr;
433 }
434 }
435 }
436
437 }
438
439 //This should never happen
440 if(bindingEnergy < 0.0)
441 G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
442 "em2050",FatalException,"Negative local energy deposit");
443
444 //bindingEnergy has been decreased
445 //by the amount of energy taken away by deexc. products
446 if (!statCode)
447 {
448 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
449 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
450 }
451 else
452 {
453 fParticleChangeForGamma->SetProposedKineticEnergy(k);
454 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
455 }
456
457 // TEST //////////////////////////
458 // if (secondaryKinetic<0) abort();
459 // if (scatteredEnergy<0) abort();
460 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
461 // if (k-scatteredEnergy<0) abort();
462 /////////////////////////////////
463
464 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
465 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
466 ionizationShell,
467 theIncomingTrack);
468 }
469}
eIonizedMolecule
@ eIonizedMolecule
Definition G4DNAChemistryManager.hh:73
G4ParticleMomentum
G4ThreeVector G4ParticleMomentum
Definition G4ParticleMomentum.hh:39
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4int
int G4int
Definition G4Types.hh:85
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4DNAChemistryManager::Instance
static G4DNAChemistryManager * Instance()
Definition G4DNAChemistryManager.cc:125
G4DNAChemistryManager::CreateWaterMolecule
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Definition G4DNAChemistryManager.cc:600
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition G4Electron.cc:86
G4Electron::Electron
static G4Electron * Electron()
Definition G4Electron.cc:91
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition G4MaterialCutsCouple.hh:166
G4VEmAngularDistribution::SampleDirectionForShell
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
Definition G4VEmAngularDistribution.cc:68
G4VEmModel::GetAngularDistribution
G4VEmAngularDistribution * GetAngularDistribution()
Definition G4VEmModel.hh:602
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition bindingEnergy.cc:36

◆ SelectFasterComputation()

void G4DNABornIonisationModel2::SelectFasterComputation ( G4bool input)
inline

Definition at line 160 of file G4DNABornIonisationModel2.hh.

161{
162 fasterCode = input;
163}

◆ SelectSPScaling()

void G4DNABornIonisationModel2::SelectSPScaling ( G4bool input)
inline

Definition at line 174 of file G4DNABornIonisationModel2.hh.

175{
176 spScaling = input;
177}

◆ SelectStationary()

void G4DNABornIonisationModel2::SelectStationary ( G4bool input)
inline

Definition at line 167 of file G4DNABornIonisationModel2.hh.

168{
169 statCode = input;
170}

◆ TransferedEnergy()

G4double G4DNABornIonisationModel2::TransferedEnergy ( G4ParticleDefinition * aParticleDefinition,
G4double incomingParticleEnergy,
G4int shell,
G4double random )

Definition at line 862 of file G4DNABornIonisationModel2.cc.

866{
867
868 G4double nrj = 0.;
869
870 G4double valueK1 = 0;
871 G4double valueK2 = 0;
872 G4double valuePROB21 = 0;
873 G4double valuePROB22 = 0;
874 G4double valuePROB12 = 0;
875 G4double valuePROB11 = 0;
876
877 G4double nrjTransf11 = 0;
878 G4double nrjTransf12 = 0;
879 G4double nrjTransf21 = 0;
880 G4double nrjTransf22 = 0;
881
882 // Protection against out of boundary access - proton case : 100 MeV
883 if (k==fTdummyVec.back()) k=k*(1.-1e-12);
884 //
885
886 // k should be in eV
887 auto k2 = std::upper_bound(fTdummyVec.begin(),
888 fTdummyVec.end(),
889 k);
890 auto k1 = k2 - 1;
891
892 /*
893 G4cout << "----> k=" << k
894 << " " << *k1
895 << " " << *k2
896 << " " << random
897 << " " << ionizationLevelIndex
898 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
899 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
900 << G4endl;
901 */
902
903 // SI : the following condition avoids situations where random >last vector element
904 if (random <= fProbaShellMap[ionizationLevelIndex][(*k1)].back()
905 && random <= fProbaShellMap[ionizationLevelIndex][(*k2)].back())
906 {
907 auto prob12 =
908 std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
909 fProbaShellMap[ionizationLevelIndex][(*k1)].end(),
910 random);
911
912 auto prob11 = prob12 - 1;
913
914 auto prob22 =
915 std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
916 fProbaShellMap[ionizationLevelIndex][(*k2)].end(),
917 random);
918
919 auto prob21 = prob22 - 1;
920
921 valueK1 = *k1;
922 valueK2 = *k2;
923 valuePROB21 = *prob21;
924 valuePROB22 = *prob22;
925 valuePROB12 = *prob12;
926 valuePROB11 = *prob11;
927
928 /*
929 G4cout << " " << random << " " << valuePROB11 << " "
930 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
931 */
932
933 nrjTransf11 = fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
934 nrjTransf12 = fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
935 nrjTransf21 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
936 nrjTransf22 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
937
938 /*
939 G4cout << " " << ionizationLevelIndex << " "
940 << random << " " <<valueK1 << " " << valueK2 << G4endl;
941
942 G4cout << " " << random << " " << nrjTransf11 << " "
943 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
944 */
945
946 }
947 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
948 if (random > fProbaShellMap[ionizationLevelIndex][(*k1)].back())
949 {
950 auto prob22 =
951 std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
952 fProbaShellMap[ionizationLevelIndex][(*k2)].end(),
953 random);
954
955 auto prob21 = prob22 - 1;
956
957 valueK1 = *k1;
958 valueK2 = *k2;
959 valuePROB21 = *prob21;
960 valuePROB22 = *prob22;
961
962 // G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
963
964 nrjTransf21 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
965 nrjTransf22 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
966
967 G4double interpolatedvalue2 = Interpolate(valuePROB21,
968 valuePROB22,
969 random,
970 nrjTransf21,
971 nrjTransf22);
972
973 // zeros are explicitly set
974
975 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
976
977 /*
978 G4cout << " " << ionizationLevelIndex << " "
979 << random << " " <<valueK1 << " " << valueK2 << G4endl;
980
981 G4cout << " " << random << " " << nrjTransf11 << " "
982 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
983
984 G4cout << "ici" << " " << value << G4endl;
985 */
986
987 return value;
988 }
989
990 // End electron and proton cases
991
992 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
993 * nrjTransf22;
994 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
995
996 if (nrjTransfProduct != 0.)
997 {
998 nrj = QuadInterpolator(valuePROB11,
999 valuePROB12,
1000 valuePROB21,
1001 valuePROB22,
1002 nrjTransf11,
1003 nrjTransf12,
1004 nrjTransf21,
1005 nrjTransf22,
1006 valueK1,
1007 valueK2,
1008 k,
1009 random);
1010 }
1011 // G4cout << nrj << endl;
1012
1013 return nrj;
1014}

Member Data Documentation

◆ fParticleChangeForGamma

G4ParticleChangeForGamma* G4DNABornIonisationModel2::fParticleChangeForGamma
protected

Definition at line 91 of file G4DNABornIonisationModel2.hh.

Referenced by G4DNABornIonisationModel2(), Initialise(), and SampleSecondaries().

The documentation for this class was generated from the following files: