| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Created on 2016/05/02
27 //
28 // Authors: D Sakata, S. Incerti
29 //
30 // This class perform electric excitation for
31 // based on Dirac B-Spline R-Matrix method wit
32 // for low energy.
33 // See following reference paper
34 // Phys.Rev.A77,062711(2008) and Phys.Rev.A78,
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........oo
46
47 G4DNADiracRMatrixExcitationModel::G4DNADiracRM
48 (const G4ParticleDefinition*,const G4String& n
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
62 }
63
64 fParticleChangeForGamma = nullptr;
65 statCode = false;
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oo
69
70 G4DNADiracRMatrixExcitationModel::~G4DNADiracR
71 {
72 delete fTableData;
73 }
74
75 //....oooOO0OOooo........oooOO0OOooo........oo
76
77 void G4DNADiracRMatrixExcitationModel::Initial
78 (const G4ParticleDefinition* particle,const G4
79 {
80
81 if (verboseLevel > 3)
82 {
83 G4cout <<
84 "Calling G4DNADiracRMatrixExcitatio
85 << G4endl;
86 }
87
88 fParticleDefinition = particle;
89
90 if(particle->GetParticleName() == "e-")
91 {
92 fTableFile = "dna/sigma_excitation_e_dirac
93 fLowEnergyLimit = 10 * eV;
94 fExperimentalEnergyLimit = 577.* eV;
95 fHighEnergyLimit = 1.0 * GeV;
96 }
97 else
98 {
99 G4Exception("G4DNADiracRMatrixExcitationMod
100 FatalException,"Not defined for other
101 return;
102 }
103
104 G4double scaleFactor = 1. * cm * cm;
105 fTableData = new G4DNACrossSectionDataSet
106 (new G4LogLogInterpolation,eV,sc
107 fTableData->LoadData(fTableFile);
108
109 if( verboseLevel>0 )
110 {
111 G4cout << "Dirac R-matrix excitation model
112 << "Energy range: "
113 << LowEnergyLimit() / eV << " eV - "<< Hig
114 << " for "<< particle->GetParticleName()
115 << G4endl;
116 }
117
118 if (isInitialised){return;}
119 fParticleChangeForGamma = GetParticleChangeF
120 isInitialised = true;
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oo
124
125 G4double G4DNADiracRMatrixExcitationModel::Cro
126 (const G4Material* ma
127 const G4ParticleDefi
128 G4double ekin,
129 G4double,
130 G4double)
131 {
132 if (verboseLevel > 3)
133 {
134 G4cout <<
135 "Calling CrossSectionPerVolume() of G4
136 << G4endl;
137 }
138
139 G4double atomicNDensity = material->GetAtomi
140
141 // Protection: for single element
142 if(material->GetNumberOfElements()>1) return
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 < fExp
154 {
155 sigma = fTableData->FindValue(ekin);
156 }
157 else if ((fExperimentalEnergyLimit <= ekin
158 {
159 sigma = GetExtendedTotalCrossSection(mat
160 }
161
162 if (verboseLevel > 2)
163 {
164 G4cout<<"_______________________________
165 G4cout<<"=== G4DNADiracRMatrixExcitation
166 G4cout<<"=== Kinetic energy (eV)=" << ek
167 <<particleDefinition->GetParticleN
168 G4cout<<"=== Cross section per atom for
169 <<sigma/cm/cm << G4endl;
170 G4cout<<"=== Cross section per atom for
171 <<sigma*atomicNDensity/(1./cm) <<
172 G4cout<<"=== G4DNADiracRMatrixExcitation
173 }
174 }
175
176 return sigma*atomicNDensity;
177 }
178
179 //....oooOO0OOooo........oooOO0OOooo........oo
180
181 void G4DNADiracRMatrixExcitationModel::SampleS
182 (std::vector<G4Dynami
183 const G4MaterialCuts
184 const G4DynamicParti
185 G4double,G4double)
186 {
187
188 if (verboseLevel > 3)
189 {
190 G4cout <<
191 "Calling SampleSecondaries() of G4DN
192 << G4endl;
193 }
194
195 G4ParticleDefinition* particle = aDynamicPar
196 G4double k = aDynamicPar
197
198 G4int level = RandomSelec
199
200 G4double excitationEnergy = ExcitationE
201 G4double newEnergy = k - excitat
202
203 if (newEnergy > 0)
204 {
205 //Energy Loss
206 fParticleChangeForGamma->ProposeMomentumDi
207 (aDynamicParticle->GetMoment
208 fParticleChangeForGamma->ProposeLocalEnerg
209 if(!statCode) fParticleChangeForGamma->Set
210 else fParticleChangeForGamma->Set
211 }
212 }
213
214 //....oooOO0OOooo........oooOO0OOooo........oo
215
216 G4double G4DNADiracRMatrixExcitationModel::Get
217 (const G4Material* m
218 const G4ParticleDef
219 G4double kineticEne
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+GetExtendedPartialCrossSectio
227
228 }
229
230 return value;
231 }
232
233 //....oooOO0OOooo........oooOO0OOooo........oo
234
235 G4double G4DNADiracRMatrixExcitationModel::Get
236 (const G4Material*,
237 G4int level,
238 const G4ParticleDefi
239 G4double kineticEner
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]+paramFun
248 /std::pow(kineticEnergy/eV-param
249 }
250 else if(level==1){
251 // y = [0]+[1]/pow(x-2,2)
252 value = paramFuncTCS_5dto6s2[0]+paramFun
253 /std::pow(kineticEnergy/eV-param
254 }
255 else if(level==2){
256 // y = [0]+[1]*log(x-2)/(x-[2])
257 value = paramFuncTCS_6sto6p1[0]+paramFun
258 *G4Log(kineticEnergy/eV-paramFun
259 /(kineticEnergy/eV-paramFuncTCS_
260 }
261 else if(level==3){
262 // y = [0]+[1]*log(x-2)/(x-[2])
263 value = paramFuncTCS_6sto6p2[0]+paramFun
264 *G4Log(kineticEnergy/eV-paramFun
265 /(kineticEnergy/eV-paramFuncTCS_
266 }
267 }
268
269 return value*cm*cm;
270 }
271
272 //....oooOO0OOooo........oooOO0OOooo........oo
273
274 G4int G4DNADiracRMatrixExcitationModel::Random
275 (const G4Material* ma
276 const G4ParticleDefi
277 G4double kineticEner
278 {
279 G4double value = 0.;
280
281 std::size_t NOfComp = fTableData->NumberOfCo
282
283 std::vector<G4double> valuesBuffer(NOfComp,
284
285 const auto n = (G4int)fTableData->NumberOfC
286
287 G4int i(n);
288
289 while (i > 0)
290 {
291 --i;
292 if
293 ((fLowEnergyLimit<=kineticEnergy)&&(kineti
294 {
295 valuesBuffer[i] = fTableData->GetCompone
296 }
297 else if
298 ((fExperimentalEnergyLimit<=kineticEnergy)
299 {
300 valuesBuffer[i]
301 = GetExtendedPartialCrossSection(m
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