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

#include <G4DNABornIonisationModel1.hh>

Inheritance diagram for G4DNABornIonisationModel1:

Public Member Functions

 G4DNABornIonisationModel1 (const G4ParticleDefinition *p=nullptr, const G4String &nam="DNABornIonisationModel")
 ~G4DNABornIonisationModel1 () override
G4DNABornIonisationModel1 & operator= (const G4DNABornIonisationModel1 &right)=delete
 G4DNABornIonisationModel1 (const G4DNABornIonisationModel1 &)=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
void StartTracking (G4Track *) 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 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 48 of file G4DNABornIonisationModel1.hh.

Constructor & Destructor Documentation

◆ G4DNABornIonisationModel1() [1/2]

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

Definition at line 46 of file G4DNABornIonisationModel1.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 SetDeexcitationFlag(true);
65 fAtomDeexcitation = nullptr;
66 fParticleChangeForGamma = nullptr;
67 fpMolWaterDensity = nullptr;
68
69 // Define default angular generator
70 SetAngularDistribution(new G4DNABornAngle());
71
72 fasterCode = G4EmParameters::Instance()->DNAFast();
73}
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4DNABornIonisationModel1::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition G4DNABornIonisationModel1.hh:94
G4EmParameters::Instance
static G4EmParameters * Instance()
Definition G4EmParameters.cc:71
G4EmParameters::DNAFast
G4bool DNAFast() const
Definition G4EmParameters.cc:439
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 G4DNABornIonisationModel1(), and operator=().

◆ ~G4DNABornIonisationModel1()

G4DNABornIonisationModel1::~G4DNABornIonisationModel1 ( )
override

Definition at line 77 of file G4DNABornIonisationModel1.cc.

78{
79 // Cross section
80
81 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
82 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
83 {
84 G4DNACrossSectionDataSet* table = pos->second;
85 delete table;
86 }
87
88 // Final state
89
90 eVecm.clear();
91 pVecm.clear();
92}

◆ G4DNABornIonisationModel1() [2/2]

G4DNABornIonisationModel1::G4DNABornIonisationModel1 ( const G4DNABornIonisationModel1 & )
delete

Member Function Documentation

◆ CrossSectionPerVolume()

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

Reimplemented from G4VEmModel.

Definition at line 319 of file G4DNABornIonisationModel1.cc.

324{
325 if (verboseLevel > 3)
326 {
327 G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel1"
328 << G4endl;
329 }
330
331 if (
332 particleDefinition != G4Proton::ProtonDefinition()
333 &&
334 particleDefinition != G4Electron::ElectronDefinition()
335 )
336
337 return 0;
338
339 // Calculate total cross section for model
340
341 G4double lowLim = 0;
342 G4double highLim = 0;
343 G4double sigma=0;
344
345 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
346
347 const G4String& particleName = particleDefinition->GetParticleName();
348
349 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
350 pos1 = lowEnergyLimit.find(particleName);
351 if (pos1 != lowEnergyLimit.end())
352 {
353 lowLim = pos1->second;
354 }
355
356 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
357 pos2 = highEnergyLimit.find(particleName);
358 if (pos2 != highEnergyLimit.end())
359 {
360 highLim = pos2->second;
361 }
362
363 if (ekin >= lowLim && ekin <= highLim)
364 {
365 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
366 pos = tableData.find(particleName);
367
368 if (pos != tableData.end())
369 {
370 G4DNACrossSectionDataSet* table = pos->second;
371 if (table != nullptr)
372 {
373 sigma = table->FindValue(ekin);
374
375 // ICRU49 electronic SP scaling - ZF, SI
376
377 if (particleDefinition == G4Proton::ProtonDefinition() && ekin < 70*MeV && spScaling)
378 {
379 G4double A = 1.39241700556072800000E-009 ;
380 G4double B = -8.52610412942622630000E-002 ;
381 sigma = sigma * G4Exp(A*(ekin/eV)+B);
382 }
383 //
384
385 }
386 }
387 else
388 {
389 G4Exception("G4DNABornIonisationModel1::CrossSectionPerVolume","em0002",
390 FatalException,"Model not applicable to particle type.");
391 }
392 }
393
394 if (verboseLevel > 2)
395 {
396 G4cout << "__________________________________" << G4endl;
397 G4cout << "G4DNABornIonisationModel1 - XS INFO START" << G4endl;
398 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
399 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
400 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
401 G4cout << "G4DNABornIonisationModel1 - XS INFO END" << G4endl;
402 }
403
404 return sigma*waterDensity;
405}
B
G4double B(G4double temperature)
Definition G4DNAElectronHoleRecombination.cc:71
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
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
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition G4Electron.cc:86
G4Material::GetIndex
std::size_t GetIndex() const
Definition G4Material.hh:239
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition G4Proton.cc:85

◆ DifferentialCrossSection()

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

Definition at line 692 of file G4DNABornIonisationModel1.cc.

696{
697 G4double sigma = 0.;
698
699 if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
700 {
701 G4double valueT1 = 0;
702 G4double valueT2 = 0;
703 G4double valueE21 = 0;
704 G4double valueE22 = 0;
705 G4double valueE12 = 0;
706 G4double valueE11 = 0;
707
708 G4double xs11 = 0;
709 G4double xs12 = 0;
710 G4double xs21 = 0;
711 G4double xs22 = 0;
712
713 if (particleDefinition == G4Electron::ElectronDefinition())
714 {
715
716 // Protection against out of boundary access - electron case : 1 MeV
717 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
718 //
719
720 // k should be in eV and energy transfer eV also
721
722 auto t2 = std::upper_bound(eTdummyVec.begin(),
723 eTdummyVec.end(),
724 k);
725
726 auto t1 = t2 - 1;
727
728 // SI : the following condition avoids situations where energyTransfer >last vector element
729 if (energyTransfer <= eVecm[(*t1)].back()
730 && energyTransfer <= eVecm[(*t2)].back())
731 {
732 auto e12 =
733 std::upper_bound(eVecm[(*t1)].begin(),
734 eVecm[(*t1)].end(),
735 energyTransfer);
736 auto e11 = e12 - 1;
737
738 auto e22 =
739 std::upper_bound(eVecm[(*t2)].begin(),
740 eVecm[(*t2)].end(),
741 energyTransfer);
742 auto e21 = e22 - 1;
743
744 valueT1 = *t1;
745 valueT2 = *t2;
746 valueE21 = *e21;
747 valueE22 = *e22;
748 valueE12 = *e12;
749 valueE11 = *e11;
750
751 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
752 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
753 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
754 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
755
756 }
757
758 }
759
760 if (particleDefinition == G4Proton::ProtonDefinition())
761 {
762 // Protection against out of boundary access - proton case : 100 MeV
763 if (k==pTdummyVec.back()) k=k*(1.-1e-12);
764 //
765
766 // k should be in eV and energy transfer eV also
767 auto t2 = std::upper_bound(pTdummyVec.begin(),
768 pTdummyVec.end(),
769 k);
770 auto t1 = t2 - 1;
771
772 auto e12 = std::upper_bound(pVecm[(*t1)].begin(),
773 pVecm[(*t1)].end(),
774 energyTransfer);
775 auto e11 = e12 - 1;
776
777 auto e22 = std::upper_bound(pVecm[(*t2)].begin(),
778 pVecm[(*t2)].end(),
779 energyTransfer);
780 auto e21 = e22 - 1;
781
782 valueT1 = *t1;
783 valueT2 = *t2;
784 valueE21 = *e21;
785 valueE22 = *e22;
786 valueE12 = *e12;
787 valueE11 = *e11;
788
789 xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
790 xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
791 xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
792 xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
793
794 }
795
796 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
797 if (xsProduct != 0.)
798 {
799 sigma = QuadInterpolator(valueE11,
800 valueE12,
801 valueE21,
802 valueE22,
803 xs11,
804 xs12,
805 xs21,
806 xs22,
807 valueT1,
808 valueT2,
809 k,
810 energyTransfer);
811 }
812 }
813
814 return sigma;
815}

◆ GetPartialCrossSection()

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

Reimplemented from G4VEmModel.

Definition at line 907 of file G4DNABornIonisationModel1.cc.

911{
912 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
913 pos = tableData.find(particle->GetParticleName());
914
915 if (pos != tableData.end())
916 {
917 G4DNACrossSectionDataSet* table = pos->second;
918 return table->GetComponent(level)->FindValue(kineticEnergy);
919 }
920
921 return 0;
922}
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition G4ParticleDefinition.hh:96

◆ Initialise()

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

Implements G4VEmModel.

Definition at line 96 of file G4DNABornIonisationModel1.cc.

98{
99
100 if (verboseLevel > 3)
101 {
102 G4cout << "Calling G4DNABornIonisationModel1::Initialise()" << G4endl;
103 }
104
105 // Energy limits
106
107 G4String fileElectron("dna/sigma_ionisation_e_born");
108 G4String fileProton("dna/sigma_ionisation_p_born");
109
110 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
111 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
112
113 G4String electron;
114 G4String proton;
115
116 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
117
118 const char *path = G4FindDataDir("G4LEDATA");
119
120 // *** ELECTRON
121
122 electron = electronDef->GetParticleName();
123
124 tableFile[electron] = fileElectron;
125
126 lowEnergyLimit[electron] = 11. * eV;
127 highEnergyLimit[electron] = 1. * MeV;
128
129 // Cross section
130
131 auto tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
132 tableE->LoadData(fileElectron);
133
134 tableData[electron] = tableE;
135
136 // Final state
137
138 std::ostringstream eFullFileName;
139
140 if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
141 if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
142
143 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
144
145 if (!eDiffCrossSection)
146 {
147 if (fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
148 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat");
149
150 if (!fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
151 FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
152 }
153
154 // Clear the arrays for re-initialization case (MT mode)
155 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
156
157 eTdummyVec.clear();
158 pTdummyVec.clear();
159
160 eVecm.clear();
161 pVecm.clear();
162
163 for (G4int j=0; j<5; j++)
164 {
165 eProbaShellMap[j].clear();
166 pProbaShellMap[j].clear();
167
168 eDiffCrossSectionData[j].clear();
169 pDiffCrossSectionData[j].clear();
170
171 eNrjTransfData[j].clear();
172 pNrjTransfData[j].clear();
173 }
174
175 //
176
177 eTdummyVec.push_back(0.);
178 while(!eDiffCrossSection.eof())
179 {
180 G4double tDummy;
181 G4double eDummy;
182 eDiffCrossSection>>tDummy>>eDummy;
183 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
184
185 G4double tmp;
186 for (G4int j=0; j<5; j++)
187 {
188 eDiffCrossSection>> tmp;
189
190 eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
191
192 if (fasterCode)
193 {
194 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
195 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
196 }
197
198 // SI - only if eof is not reached
199 if (!eDiffCrossSection.eof() && !fasterCode) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
200
201 if (!fasterCode) eVecm[tDummy].push_back(eDummy);
202
203 }
204 }
205
206 // *** PROTON
207
208 proton = protonDef->GetParticleName();
209
210 tableFile[proton] = fileProton;
211
212 lowEnergyLimit[proton] = 500. * keV;
213 highEnergyLimit[proton] = 100. * MeV;
214
215 // Cross section
216
217 auto tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
218 tableP->LoadData(fileProton);
219
220 tableData[proton] = tableP;
221
222 // Final state
223
224 std::ostringstream pFullFileName;
225
226 if (fasterCode) pFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
227
228 if (!fasterCode) pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
229
230 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
231
232 if (!pDiffCrossSection)
233 {
234 if (fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
235 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat");
236
237 if (!fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
238 FatalException,"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
239 }
240
241 pTdummyVec.push_back(0.);
242 while(!pDiffCrossSection.eof())
243 {
244 G4double tDummy;
245 G4double eDummy;
246 pDiffCrossSection>>tDummy>>eDummy;
247 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
248 for (G4int j=0; j<5; j++)
249 {
250 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
251
252 if (fasterCode)
253 {
254 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
255 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
256 }
257
258 // SI - only if eof is not reached !
259 if (!pDiffCrossSection.eof() && !fasterCode) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
260
261 if (!fasterCode) pVecm[tDummy].push_back(eDummy);
262 }
263 }
264
265 //
266
267 if (particle==electronDef)
268 {
269 SetLowEnergyLimit(lowEnergyLimit[electron]);
270 SetHighEnergyLimit(highEnergyLimit[electron]);
271 }
272
273 if (particle==protonDef)
274 {
275 SetLowEnergyLimit(lowEnergyLimit[proton]);
276 SetHighEnergyLimit(highEnergyLimit[proton]);
277 }
278
279 if( verboseLevel>0 )
280 {
281 G4cout << "Born ionisation model is initialized " << G4endl
282 << "Energy range: "
283 << LowEnergyLimit() / eV << " eV - "
284 << HighEnergyLimit() / keV << " keV for "
285 << particle->GetParticleName()
286 << G4endl;
287 }
288
289 if (isInitialised) { return; }
290 fParticleChangeForGamma = GetParticleChangeForGamma();
291
292 // Initialize water density pointer
293 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
294 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
295
296 statCode = G4EmParameters::Instance()->DNAStationary();
297
298 // AD
299 if (!statCode)
300 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
301
302 // chemistry
303 auto chem = G4DNAChemistryManager::Instance();
304 if (chem->IsChemistryActivated()) {
305 fChemistry = chem;
306 }
307 isInitialised = true;
308}
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition G4FindDataDir.cc:94
G4int
int G4int
Definition G4Types.hh:85
G4DNAChemistryManager::Instance
static G4DNAChemistryManager * Instance()
Definition G4DNAChemistryManager.cc:125
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition G4DNAMolecularMaterial.cc:80
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
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=()

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

◆ SampleSecondaries()

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

Implements G4VEmModel.

Definition at line 409 of file G4DNABornIonisationModel1.cc.

414{
415
416 if (verboseLevel > 3)
417 {
418 G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel1"
419 << G4endl;
420 }
421
422 G4double lowLim = 0;
423 G4double highLim = 0;
424
425 G4double k = particle->GetKineticEnergy();
426
427 const G4String& particleName = particle->GetDefinition()->GetParticleName();
428
429 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
430 pos1 = lowEnergyLimit.find(particleName);
431
432 if (pos1 != lowEnergyLimit.end())
433 {
434 lowLim = pos1->second;
435 }
436
437 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
438 pos2 = highEnergyLimit.find(particleName);
439
440 if (pos2 != highEnergyLimit.end())
441 {
442 highLim = pos2->second;
443 }
444
445 if (k >= lowLim && k <= highLim)
446 {
447 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
448 G4double particleMass = particle->GetDefinition()->GetPDGMass();
449 G4double totalEnergy = k + particleMass;
450 G4double pSquare = k * (totalEnergy + particleMass);
451 G4double totalMomentum = std::sqrt(pSquare);
452
453 G4int ionizationShell = 0;
454
455 // SI: The following protection is necessary to avoid infinite loops :
456 // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
457 // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
458 // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
459 // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
460
461 if (!fasterCode) {
462 ionizationShell = RandomSelect(k,particleName);
463 } else {
464 do {
465 ionizationShell = RandomSelect(k,particleName);
466 } while (k<19*eV && ionizationShell==2 && particle->GetDefinition()==G4Electron::ElectronDefinition());
467 }
468 G4double bindingEnergy = 0;
469 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
470
471 // SI: additional protection if tcs interpolation method is modified
472 if (k<bindingEnergy) return;
473 //
474
475 G4double secondaryKinetic=-1000*eV;
476
477 if (!fasterCode)
478 {
479 secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
480 }
481 else
482 {
483 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
484 }
485 G4int Z = 8;
486
487 G4ThreeVector deltaDirection =
488 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
489 Z, ionizationShell,
490 couple->GetMaterial());
491
492 if (secondaryKinetic>0)
493 {
494 auto dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
495 fvect->push_back(dp);
496 }
497
498 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
499 {
500 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
501
502 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
503 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
504 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
505 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
506 finalPx /= finalMomentum;
507 finalPy /= finalMomentum;
508 finalPz /= finalMomentum;
509
510 G4ThreeVector direction;
511 direction.set(finalPx,finalPy,finalPz);
512
513 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
514 }
515
516 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
517
518 // AM: sample deexcitation
519 // here we assume that H_{2}O electronic levels are the same as Oxygen.
520 // this can be considered true with a rough 10% error in energy on K-shell,
521
522 std::size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
523 std::size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
524
525 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
526
527 // SI: only atomic deexcitation from K shell is considered
528 if((fAtomDeexcitation != nullptr) && ionizationShell == 4)
529 {
530 const G4AtomicShell* shell =
531 fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
532 secNumberInit = fvect->size();
533 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
534 secNumberFinal = fvect->size();
535
536 //TEST
537 //G4cout << "ionizationShell=" << ionizationShell<< G4endl;
538 //G4cout << "bindingEnergy=" << bindingEnergy/eV<< G4endl;
539
540 if(secNumberFinal > secNumberInit)
541 {
542 for (std::size_t i=secNumberInit; i<secNumberFinal; ++i)
543 {
544 //Check if there is enough residual energy
545 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
546 {
547 //Ok, this is a valid secondary: keep it
548 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
549 //G4cout << "--deex nrj=" << ((*fvect)[i])->GetKineticEnergy()/eV
550 //<< G4endl;
551 }
552 else
553 {
554 //Invalid secondary: not enough energy to create it!
555 //Keep its energy in the local deposit
556 delete (*fvect)[i];
557 (*fvect)[i]=nullptr;
558 }
559 }
560 }
561
562 //TEST
563 //G4cout << "k=" << k/eV<< G4endl;
564 //G4cout << "secondaryKinetic=" << secondaryKinetic/eV<< G4endl;
565 //G4cout << "scatteredEnergy=" << scatteredEnergy/eV<< G4endl;
566 //G4cout << "deposited energy=" << bindingEnergy/eV<< G4endl;
567 //
568
569 }
570
571 //This should never happen
572 if(bindingEnergy < 0.0)
573 G4Exception("G4DNABornIonisatioModel1::SampleSecondaries()",
574 "em2050",FatalException,"Negative local energy deposit");
575
576 //bindingEnergy has been decreased
577 //by the amount of energy taken away by deexc. products
578 if (!statCode)
579 {
580 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
581 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
582 }
583 else
584 {
585 fParticleChangeForGamma->SetProposedKineticEnergy(k);
586 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
587 }
588
589 // create radical
590 if (nullptr != fChemistry) {
591 fChemistry->CreateWaterMolecule(eIonizedMolecule, ionizationShell, fTrack);
592 }
593 }
594}
eIonizedMolecule
@ eIonizedMolecule
Definition G4DNAChemistryManager.hh:73
G4ParticleMomentum
G4ThreeVector G4ParticleMomentum
Definition G4ParticleMomentum.hh:39
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
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)
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
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 G4DNABornIonisationModel1::SelectFasterComputation ( G4bool input)
inline

Definition at line 171 of file G4DNABornIonisationModel1.hh.

172{
173 fasterCode = input;
174}

◆ SelectSPScaling()

void G4DNABornIonisationModel1::SelectSPScaling ( G4bool input)
inline

Definition at line 185 of file G4DNABornIonisationModel1.hh.

186{
187 spScaling = input;
188}

◆ SelectStationary()

void G4DNABornIonisationModel1::SelectStationary ( G4bool input)
inline

Definition at line 178 of file G4DNABornIonisationModel1.hh.

179{
180 statCode = input;
181}

◆ StartTracking()

void G4DNABornIonisationModel1::StartTracking ( G4Track * track)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 312 of file G4DNABornIonisationModel1.cc.

313{
314 fTrack = track;
315}

◆ TransferedEnergy()

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

Definition at line 1012 of file G4DNABornIonisationModel1.cc.

1016{
1017 G4double nrj = 0.;
1018
1019 G4double valueK1 = 0;
1020 G4double valueK2 = 0;
1021 G4double valuePROB21 = 0;
1022 G4double valuePROB22 = 0;
1023 G4double valuePROB12 = 0;
1024 G4double valuePROB11 = 0;
1025
1026 G4double nrjTransf11 = 0;
1027 G4double nrjTransf12 = 0;
1028 G4double nrjTransf21 = 0;
1029 G4double nrjTransf22 = 0;
1030
1031 if (particleDefinition == G4Electron::ElectronDefinition())
1032 {
1033 // Protection against out of boundary access - electron case : 1 MeV
1034 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
1035 //
1036
1037 // k should be in eV
1038 auto k2 = std::upper_bound(eTdummyVec.begin(),
1039 eTdummyVec.end(),
1040 k);
1041 auto k1 = k2 - 1;
1042
1043 /*
1044 G4cout << "----> k=" << k
1045 << " " << *k1
1046 << " " << *k2
1047 << " " << random
1048 << " " << ionizationLevelIndex
1049 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1050 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
1051 << G4endl;
1052 */
1053
1054 // SI : the following condition avoids situations where random >last vector element
1055 if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1056 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
1057 {
1058 auto prob12 =
1059 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1060 eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1061 random);
1062
1063 auto prob11 = prob12 - 1;
1064
1065 auto prob22 =
1066 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1067 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1068 random);
1069
1070 auto prob21 = prob22 - 1;
1071
1072 valueK1 = *k1;
1073 valueK2 = *k2;
1074 valuePROB21 = *prob21;
1075 valuePROB22 = *prob22;
1076 valuePROB12 = *prob12;
1077 valuePROB11 = *prob11;
1078
1079 /*
1080 G4cout << " " << random << " " << valuePROB11 << " "
1081 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1082 */
1083
1084 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1085 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1086 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1087 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1088
1089 /*
1090 G4cout << " " << ionizationLevelIndex << " "
1091 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1092
1093 G4cout << " " << random << " " << nrjTransf11 << " "
1094 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1095 */
1096
1097 }
1098
1099 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1100 if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
1101 {
1102 auto prob22 =
1103 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1104 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1105 random);
1106
1107 auto prob21 = prob22 - 1;
1108
1109 valueK1 = *k1;
1110 valueK2 = *k2;
1111 valuePROB21 = *prob21;
1112 valuePROB22 = *prob22;
1113
1114 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1115
1116 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1117 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1118
1119 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1120 valuePROB22,
1121 random,
1122 nrjTransf21,
1123 nrjTransf22);
1124
1125 // zeros are explicitly set
1126
1127 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1128
1129 /*
1130 G4cout << " " << ionizationLevelIndex << " "
1131 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1132
1133 G4cout << " " << random << " " << nrjTransf11 << " "
1134 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1135
1136 G4cout << "ici" << " " << value << G4endl;
1137 */
1138
1139 return value;
1140 }
1141 }
1142 //
1143 else if (particleDefinition == G4Proton::ProtonDefinition())
1144 {
1145 // Protection against out of boundary access - proton case : 100 MeV
1146 if (k==pTdummyVec.back()) k=k*(1.-1e-12);
1147 //
1148
1149 // k should be in eV
1150
1151 auto k2 = std::upper_bound(pTdummyVec.begin(),
1152 pTdummyVec.end(),
1153 k);
1154
1155 auto k1 = k2 - 1;
1156
1157 /*
1158 G4cout << "----> k=" << k
1159 << " " << *k1
1160 << " " << *k2
1161 << " " << random
1162 << " " << ionizationLevelIndex
1163 << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1164 << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1165 << G4endl;
1166 */
1167
1168 // SI : the following condition avoids situations where random > last vector element,
1169 // for eg. when the last element is zero
1170 if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1171 && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1172 {
1173 auto prob12 =
1174 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1175 pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1176 random);
1177
1178 auto prob11 = prob12 - 1;
1179
1180 auto prob22 =
1181 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1182 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1183 random);
1184
1185 auto prob21 = prob22 - 1;
1186
1187 valueK1 = *k1;
1188 valueK2 = *k2;
1189 valuePROB21 = *prob21;
1190 valuePROB22 = *prob22;
1191 valuePROB12 = *prob12;
1192 valuePROB11 = *prob11;
1193
1194 /*
1195 G4cout << " " << random << " " << valuePROB11 << " "
1196 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1197 */
1198
1199 nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1200 nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1201 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1202 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1203
1204 /*
1205 G4cout << " " << ionizationLevelIndex << " "
1206 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1207
1208 G4cout << " " << random << " " << nrjTransf11 << " "
1209 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1210 */
1211 }
1212
1213 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1214
1215 if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1216 {
1217 auto prob22 =
1218 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1219 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1220 random);
1221
1222 auto prob21 = prob22 - 1;
1223
1224 valueK1 = *k1;
1225 valueK2 = *k2;
1226 valuePROB21 = *prob21;
1227 valuePROB22 = *prob22;
1228
1229 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1230
1231 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1232 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1233
1234 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1235 valuePROB22,
1236 random,
1237 nrjTransf21,
1238 nrjTransf22);
1239
1240 // zeros are explicitly set
1241
1242 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1243
1244 /*
1245 G4cout << " " << ionizationLevelIndex << " "
1246 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1247
1248 G4cout << " " << random << " " << nrjTransf11 << " "
1249 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1250
1251 G4cout << "ici" << " " << value << G4endl;
1252 */
1253
1254 return value;
1255 }
1256 }
1257 // End electron and proton cases
1258
1259 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1260 * nrjTransf22;
1261 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1262
1263 if (nrjTransfProduct != 0.)
1264 {
1265 nrj = QuadInterpolator(valuePROB11,
1266 valuePROB12,
1267 valuePROB21,
1268 valuePROB22,
1269 nrjTransf11,
1270 nrjTransf12,
1271 nrjTransf21,
1272 nrjTransf22,
1273 valueK1,
1274 valueK2,
1275 k,
1276 random);
1277 }
1278 //G4cout << nrj << endl;
1279
1280 return nrj;
1281}

Member Data Documentation

◆ fParticleChangeForGamma

G4ParticleChangeForGamma* G4DNABornIonisationModel1::fParticleChangeForGamma
protected

Definition at line 94 of file G4DNABornIonisationModel1.hh.

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

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