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