| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Created on 2016/05/02
27 //
28 // Authors: D Sakata, S. Incerti
29 //
30 // This class perform electric excitation for electron transportation in gold,
31 // based on Dirac B-Spline R-Matrix method with scaled experimental data
32 // for low energy.
33 // See following reference paper
34 // Phys.Rev.A77,062711(2008) and Phys.Rev.A78,042713(2008)
35
36 #include "G4DNADiracRMatrixExcitationModel.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4UAtomicDeexcitation.hh"
39 #include "G4LossTableManager.hh"
40 #include "G4Gamma.hh"
41 #include "G4RandomDirection.hh"
42
43 #include <vector>
44
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46
47 G4DNADiracRMatrixExcitationModel::G4DNADiracRMatrixExcitationModel
48 (const G4ParticleDefinition*,const G4String& nam) :
49 G4VEmModel(nam)
50 {
51 fpMaterialDensity = nullptr;
52 fHighEnergyLimit = 0;
53 fExperimentalEnergyLimit= 0;
54 fLowEnergyLimit = 0;
55 fParticleDefinition = nullptr;
56
57 verboseLevel = 0;
58
59 if (verboseLevel > 0)
60 {
61 G4cout << "Dirac R-matrix excitation model is constructed " << G4endl;
62 }
63
64 fParticleChangeForGamma = nullptr;
65 statCode = false;
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
70 G4DNADiracRMatrixExcitationModel::~G4DNADiracRMatrixExcitationModel()
71 {
72 delete fTableData;
73 }
74
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77 void G4DNADiracRMatrixExcitationModel::Initialise
78 (const G4ParticleDefinition* particle,const G4DataVector& /*cuts*/)
79 {
80
81 if (verboseLevel > 3)
82 {
83 G4cout <<
84 "Calling G4DNADiracRMatrixExcitationModel::Initialise()"
85 << G4endl;
86 }
87
88 fParticleDefinition = particle;
89
90 if(particle->GetParticleName() == "e-")
91 {
92 fTableFile = "dna/sigma_excitation_e_diracrmatrix_Z79";
93 fLowEnergyLimit = 10 * eV;
94 fExperimentalEnergyLimit = 577.* eV;
95 fHighEnergyLimit = 1.0 * GeV;
96 }
97 else
98 {
99 G4Exception("G4DNADiracRMatrixExcitationModel::Initialise","em0001",
100 FatalException,"Not defined for other particles than electrons.");
101 return;
102 }
103
104 G4double scaleFactor = 1. * cm * cm;
105 fTableData = new G4DNACrossSectionDataSet
106 (new G4LogLogInterpolation,eV,scaleFactor );
107 fTableData->LoadData(fTableFile);
108
109 if( verboseLevel>0 )
110 {
111 G4cout << "Dirac R-matrix excitation model is initialized " << G4endl
112 << "Energy range: "
113 << LowEnergyLimit() / eV << " eV - "<< HighEnergyLimit() / keV << " keV "
114 << " for "<< particle->GetParticleName()
115 << G4endl;
116 }
117
118 if (isInitialised){return;}
119 fParticleChangeForGamma = GetParticleChangeForGamma();
120 isInitialised = true;
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124
125 G4double G4DNADiracRMatrixExcitationModel::CrossSectionPerVolume
126 (const G4Material* material,
127 const G4ParticleDefinition* particleDefinition,
128 G4double ekin,
129 G4double,
130 G4double)
131 {
132 if (verboseLevel > 3)
133 {
134 G4cout <<
135 "Calling CrossSectionPerVolume() of G4DNADiracRMatrixExcitationModel"
136 << G4endl;
137 }
138
139 G4double atomicNDensity = material->GetAtomicNumDensityVector()[0];
140
141 // Protection: for single element
142 if(material->GetNumberOfElements()>1) return 0.;
143
144 G4double z = material->GetZ();
145
146 // Protection: for Gold
147 if(z!=79){return 0.;}
148
149 G4double sigma=0.;
150
151 if(atomicNDensity!= 0.0)
152 {
153 if (ekin >= fLowEnergyLimit && ekin < fExperimentalEnergyLimit)
154 {
155 sigma = fTableData->FindValue(ekin);
156 }
157 else if ((fExperimentalEnergyLimit <= ekin) && (ekin < fHighEnergyLimit))
158 {
159 sigma = GetExtendedTotalCrossSection(material,particleDefinition,ekin);
160 }
161
162 if (verboseLevel > 2)
163 {
164 G4cout<<"__________________________________" << G4endl;
165 G4cout<<"=== G4DNADiracRMatrixExcitationModel - XS INFO START"<<G4endl;
166 G4cout<<"=== Kinetic energy (eV)=" << ekin/eV << " particle : "
167 <<particleDefinition->GetParticleName() << G4endl;
168 G4cout<<"=== Cross section per atom for Z="<<z<<" is (cm^2)"
169 <<sigma/cm/cm << G4endl;
170 G4cout<<"=== Cross section per atom for Z="<<z<<" is (cm^-1)="
171 <<sigma*atomicNDensity/(1./cm) << G4endl;
172 G4cout<<"=== G4DNADiracRMatrixExcitationModel - XS INFO END"<<G4endl;
173 }
174 }
175
176 return sigma*atomicNDensity;
177 }
178
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180
181 void G4DNADiracRMatrixExcitationModel::SampleSecondaries
182 (std::vector<G4DynamicParticle*>* /*fvect*/,
183 const G4MaterialCutsCouple* couple,
184 const G4DynamicParticle* aDynamicParticle,
185 G4double,G4double)
186 {
187
188 if (verboseLevel > 3)
189 {
190 G4cout <<
191 "Calling SampleSecondaries() of G4DNADiracRMatrixExcitationModel"
192 << G4endl;
193 }
194
195 G4ParticleDefinition* particle = aDynamicParticle->GetDefinition();
196 G4double k = aDynamicParticle->GetKineticEnergy();
197
198 G4int level = RandomSelect(couple->GetMaterial(),particle,
199 k);
200 G4double excitationEnergy = ExcitationEnergyAu[level]*eV;
201 G4double newEnergy = k - excitationEnergy;
202
203 if (newEnergy > 0)
204 {
205 //Energy Loss
206 fParticleChangeForGamma->ProposeMomentumDirection
207 (aDynamicParticle->GetMomentumDirection());
208 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
209 if(!statCode) fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
210 else fParticleChangeForGamma->SetProposedKineticEnergy(k);
211 }
212 }
213
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
215
216 G4double G4DNADiracRMatrixExcitationModel::GetExtendedTotalCrossSection
217 (const G4Material* material,
218 const G4ParticleDefinition* particle,
219 G4double kineticEnergy)
220 {
221 G4double value=0;
222
223 size_t N=fTableData->NumberOfComponents();
224
225 for(int i=0;i<(int)N;i++){
226 value = value+GetExtendedPartialCrossSection(material,i,particle,
227 kineticEnergy);
228 }
229
230 return value;
231 }
232
233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
234
235 G4double G4DNADiracRMatrixExcitationModel::GetExtendedPartialCrossSection
236 (const G4Material*,
237 G4int level,
238 const G4ParticleDefinition* particle,
239 G4double kineticEnergy)
240 {
241 G4double value=0;
242
243 if(particle->GetParticleName()=="e-"){
244
245 if(level==0){
246 // y = [0]+[1]/pow(x-2,2)
247 value = paramFuncTCS_5dto6s1[0]+paramFuncTCS_5dto6s1[1]
248 /std::pow(kineticEnergy/eV-paramFuncTCS_5dto6s1[2],2);
249 }
250 else if(level==1){
251 // y = [0]+[1]/pow(x-2,2)
252 value = paramFuncTCS_5dto6s2[0]+paramFuncTCS_5dto6s2[1]
253 /std::pow(kineticEnergy/eV-paramFuncTCS_5dto6s2[2],2);
254 }
255 else if(level==2){
256 // y = [0]+[1]*log(x-2)/(x-[2])
257 value = paramFuncTCS_6sto6p1[0]+paramFuncTCS_6sto6p1[1]
258 *G4Log(kineticEnergy/eV-paramFuncTCS_6sto6p1[2])
259 /(kineticEnergy/eV-paramFuncTCS_6sto6p1[2]);
260 }
261 else if(level==3){
262 // y = [0]+[1]*log(x-2)/(x-[2])
263 value = paramFuncTCS_6sto6p2[0]+paramFuncTCS_6sto6p2[1]
264 *G4Log(kineticEnergy/eV-paramFuncTCS_6sto6p2[2])
265 /(kineticEnergy/eV-paramFuncTCS_6sto6p2[2]);
266 }
267 }
268
269 return value*cm*cm;
270 }
271
272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
273
274 G4int G4DNADiracRMatrixExcitationModel::RandomSelect
275 (const G4Material* material,
276 const G4ParticleDefinition* particle,
277 G4double kineticEnergy)
278 {
279 G4double value = 0.;
280
281 std::size_t NOfComp = fTableData->NumberOfComponents();
282
283 std::vector<G4double> valuesBuffer(NOfComp, 0.0);
284
285 const auto n = (G4int)fTableData->NumberOfComponents();
286
287 G4int i(n);
288
289 while (i > 0)
290 {
291 --i;
292 if
293 ((fLowEnergyLimit<=kineticEnergy)&&(kineticEnergy<fExperimentalEnergyLimit))
294 {
295 valuesBuffer[i] = fTableData->GetComponent(i)->FindValue(kineticEnergy);
296 }
297 else if
298 ((fExperimentalEnergyLimit<=kineticEnergy)&&(kineticEnergy<fHighEnergyLimit))
299 {
300 valuesBuffer[i]
301 = GetExtendedPartialCrossSection(material,i,particle,kineticEnergy);
302 }
303 value += valuesBuffer[i];
304 }
305 value *= G4UniformRand();
306 i = n;
307 while (i > 0)
308 {
309 --i;
310 if (valuesBuffer[i] > value)
311 {
312 return i;
313 }
314 value -= valuesBuffer[i];
315 }
316 return 9999;
317 }
318
319