Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Protected Attributes | List of all members
G4DNABornIonisationModel1 Class Reference

#include <G4DNABornIonisationModel1.hh>

+ Inheritance diagram for G4DNABornIonisationModel1:

Public Member Functions

 G4DNABornIonisationModel1 (const G4ParticleDefinition *p=0, const G4String &nam="DNABornIonisationModel")
 
virtual ~G4DNABornIonisationModel1 ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &= *(new G4DataVector()))
 
virtual G4double CrossSectionPerVolume (const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int, const G4ParticleDefinition *, G4double)
 
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 Initialise (const G4ParticleDefinition *, const G4DataVector &)=0
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
 
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 CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
 
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 G4MaterialCutsCouple *, const G4DynamicParticle *, 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 LPMFlag () 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 SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetFluctuationFlag (G4bool val)
 
void SetMasterThread (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)
 
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
 
size_t currentCoupleIndex = 0
 
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 G4DNABornIonisationModel1.hh.

Constructor & Destructor Documentation

◆ G4DNABornIonisationModel1()

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

Definition at line 45 of file G4DNABornIonisationModel1.cc.

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

◆ ~G4DNABornIonisationModel1()

G4DNABornIonisationModel1::~G4DNABornIonisationModel1 ( )
virtual

Definition at line 86 of file G4DNABornIonisationModel1.cc.

87{
88 // Cross section
89
90 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
91 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
92 {
93 G4DNACrossSectionDataSet* table = pos->second;
94 delete table;
95 }
96
97 // Final state
98
99 eVecm.clear();
100 pVecm.clear();
101}

Member Function Documentation

◆ CrossSectionPerVolume()

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

Reimplemented from G4VEmModel.

Definition at line 317 of file G4DNABornIonisationModel1.cc.

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

◆ DifferentialCrossSection()

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

Definition at line 744 of file G4DNABornIonisationModel1.cc.

748{
749 G4double sigma = 0.;
750
751 if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
752 {
753 G4double valueT1 = 0;
754 G4double valueT2 = 0;
755 G4double valueE21 = 0;
756 G4double valueE22 = 0;
757 G4double valueE12 = 0;
758 G4double valueE11 = 0;
759
760 G4double xs11 = 0;
761 G4double xs12 = 0;
762 G4double xs21 = 0;
763 G4double xs22 = 0;
764
765 if (particleDefinition == G4Electron::ElectronDefinition())
766 {
767
768 // Protection against out of boundary access - electron case : 1 MeV
769 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
770 //
771
772 // k should be in eV and energy transfer eV also
773
774 std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
775 eTdummyVec.end(),
776 k);
777
778 std::vector<G4double>::iterator t1 = t2 - 1;
779
780 // SI : the following condition avoids situations where energyTransfer >last vector element
781 if (energyTransfer <= eVecm[(*t1)].back()
782 && energyTransfer <= eVecm[(*t2)].back())
783 {
784 std::vector<G4double>::iterator e12 =
785 std::upper_bound(eVecm[(*t1)].begin(),
786 eVecm[(*t1)].end(),
787 energyTransfer);
788 std::vector<G4double>::iterator e11 = e12 - 1;
789
790 std::vector<G4double>::iterator e22 =
791 std::upper_bound(eVecm[(*t2)].begin(),
792 eVecm[(*t2)].end(),
793 energyTransfer);
794 std::vector<G4double>::iterator e21 = e22 - 1;
795
796 valueT1 = *t1;
797 valueT2 = *t2;
798 valueE21 = *e21;
799 valueE22 = *e22;
800 valueE12 = *e12;
801 valueE11 = *e11;
802
803 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
804 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
805 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
806 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
807
808 }
809
810 }
811
812 if (particleDefinition == G4Proton::ProtonDefinition())
813 {
814 // Protection against out of boundary access - proton case : 100 MeV
815 if (k==pTdummyVec.back()) k=k*(1.-1e-12);
816 //
817
818 // k should be in eV and energy transfer eV also
819 std::vector<G4double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),
820 pTdummyVec.end(),
821 k);
822 std::vector<G4double>::iterator t1 = t2 - 1;
823
824 std::vector<G4double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),
825 pVecm[(*t1)].end(),
826 energyTransfer);
827 std::vector<G4double>::iterator e11 = e12 - 1;
828
829 std::vector<G4double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),
830 pVecm[(*t2)].end(),
831 energyTransfer);
832 std::vector<G4double>::iterator e21 = e22 - 1;
833
834 valueT1 = *t1;
835 valueT2 = *t2;
836 valueE21 = *e21;
837 valueE22 = *e22;
838 valueE12 = *e12;
839 valueE11 = *e11;
840
841 xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
842 xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
843 xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
844 xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
845
846 }
847
848 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
849 if (xsProduct != 0.)
850 {
851 sigma = QuadInterpolator(valueE11,
852 valueE12,
853 valueE21,
854 valueE22,
855 xs11,
856 xs12,
857 xs21,
858 xs22,
859 valueT1,
860 valueT2,
861 k,
862 energyTransfer);
863 }
864 }
865
866 return sigma;
867}

◆ GetPartialCrossSection()

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

Reimplemented from G4VEmModel.

Definition at line 959 of file G4DNABornIonisationModel1.cc.

963{
964 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
965 pos = tableData.find(particle->GetParticleName());
966
967 if (pos != tableData.end())
968 {
969 G4DNACrossSectionDataSet* table = pos->second;
970 return table->GetComponent(level)->FindValue(kineticEnergy);
971 }
972
973 return 0;
974}
G4DNACrossSectionDataSet::GetComponent
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
Definition: G4DNACrossSectionDataSet.hh:68
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:109
G4VEMDataSet::FindValue
virtual G4double FindValue(G4double x, G4int componentId=0) const =0

◆ Initialise()

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

Implements G4VEmModel.

Definition at line 105 of file G4DNABornIonisationModel1.cc.

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

◆ SampleSecondaries()

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

Implements G4VEmModel.

Definition at line 407 of file G4DNABornIonisationModel1.cc.

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

◆ SelectFasterComputation()

void G4DNABornIonisationModel1::SelectFasterComputation ( G4bool  input)
inline

Definition at line 174 of file G4DNABornIonisationModel1.hh.

175{
176 fasterCode = input;
177}

◆ SelectSPScaling()

void G4DNABornIonisationModel1::SelectSPScaling ( G4bool  input)
inline

Definition at line 188 of file G4DNABornIonisationModel1.hh.

189{
190 spScaling = input;
191}

◆ SelectStationary()

void G4DNABornIonisationModel1::SelectStationary ( G4bool  input)
inline

Definition at line 181 of file G4DNABornIonisationModel1.hh.

182{
183 statCode = input;
184}

◆ TransferedEnergy()

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

Definition at line 1064 of file G4DNABornIonisationModel1.cc.

1068{
1069 G4double nrj = 0.;
1070
1071 G4double valueK1 = 0;
1072 G4double valueK2 = 0;
1073 G4double valuePROB21 = 0;
1074 G4double valuePROB22 = 0;
1075 G4double valuePROB12 = 0;
1076 G4double valuePROB11 = 0;
1077
1078 G4double nrjTransf11 = 0;
1079 G4double nrjTransf12 = 0;
1080 G4double nrjTransf21 = 0;
1081 G4double nrjTransf22 = 0;
1082
1083 if (particleDefinition == G4Electron::ElectronDefinition())
1084 {
1085 // Protection against out of boundary access - electron case : 1 MeV
1086 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
1087 //
1088
1089 // k should be in eV
1090 std::vector<G4double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),
1091 eTdummyVec.end(),
1092 k);
1093 std::vector<G4double>::iterator k1 = k2 - 1;
1094
1095 /*
1096 G4cout << "----> k=" << k
1097 << " " << *k1
1098 << " " << *k2
1099 << " " << random
1100 << " " << ionizationLevelIndex
1101 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1102 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
1103 << G4endl;
1104 */
1105
1106 // SI : the following condition avoids situations where random >last vector element
1107 if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1108 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
1109 {
1110 std::vector<G4double>::iterator prob12 =
1111 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1112 eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1113 random);
1114
1115 std::vector<G4double>::iterator prob11 = prob12 - 1;
1116
1117 std::vector<G4double>::iterator prob22 =
1118 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1119 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1120 random);
1121
1122 std::vector<G4double>::iterator prob21 = prob22 - 1;
1123
1124 valueK1 = *k1;
1125 valueK2 = *k2;
1126 valuePROB21 = *prob21;
1127 valuePROB22 = *prob22;
1128 valuePROB12 = *prob12;
1129 valuePROB11 = *prob11;
1130
1131 /*
1132 G4cout << " " << random << " " << valuePROB11 << " "
1133 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1134 */
1135
1136 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1137 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1138 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1139 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1140
1141 /*
1142 G4cout << " " << ionizationLevelIndex << " "
1143 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1144
1145 G4cout << " " << random << " " << nrjTransf11 << " "
1146 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1147 */
1148
1149 }
1150
1151 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1152 if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
1153 {
1154 std::vector<G4double>::iterator prob22 =
1155 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1156 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1157 random);
1158
1159 std::vector<G4double>::iterator prob21 = prob22 - 1;
1160
1161 valueK1 = *k1;
1162 valueK2 = *k2;
1163 valuePROB21 = *prob21;
1164 valuePROB22 = *prob22;
1165
1166 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1167
1168 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1169 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1170
1171 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1172 valuePROB22,
1173 random,
1174 nrjTransf21,
1175 nrjTransf22);
1176
1177 // zeros are explicitly set
1178
1179 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1180
1181 /*
1182 G4cout << " " << ionizationLevelIndex << " "
1183 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1184
1185 G4cout << " " << random << " " << nrjTransf11 << " "
1186 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1187
1188 G4cout << "ici" << " " << value << G4endl;
1189 */
1190
1191 return value;
1192 }
1193 }
1194 //
1195 else if (particleDefinition == G4Proton::ProtonDefinition())
1196 {
1197 // Protection against out of boundary access - proton case : 100 MeV
1198 if (k==pTdummyVec.back()) k=k*(1.-1e-12);
1199 //
1200
1201 // k should be in eV
1202
1203 std::vector<G4double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),
1204 pTdummyVec.end(),
1205 k);
1206
1207 std::vector<G4double>::iterator k1 = k2 - 1;
1208
1209 /*
1210 G4cout << "----> k=" << k
1211 << " " << *k1
1212 << " " << *k2
1213 << " " << random
1214 << " " << ionizationLevelIndex
1215 << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1216 << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1217 << G4endl;
1218 */
1219
1220 // SI : the following condition avoids situations where random > last vector element,
1221 // for eg. when the last element is zero
1222 if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1223 && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1224 {
1225 std::vector<G4double>::iterator prob12 =
1226 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1227 pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1228 random);
1229
1230 std::vector<G4double>::iterator prob11 = prob12 - 1;
1231
1232 std::vector<G4double>::iterator prob22 =
1233 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1234 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1235 random);
1236
1237 std::vector<G4double>::iterator prob21 = prob22 - 1;
1238
1239 valueK1 = *k1;
1240 valueK2 = *k2;
1241 valuePROB21 = *prob21;
1242 valuePROB22 = *prob22;
1243 valuePROB12 = *prob12;
1244 valuePROB11 = *prob11;
1245
1246 /*
1247 G4cout << " " << random << " " << valuePROB11 << " "
1248 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1249 */
1250
1251 nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1252 nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1253 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1254 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1255
1256 /*
1257 G4cout << " " << ionizationLevelIndex << " "
1258 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1259
1260 G4cout << " " << random << " " << nrjTransf11 << " "
1261 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1262 */
1263 }
1264
1265 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1266
1267 if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1268 {
1269 std::vector<G4double>::iterator prob22 =
1270 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1271 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1272 random);
1273
1274 std::vector<G4double>::iterator prob21 = prob22 - 1;
1275
1276 valueK1 = *k1;
1277 valueK2 = *k2;
1278 valuePROB21 = *prob21;
1279 valuePROB22 = *prob22;
1280
1281 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1282
1283 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1284 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1285
1286 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1287 valuePROB22,
1288 random,
1289 nrjTransf21,
1290 nrjTransf22);
1291
1292 // zeros are explicitly set
1293
1294 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1295
1296 /*
1297 G4cout << " " << ionizationLevelIndex << " "
1298 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1299
1300 G4cout << " " << random << " " << nrjTransf11 << " "
1301 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1302
1303 G4cout << "ici" << " " << value << G4endl;
1304 */
1305
1306 return value;
1307 }
1308 }
1309 // End electron and proton cases
1310
1311 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1312 * nrjTransf22;
1313 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1314
1315 if (nrjTransfProduct != 0.)
1316 {
1317 nrj = QuadInterpolator(valuePROB11,
1318 valuePROB12,
1319 valuePROB21,
1320 valuePROB22,
1321 nrjTransf11,
1322 nrjTransf12,
1323 nrjTransf21,
1324 nrjTransf22,
1325 valueK1,
1326 valueK2,
1327 k,
1328 random);
1329 }
1330 //G4cout << nrj << endl;
1331
1332 return nrj;
1333}

Member Data Documentation

◆ fParticleChangeForGamma

G4ParticleChangeForGamma* G4DNABornIonisationModel1::fParticleChangeForGamma
protected

Definition at line 89 of file G4DNABornIonisationModel1.hh.

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

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