| 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 //
27 // -------------------------------------------
28 //
29 // Author: E.Mendoza
30 //
31 // Creation date: May 2024
32 //
33 // Modifications:
34 //
35 // -------------------------------------------
36 //
37 // NuDEX code (https://doi.org/10.1016/j.nima
38 //
39
40
41
42 #include "G4NuDEXInternalConversion.hh"
43
44
45
46 //If alpha>0, use that value
47 G4bool G4NuDEXInternalConversion::SampleIntern
48
49 if(theZ<MINZINTABLES){ //then we have no inf
50 if(alpha<0){
51 Ne=0;
52 Ng=0;
53 return false;
54 }
55 else{
56 G4double rand=theRandom4->Uniform(0,alph
57 if(rand<alpha){ //then electron conversi
58 Ne=1;
59 Ng=0;
60 Eele[0]=Ene; //which is not correct, but we
61 return true;
62 }
63 return false;
64 }
65 }
66
67
68 Ne=0;
69 Ng=0;
70
71 if(multipolarity==0){ //maybe it is better t
72 //return true;
73 if(alpha<=0){
74 return false;
75 }
76 //NuDEXException(__FILE__,std::to_string(_
77 }
78
79 G4bool usegivenalpha=true;
80 if(NShells==0 || std::abs(multipolarity)>ICC
81 if(alpha<0){
82 usegivenalpha=false;
83 alpha=GetICC(Ene,multipolarity);
84 }
85
86 G4double rand=theRandom4->Uniform(0,alpha+1)
87 if(rand<alpha){ //then electron conversion
88 if(!CalculateProducts){return true;}
89 //Select the orbital:
90 if(usegivenalpha){rand=rand*GetICC(Ene,mul
91 G4double cumul=0;
92 for(G4int i=1;i<NShells;i++){
93 cumul+=GetICC(Ene,multipolarity,i);
94 //std::cout<<Ene<<" "<<multipolarity<<"
95 if(cumul>=rand || multipolarity==0){ //t
96 Ne=1;
97 Eele[0]=Ene-BindingEnergy[i];
98 FillElectronHole(i); //now there is a hole t
99 if(Eele[0]<0){
100 std::cout<<" For Z = "<<theZ<<" and orbita
101 std::cout<<" Given alpha is "<<alpha<<" ("
102 for(G4int j=1;j<=NShells;j++){
103 std::cout<<j<<" "<<GetICC(Ene,multipola
104 }
105 Eele[0]=0;
106 }
107 return true;
108 }
109 }
110 std::cout<<" ############ Warning in "<<__
111 std::cout<<" Given alpha is "<<alpha<<" an
112 for(G4int i=1;i<=NShells;i++){
113 std::cout<<i<<" "<<GetICC(Ene,multipola
114 }
115 Ne=1;
116 Eele[0]=Ene-BindingEnergy[NShells-1];
117 return true;
118 }
119
120 return false;
121 }
122
123
124 void G4NuDEXInternalConversion::FillElectronHo
125
126 //A very simplified version of the process (
127
128 G4double fluoyield=0;
129 if(i_shell==1){ //K-shell
130 //Hubbell et al. (1994) formula for the fl
131 G4double C0=0.0370,C1=0.03112,C2=5.44e-5,C
132 G4double w_fac=std::pow(C0+C1*theZ+C2*theZ
133 fluoyield=w_fac/(1.+w_fac);
134 }
135 else if(i_shell>=2 && i_shell<=4){ //L-shell
136 //Hubbell et al. (1994) formula for the fl
137 if(theZ>=3 && theZ<=36){
138 fluoyield=1.939e-8*std::pow(theZ,3.8874)
139 }
140 else if(theZ>36){
141 G4double C0=0.17765,C1=0.00298937,C2=8.9
142 G4double w_fac=std::pow(C0+C1*theZ+C2*th
143 fluoyield=w_fac/(1.+w_fac);
144 }
145 }
146
147
148 G4double rand=theRandom4->Uniform(0,1);
149 if(rand<fluoyield){ //gamma emission
150 Egam[Ng]=BindingEnergy[i_shell];
151 Ng++;
152 }
153 else{ //electron emission
154 Eele[Ne]=BindingEnergy[i_shell];
155 Ne++;
156 }
157
158
159 }
160
161
162
163
164
165
166
167 //If i_shell<0 --> the total alpha
168 G4double G4NuDEXInternalConversion::GetICC(G4d
169
170 if(theZ<MINZINTABLES){ //then we have no inf
171 return 0;
172 }
173
174 if(NShells==0 || std::abs(multipolarity)>ICC
175
176 //-----------------------------------------
177 //Total:
178 //The following line does not work, due to i
179 //if(i_shell<0){i_shell=NShells;}
180
181 if(i_shell<0){
182 G4double result=0;
183 for(G4int i=1;i<NShells;i++){
184 result+=GetICC(Ene,multipolarity,i);
185 }
186 return result;
187 }
188 //-----------------------------------------
189
190
191 if(Ene<BindingEnergy[i_shell]){return 0;}
192
193 if(np[i_shell]==0){
194 std::cout<<" shell "<<i_shell<<" has not b
195 NuDEXException(__FILE__,std::to_string(__L
196 }
197
198 if(i_shell==NShells && Ene<Eg[i_shell][0]){
199 G4double total=0;
200 for(G4int i=1;i<NShells;i++){total+=GetICC
201 return total;
202 }
203
204 if(multipolarity>0){
205 return Interpolate(Ene,np[i_shell],Eg[i_sh
206 }
207 else if(multipolarity<0){
208 return Interpolate(Ene,np[i_shell],Eg[i_sh
209 }
210 return 0;
211 }
212
213
214 G4NuDEXInternalConversion::G4NuDEXInternalConv
215 theZ=Z;
216 NShells=0;
217 Ne=Ng=0;
218 for(G4int i=0;i<ICC_MAXNSHELLS;i++){
219 Eg[i]=0; np[i]=0; BindingEnergy[i]=0;
220 for(G4int j=0;j<ICC_NMULTIP;j++){
221 Icc_E[j][i]=0; Icc_M[j][i]=0;
222 }
223 }
224 theRandom4= new G4NuDEXRandom(1234567);
225 }
226
227
228 G4NuDEXInternalConversion::~G4NuDEXInternalCon
229 for(G4int i=0;i<ICC_MAXNSHELLS;i++){
230 if(Eg[i]!=0){delete [] Eg[i];}
231 for(G4int j=0;j<ICC_NMULTIP;j++){
232 if(Icc_E[j][i]!=0){delete [] Icc_E[j][i]
233 if(Icc_M[j][i]!=0){delete [] Icc_M[j][i]
234 }
235 }
236 delete theRandom4;
237 }
238
239 void G4NuDEXInternalConversion::PrintICC(std::
240
241 char word[1000];
242 out<<" #####################################
243 out<<" ICC"<<std::endl;
244 out<<" Z = "<<theZ<<std::endl;
245 out<<" NShells = "<<NShells<<std::endl;
246 out<<" -------------------------------------
247 out<<" Total calculated from the sum of the
248 out<<" E_g E1 E2
249 for(G4int j=0;j<np[NShells];j++){
250 snprintf(word,1000,"%10.4g",Eg[NShells][j]
251 for(G4int k=0;k<ICC_NMULTIP;k++){
252 snprintf(word,1000," %10.4g",Icc_E[k][N
253 }
254 for(G4int k=0;k<ICC_NMULTIP;k++){
255 snprintf(word,1000," %10.4g",Icc_M[k][N
256 }
257 out<<std::endl;
258 }
259 out<<" -------------------------------------
260 for(G4int i=0;i<NShells;i++){
261 out<<" -----------------------------------
262 out<<" Binding energy = "<<BindingEnergy[i
263 out<<" E_g E1 E2
264 for(G4int j=0;j<np[i];j++){
265 snprintf(word,1000,"%10.4g",Eg[i][j]); o
266 for(G4int k=0;k<ICC_NMULTIP;k++){
267 snprintf(word,1000," %10.4g",Icc_E[k][i][j]
268 }
269 for(G4int k=0;k<ICC_NMULTIP;k++){
270 snprintf(word,1000," %10.4g",Icc_M[k][i][j]
271 }
272 out<<std::endl;
273 }
274 out<<" -----------------------------------
275 }
276 out<<" #####################################
277 }
278
279 void G4NuDEXInternalConversion::Init(const cha
280
281 if(theZ<MINZINTABLES){ //then we have no inf
282 return;
283 }
284
285 if(NShells!=0){ //Init only once
286 NuDEXException(__FILE__,std::to_string(__L
287 }
288
289 std::ifstream in(fname);
290 if(!in.good()){
291 std::cout<<" ################ Error openin
292 NuDEXException(__FILE__,std::to_string(__L
293 }
294 std::string word;
295 NShells=1;
296 while(in>>word){
297 if(word.c_str()[0]=='Z' && word.c_str()[1]
298 if(std::atoi(&(word.c_str()[2]))==theZ){
299 in>>word>>word;
300 G4int orbindex;
301 if(word==std::string("Total")){
302 in.ignore(1000,'\n');
303 in.ignore(1000,'\n');
304 orbindex=0;
305 }
306 else{
307 orbindex=NShells;
308 in>>word>>word>>BindingEnergy[NShells]>>wo
309 BindingEnergy[NShells]*=1.e-6; // all in M
310 in.ignore(1000,'\n');
311 in.ignore(1000,'\n');
312 }
313 //------------------------------------------
314 size_t sz,sz2;
315 G4int np_tmp=0;
316 G4double Eg_tmp[1000],Icc_E_tmp[ICC_NMULTIP]
317 while(getline(in,word)){
318 if(word.size()<100){
319 np[orbindex]=np_tmp;
320 Eg[orbindex]=new G4double[np_tmp];
321 for(G4int j=0;j<np_tmp;j++){
322 Eg[orbindex][j]=Eg_tmp[j];
323 }
324 for(G4int i=0;i<ICC_NMULTIP;i++){
325 Icc_E[i][orbindex]=new G4double[np_tmp
326 Icc_M[i][orbindex]=new G4double[np_tmp
327 }
328 for(G4int i=0;i<ICC_NMULTIP;i++){
329 for(G4int j=0;j<np_tmp;j++){
330 Icc_E[i][orbindex][j]=Icc_E_tmp[i][j];
331 Icc_M[i][orbindex][j]=Icc_M_tmp[i][j];
332 }
333 }
334 if(orbindex!=0){NShells++;}
335 break;
336 }
337 else{
338 sz=0;
339 Eg_tmp[np_tmp]=std::stof(word,&sz2);
340 Eg_tmp[np_tmp]*=1.e-3; //all to MeV
341 sz+=sz2;
342 for(G4int i=0;i<ICC_NMULTIP;i++){
343 Icc_E_tmp[i][np_tmp]=std::stof(word.su
344 }
345 for(G4int i=0;i<ICC_NMULTIP;i++){
346 Icc_M_tmp[i][np_tmp]=std::stof(word.su
347 }
348 if((G4int)(std::stof(word.substr(sz),&sz
349 NuDEXException(__FILE__,std::to_string
350 }
351 sz+=sz2;
352 sz2=word.find_first_not_of(' ',sz);
353 if(np_tmp==0){OrbitalName[orbindex]=word
354 np_tmp++;
355 }
356 }
357 if(orbindex==0){break;}
358 //------------------------------------------
359 }
360 }
361 }
362 if(!in.good()){
363 NuDEXException(__FILE__,std::to_string(__L
364 }
365 in.close();
366
367 MakeTotal();
368 }
369
370
371
372 // Total Icc goes to nShell=NShells
373 void G4NuDEXInternalConversion::MakeTotal(){
374
375 if(np[0]==0 || Eg[0]==0){
376 NuDEXException(__FILE__,std::to_string(__L
377 }
378
379 //We evaluate it in the same frame as the to
380 BindingEnergy[NShells]=0;
381 np[NShells]=np[0];
382 Eg[NShells]=new G4double[np[NShells]];
383 for(G4int i=0;i<ICC_NMULTIP;i++){
384 Icc_E[i][NShells]=new G4double[np[NShells]
385 Icc_M[i][NShells]=new G4double[np[NShells]
386 }
387 for(G4int k=0;k<np[NShells];k++){
388 for(G4int j=0;j<ICC_NMULTIP;j++){
389 Icc_E[j][NShells][k]=0;
390 Icc_M[j][NShells][k]=0;
391 }
392 }
393
394
395 for(G4int k=0;k<np[NShells];k++){
396 Eg[NShells][k]=Eg[0][k];
397 for(G4int i=1;i<NShells;i++){
398 for(G4int j=0;j<ICC_NMULTIP;j++){
399 Icc_E[j][NShells][k]+=GetICC(Eg[NShells][k],
400 Icc_M[j][NShells][k]+=GetICC(Eg[NShells][k],
401 }
402 }
403 }
404
405 }
406
407 //if val>x[npmax] then ---> return 0
408 G4double G4NuDEXInternalConversion::Interpolat
409
410 G4int i_interp=-1;
411 for(G4int i=1;i<npoints;i++){
412 if(x[i]>=val){i_interp=i-1; break;}
413 }
414 if(i_interp<0){return 0;}
415
416
417 /*
418 //y=a0+a1*x
419 G4double a1=(y[i_interp+1]-y[i_interp])/(x[i
420 G4double a0=y[i_interp]-a1*x[i_interp];
421
422 return (a0+a1*val);
423 */
424
425 //It is better a log-log interpolation:
426 if(y[i_interp+1]<=0 || y[i_interp]<=0 || x[i
427 //y=a0+a1*x
428 G4double a1=(y[i_interp+1]-y[i_interp])/(x
429 G4double a0=y[i_interp]-a1*x[i_interp];
430
431 return (a0+a1*val);
432 }
433
434 //log(y)=a0+a1*log(x)
435 G4double a1=std::log(y[i_interp+1]/y[i_inter
436 G4double a0=std::log(y[i_interp])-a1*std::lo
437
438 G4double result=std::exp(a0+a1*std::log(val)
439 return result;
440 }
441