| 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 //
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.2022.167894)
38 //
39
40
41
42 #include "G4NuDEXInternalConversion.hh"
43
44
45
46 //If alpha>0, use that value
47 G4bool G4NuDEXInternalConversion::SampleInternalConversion(G4double Ene,G4int multipolarity,G4double alpha,G4bool CalculateProducts){
48
49 if(theZ<MINZINTABLES){ //then we have no info
50 if(alpha<0){
51 Ne=0;
52 Ng=0;
53 return false;
54 }
55 else{
56 G4double rand=theRandom4->Uniform(0,alpha+1);
57 if(rand<alpha){ //then electron conversion
58 Ne=1;
59 Ng=0;
60 Eele[0]=Ene; //which is not correct, but we don't know the binding energy
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 to return true ... ?? --> no
72 //return true;
73 if(alpha<=0){
74 return false;
75 }
76 //NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
77 }
78
79 G4bool usegivenalpha=true;
80 if(NShells==0 || std::abs(multipolarity)>ICC_NMULTIP){return false;}
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,multipolarity)/alpha;} //renormalize rand to our alpha
91 G4double cumul=0;
92 for(G4int i=1;i<NShells;i++){
93 cumul+=GetICC(Ene,multipolarity,i);
94 //std::cout<<Ene<<" "<<multipolarity<<" "<<i<<" "<<GetICC(Ene,multipolarity,i)<<" "<<rand-1<<std::endl;
95 if(cumul>=rand || multipolarity==0){ //then is this orbital
96 Ne=1;
97 Eele[0]=Ene-BindingEnergy[i];
98 FillElectronHole(i); //now there is a hole there, in the filling procedure we emitt gammas and/or electrons
99 if(Eele[0]<0){
100 std::cout<<" For Z = "<<theZ<<" and orbital "<<OrbitalName[i]<<" --> Ene = "<<Ene<<" and BindingEnergy = "<<BindingEnergy[i]<<std::endl;
101 std::cout<<" Given alpha is "<<alpha<<" ("<<usegivenalpha<<"), rand = "<<rand<<" and tabulated alpha for Ene = "<<Ene<<" and mult = "<<multipolarity<<" is "<<GetICC(Ene,multipolarity)<<" -- cumul = "<<cumul<<std::endl;
102 for(G4int j=1;j<=NShells;j++){
103 std::cout<<j<<" "<<GetICC(Ene,multipolarity,j)<<std::endl;
104 }
105 Eele[0]=0;
106 }
107 return true;
108 }
109 }
110 std::cout<<" ############ Warning in "<<__FILE__<<", line "<<__LINE__<<" ############"<<std::endl;
111 std::cout<<" Given alpha is "<<alpha<<" and tabulated alpha for Ene = "<<Ene<<" and mult = "<<multipolarity<<" is "<<GetICC(Ene,multipolarity)<<" -- cumul = "<<cumul<<std::endl;
112 for(G4int i=1;i<=NShells;i++){
113 std::cout<<i<<" "<<GetICC(Ene,multipolarity,i)<<std::endl;
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::FillElectronHole(G4int i_shell){
125
126 //A very simplified version of the process (... and false). It can be done with accuracy with G4AtomicTransitionManager
127
128 G4double fluoyield=0;
129 if(i_shell==1){ //K-shell
130 //Hubbell et al. (1994) formula for the fluorescence yield:
131 G4double C0=0.0370,C1=0.03112,C2=5.44e-5,C3=-1.25e-6;
132 G4double w_fac=std::pow(C0+C1*theZ+C2*theZ*theZ+C3*theZ*theZ*theZ,4);
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 fluorescence yield:
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.91297e-5,C3=-2.67184e-7;
142 G4double w_fac=std::pow(C0+C1*theZ+C2*theZ*theZ+C3*theZ*theZ*theZ,4);
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(G4double Ene,G4int multipolarity,G4int i_shell){
169
170 if(theZ<MINZINTABLES){ //then we have no info
171 return 0;
172 }
173
174 if(NShells==0 || std::abs(multipolarity)>ICC_NMULTIP){return 0;}
175
176 //-----------------------------------------
177 //Total:
178 //The following line does not work, due to interpolation below binding energies:
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 been initialized"<<std::endl;
195 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
196 }
197
198 if(i_shell==NShells && Ene<Eg[i_shell][0]){ //then we cannot extrapolate, because of the binding energies of the different shells
199 G4double total=0;
200 for(G4int i=1;i<NShells;i++){total+=GetICC(Ene,multipolarity,i);}
201 return total;
202 }
203
204 if(multipolarity>0){
205 return Interpolate(Ene,np[i_shell],Eg[i_shell],Icc_E[multipolarity-1][i_shell]);
206 }
207 else if(multipolarity<0){
208 return Interpolate(Ene,np[i_shell],Eg[i_shell],Icc_M[(-multipolarity)-1][i_shell]);
209 }
210 return 0;
211 }
212
213
214 G4NuDEXInternalConversion::G4NuDEXInternalConversion(G4int Z){
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::~G4NuDEXInternalConversion(){
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::ostream &out){
240
241 char word[1000];
242 out<<" ######################################################################################################################################### "<<std::endl;
243 out<<" ICC"<<std::endl;
244 out<<" Z = "<<theZ<<std::endl;
245 out<<" NShells = "<<NShells<<std::endl;
246 out<<" ----------------------------------------------------------------------------------------------------------------------------------------"<<std::endl;
247 out<<" Total calculated from the sum of the partials:"<<std::endl;
248 out<<" E_g E1 E2 E3 E4 E5 M1 M2 M3 M4 M5 "<<std::endl;
249 for(G4int j=0;j<np[NShells];j++){
250 snprintf(word,1000,"%10.4g",Eg[NShells][j]); out<<word;
251 for(G4int k=0;k<ICC_NMULTIP;k++){
252 snprintf(word,1000," %10.4g",Icc_E[k][NShells][j]); out<<word;
253 }
254 for(G4int k=0;k<ICC_NMULTIP;k++){
255 snprintf(word,1000," %10.4g",Icc_M[k][NShells][j]); out<<word;
256 }
257 out<<std::endl;
258 }
259 out<<" ----------------------------------------------------------------------------------------------------------------------------------------"<<std::endl;
260 for(G4int i=0;i<NShells;i++){
261 out<<" ----------------------------------------------------------------------------------------------------------------------------------------"<<std::endl;
262 out<<" Binding energy = "<<BindingEnergy[i]<<" MeV - OrbitalName = "<<OrbitalName[i]<<" - np = "<<np[i]<<std::endl;
263 out<<" E_g E1 E2 E3 E4 E5 M1 M2 M3 M4 M5 "<<std::endl;
264 for(G4int j=0;j<np[i];j++){
265 snprintf(word,1000,"%10.4g",Eg[i][j]); out<<word;
266 for(G4int k=0;k<ICC_NMULTIP;k++){
267 snprintf(word,1000," %10.4g",Icc_E[k][i][j]); out<<word;
268 }
269 for(G4int k=0;k<ICC_NMULTIP;k++){
270 snprintf(word,1000," %10.4g",Icc_M[k][i][j]); out<<word;
271 }
272 out<<std::endl;
273 }
274 out<<" ----------------------------------------------------------------------------------------------------------------------------------------"<<std::endl;
275 }
276 out<<" ########################################################################################################################################## "<<std::endl;
277 }
278
279 void G4NuDEXInternalConversion::Init(const char* fname){
280
281 if(theZ<MINZINTABLES){ //then we have no info
282 return;
283 }
284
285 if(NShells!=0){ //Init only once
286 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
287 }
288
289 std::ifstream in(fname);
290 if(!in.good()){
291 std::cout<<" ################ Error opening "<<fname<<" ################"<<std::endl;
292 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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]>>word;
309 BindingEnergy[NShells]*=1.e-6; // all in MeV
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][100],Icc_M_tmp[ICC_NMULTIP][100];
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.substr(sz),&sz2); sz+=sz2;
344 }
345 for(G4int i=0;i<ICC_NMULTIP;i++){
346 Icc_M_tmp[i][np_tmp]=std::stof(word.substr(sz),&sz2); sz+=sz2;
347 }
348 if((G4int)(std::stof(word.substr(sz),&sz2)+0.01)!=theZ){
349 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
350 }
351 sz+=sz2;
352 sz2=word.find_first_not_of(' ',sz);
353 if(np_tmp==0){OrbitalName[orbindex]=word.substr(sz2,word.substr(sz2).size()-1);}
354 np_tmp++;
355 }
356 }
357 if(orbindex==0){break;}
358 //--------------------------------------------------------------------------------
359 }
360 }
361 }
362 if(!in.good()){
363 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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(__LINE__).c_str(),"##### Error in NuDEX #####");
377 }
378
379 //We evaluate it in the same frame as the total given by the data:
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],j+1,i);
400 Icc_M[j][NShells][k]+=GetICC(Eg[NShells][k],-j-1,i);
401 }
402 }
403 }
404
405 }
406
407 //if val>x[npmax] then ---> return 0
408 G4double G4NuDEXInternalConversion::Interpolate(G4double val,G4int npoints,G4double* x,G4double* y){
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_interp+1]-x[i_interp]);
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_interp+1]<=0 || x[i_interp]<=0){
427 //y=a0+a1*x
428 G4double a1=(y[i_interp+1]-y[i_interp])/(x[i_interp+1]-x[i_interp]);
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_interp])/std::log(x[i_interp+1]/x[i_interp]);
436 G4double a0=std::log(y[i_interp])-a1*std::log(x[i_interp]);
437
438 G4double result=std::exp(a0+a1*std::log(val));
439 return result;
440 }
441