| 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
43 #include "G4NuDEXRandom.hh"
44 #include "G4NuDEXLevelDensity.hh"
45
46
47
48 G4NuDEXLevelDensity::G4NuDEXLevelDensity(G4int aZ,G4int aA,G4int ldtype){
49
50 Z_Int=aZ;
51 A_Int=aA;
52 LDType=ldtype;
53 if(LDType<0){LDType=DEFAULTLDTYPE;}
54 A_mass=A_Int;
55 if(LDType!=1 && LDType!=2 && LDType!=3){
56 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
57 }
58 HasData=false;
59
60 Sn=-1; D0=-1; I0=-1000;
61 Ed=0;
62 ainf_ldpar=0; gamma_ldpar=0; dW_ldpar=0; Delta_ldpar=0; T_ldpar=0; E0_ldpar=0; Ex_ldpar=0;
63 }
64
65
66 G4int G4NuDEXLevelDensity::ReadLDParameters(const char* dirname,const char* inputfname,const char* defaultinputfname){
67
68 char fname[100];
69 if(LDType==1 || LDType==3){ // Back-Shifted-Fermi-Gas model
70 snprintf(fname,100,"%s/LevelDensities/level-densities-bfmeff.dat",dirname);
71 }
72 else{ // Constant Temperature
73 snprintf(fname,100,"%s/LevelDensities/level-densities-ctmeff.dat",dirname);
74 }
75 G4double EL=0,EU=0;
76
77 std::ifstream in(fname);
78 if(!in.good()){
79 std::cout<<" ######## Error opening file "<<fname<<" ########"<<std::endl;
80 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
81 }
82 G4int aZ,aA;
83 char word[200];
84 in.ignore(10000,'\n');
85
86 //std::cout<<" LDType="<<LDType<<" "<<fname<<" "<<Z_Int<<" "<<A_Int<<std::endl;
87
88 while(in>>aZ>>aA){
89 if(aZ==Z_Int && aA==A_Int){
90 if(LDType==1 || LDType==3){
91 in>>word>>I0>>Sn>>D0>>word>>word>>EL>>word>>EU>>dW_ldpar>>gamma_ldpar>>ainf_ldpar>>word>>Delta_ldpar;
92 Ex_ldpar=0; E0_ldpar=0; T_ldpar=0;
93 Ed=(EL+EU)/2.;
94 }
95 else if(LDType==2){
96 in>>word>>I0>>Sn>>D0>>word>>word>>EL>>word>>EU>>dW_ldpar>>gamma_ldpar>>ainf_ldpar>>word>>Delta_ldpar>>Ex_ldpar>>E0_ldpar>>T_ldpar;
97 Ed=(EL+EU)/2.;
98 }
99 else{
100 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
101 }
102 if(in.good()){
103 HasData=true;
104 break;
105 }
106 }
107 in.ignore(10000,'\n');
108 }
109 in.close();
110
111
112 //Re-write some parameters if inputfname!=0
113 if(defaultinputfname!=0){
114 SearchLDParametersInInputFile(defaultinputfname);
115 }
116 if(inputfname!=0){
117 SearchLDParametersInInputFile(inputfname);
118 }
119
120 if(!HasData){ //no data
121 G4int check=CalculateLDParameters_BSFG(dirname);
122 if(check==0){
123 HasData=true;
124 if(LDType==2){
125 LDType=1;
126 std::cout<<" ##### WARNING: level density model for ZA="<<Z_Int*1000+A_Int<<" changed to Back-Shifted-Fermi-Gas model #####"<<std::endl;
127 }
128 }
129 }
130
131 if(HasData){return 0;}
132
133
134 //else, some problem reading ...
135 return -1;
136 }
137
138
139 G4int G4NuDEXLevelDensity::CalculateLDParameters_BSFG(const char* dirname){
140
141 //Eq. 61 of RIPL-3:
142 G4double alpha=0.0722396; //MeV^{-1}
143 G4double beta= 0.195267; //MeV^{-1}
144 G4double gamma0=0.410289; //MeV^{-1}
145 G4double delta=0.173015; //MeV
146
147 //Delta_ldpar: Eq. 50 of RIPL-3:
148 G4double n_par=0;
149 if((Z_Int%2)==1 && ((A_Int-Z_Int)%2)==1){n_par=-1;} //odd-odd (impar-impar)
150 if((Z_Int%2)==0 && ((A_Int-Z_Int)%2)==0){n_par=1;} //even-even (par-par)
151 Delta_ldpar=n_par*12/std::sqrt(A_mass)+delta;
152
153 //ainf_ldpar: Eq. 52 of RIPL-3:
154 ainf_ldpar=alpha*A_Int+beta*std::pow(A_mass,2./3.);
155
156 //gamma_ldpar: Eq. 53 of RIPL-3:
157 gamma_ldpar=gamma0/std::pow(A_mass,1./3.);
158
159 //dW_ldpar --> from data file
160 char fname[100];
161 snprintf(fname,100,"%s/LevelDensities/shellcor-ms.dat",dirname);
162 std::ifstream in(fname);
163 if(!in.good()){
164 std::cout<<" ######## Error opening file "<<fname<<" ########"<<std::endl;
165 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
166 }
167 G4int aZ,aA;
168 char word[200];
169 in.ignore(10000,'\n');
170 in.ignore(10000,'\n');
171 in.ignore(10000,'\n');
172 in.ignore(10000,'\n');
173 while(in>>aZ>>aA){
174 if(aZ==Z_Int && aA==A_Int){
175 in>>word>>dW_ldpar;
176 if(in.good()){break;}
177 }
178 in.ignore(10000,'\n');
179 }
180 if(!in.good()){//no data found
181 return -1;
182 }
183 in.close();
184
185 Ex_ldpar=0; E0_ldpar=0; T_ldpar=0;
186 Ed=0;
187
188
189 return 0;
190 }
191
192
193 G4int G4NuDEXLevelDensity::SearchLDParametersInInputFile(const char* inputfname){
194
195 if(inputfname!=0){
196 std::ifstream in2(inputfname);
197 if(!in2.good()){
198 std::cout<<" ############## Error opening "<<inputfname<<" ##############"<<std::endl;
199 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
200 }
201 std::string word_tmp;
202 while(in2>>word_tmp){
203 if(word_tmp[0]=='#'){in2.ignore(10000,'\n');}
204 if(word_tmp==std::string("END")){break;}
205 if(word_tmp==std::string("LDPARAMETERS")){
206 in2>>LDType;
207 if(LDType==1){
208 in2>>dW_ldpar>>gamma_ldpar>>ainf_ldpar>>Delta_ldpar;
209 }
210 else if(LDType==2){
211 in2>>dW_ldpar>>gamma_ldpar>>ainf_ldpar>>Delta_ldpar>>Ex_ldpar>>E0_ldpar>>T_ldpar;
212 }
213 else if(LDType==3){
214 in2>>ainf_ldpar>>Delta_ldpar;
215 }
216 else{
217 std::cout<<" ############## Error: Unknown LDType="<<LDType<<" in "<<inputfname<<" ##############"<<std::endl;
218 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
219 }
220 if(!in2.good()){
221 std::cout<<" ############## Error reading "<<inputfname<<" ##############"<<std::endl;
222 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
223 }
224 HasData=true;
225 break;
226 }
227 }
228 in2.close();
229 }
230
231 return 0;
232 }
233
234 void G4NuDEXLevelDensity::PrintParametersInInputFileFormat(std::ostream &out){
235
236 out<<"LDPARAMETERS"<<std::endl;
237 out<<LDType<<std::endl;
238 G4long oldprc = out.precision(15);
239 if(LDType==1){
240 out<<dW_ldpar<<" "<<gamma_ldpar<<" "<<ainf_ldpar<<" "<<Delta_ldpar<<std::endl;
241 }
242 else if(LDType==2){
243 out<<dW_ldpar<<" "<<gamma_ldpar<<" "<<ainf_ldpar<<" "<<Delta_ldpar<<" "<<Ex_ldpar<<" "<<E0_ldpar<<" "<<T_ldpar<<std::endl;
244 }
245 else if(LDType==3){
246 out<<ainf_ldpar<<" "<<Delta_ldpar<<std::endl;
247 }
248 out.precision(oldprc);
249 out<<std::endl;
250
251 }
252
253
254 G4double G4NuDEXLevelDensity::GetNucleusTemperature(G4double ExcEnergy){
255
256 if(!HasData){
257 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
258 }
259
260 if(ExcEnergy<Ex_ldpar && LDType==2){
261 return T_ldpar;
262 }
263
264 G4double Uval=ExcEnergy-Delta_ldpar;
265 if(Uval<=0){return 0;}
266 G4double a_ldpar=ainf_ldpar*(1.+dW_ldpar/Uval*(1.-std::exp(-gamma_ldpar*Uval)));
267 if(LDType==3){
268 a_ldpar=ainf_ldpar;
269 }
270 return std::sqrt(Uval/a_ldpar);
271
272
273 }
274
275
276 //Gilbert-Cameron:
277 G4double G4NuDEXLevelDensity::GetLevelDensity(G4double ExcEnergy,G4double spin,G4bool ,G4bool TotalLevelDensity){
278
279 if(!HasData){
280 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
281 }
282
283 //If A_Int even/odd --> spinx2 (spin_val*2) should be even/odd
284 if(((A_Int+(G4int)(spin*2+0.01))%2)!=0 && (TotalLevelDensity==false)){
285 return 0;
286 }
287
288 G4double Uval=ExcEnergy-Delta_ldpar;
289 if(Uval<0){Uval=1.e-6;}
290
291 //----------------------------------------------------------------
292 // Back shifted: von Egidy et al., NP A481 (1988) 189
293 if(LDType==3){
294 G4double sig2=0.0888*std::pow(A_mass,2./3.)*std::sqrt(ainf_ldpar*Uval);
295 G4double sig=std::sqrt(sig2);
296 G4double rho=0.05893*std::exp(2.*std::sqrt(ainf_ldpar*Uval))/sig/std::pow(ainf_ldpar,0.25)/std::pow(Uval,1.25);
297 G4double xj2=(spin+0.5)*(spin+0.5);
298 G4double fj=(2.*spin+1.)/2./sig2*std::exp(-xj2/2./sig2);
299 return 0.5*fj*rho;
300 }
301 //----------------------------------------------------------------
302
303
304 //--------------------------------------------------------------------------------
305 //statistical factor from eq. 39 of RIPL-3 manual, and sigma2 from eqs. 57-60
306 G4double Uval_Sn=Sn-Delta_ldpar;
307 G4double a_ldpar=ainf_ldpar*(1.+dW_ldpar/Uval*(1.-std::exp(-gamma_ldpar*Uval)));
308 G4double a_ldpar_Sn=ainf_ldpar*(1.+dW_ldpar/Uval_Sn*(1.-std::exp(-gamma_ldpar*Uval_Sn)));
309 G4double sigma2_f=0.01389*std::pow(A_mass,5./3.)/ainf_ldpar*std::sqrt(a_ldpar*Uval);
310 G4double sigma2_f_Sn=0.01389*std::pow(A_mass,5./3.)/ainf_ldpar*std::sqrt(a_ldpar_Sn*Uval);
311 G4double sigma2_d=(0.83*std::pow(A_mass,0.26))*(0.83*std::pow(A_mass,0.26));
312
313 G4double sigma2;
314 if(ExcEnergy<=Ed){
315 sigma2=sigma2_d;//if ExcEnergy<Ed
316 }
317 else if(ExcEnergy<=Sn){
318 sigma2=sigma2_d+(ExcEnergy-Ed)/(Sn-Ed)*(sigma2_f_Sn-sigma2_d);
319 }
320 else{
321 sigma2=sigma2_f;
322 }
323 G4double statfactor=1./2.*(2.*spin+1.)/(2.*sigma2)*std::exp(-(spin+1/2.)*(spin+1/2.)/2./sigma2);
324 if(TotalLevelDensity==true){
325 statfactor=1;
326 }
327 //--------------------------------------------------------------------------------
328
329 //CT + BSFG: Gilbert & Cameron, Can.J.Phys. 43 (1965) 1446
330 if(LDType==2 && ExcEnergy<Ex_ldpar){
331 G4double totalrho=std::exp((ExcEnergy-E0_ldpar)/T_ldpar)/T_ldpar;
332 return totalrho*statfactor;
333 }
334
335 //Else: BSFGM (LDType==1 or ExcEnergy>Ex_ldpar)
336 G4double rhotot_f=1./std::sqrt(2.*sigma2)/12.*std::exp(2.*std::sqrt(a_ldpar*Uval))/std::pow(a_ldpar,1./4.)/std::pow(Uval,5./4.);
337 G4double rhotot_0=std::exp(1.)*a_ldpar/12./std::sqrt(sigma2)*std::exp(a_ldpar*Uval);
338 G4double totalrho=1./(1./rhotot_f+1./rhotot_0);
339
340 return totalrho*statfactor;
341
342 }
343
344
345 G4double G4NuDEXLevelDensity::EstimateInverse(G4double LevDen_iMeV,G4double spin,G4bool parity){
346
347 //We assume that rho is a monotonically increasing function
348
349 G4double tolerance=0.001; //the result will have this relative tolerance. 0.01 means 1%
350
351 G4double xmin=0;
352 G4double xmax=1;
353 while(GetLevelDensity(xmax,spin,parity)<LevDen_iMeV){
354 xmax*=2;
355 }
356
357 while(xmin/xmax<1-tolerance){
358 G4double x0=(xmin+xmax)/2.;
359 if(GetLevelDensity(x0,spin,parity)<LevDen_iMeV){
360 xmin=x0;
361 }
362 else{
363 xmax=x0;
364 }
365 }
366
367 return (xmin+xmax)/2.;
368
369 }
370
371
372 G4double G4NuDEXLevelDensity::Integrate(G4double Emin,G4double Emax,G4double spin,G4bool parity){
373
374 G4int nb=1000;
375 G4double Integral=0,x1,x2,y1,y2;
376 for(G4int i=0;i<nb;i++){
377 x1=Emin+(Emax-Emin)*i/(G4double)(nb-1.);
378 x2=Emin+(Emax-Emin)*(i+1.)/(G4double)(nb-1.);
379 y1=GetLevelDensity(x1,spin,parity);
380 y2=GetLevelDensity(x2,spin,parity);
381 Integral+=(y1+y2)/2.*(x2-x1);
382 }
383
384 return Integral;
385 }
386
387 void G4NuDEXLevelDensity::PrintParameters(std::ostream &out){
388
389 out<<" Level density type: "<<LDType<<std::endl;
390 if(LDType==1){ // Back-Shifted-Fermi-Gas model
391 out<<" ainf = "<<ainf_ldpar<<" gamma = "<<gamma_ldpar<<" dW = "<<dW_ldpar<<" Delta = "<<Delta_ldpar<<std::endl;
392 }
393 else{
394 out<<" ainf = "<<ainf_ldpar<<" gamma = "<<gamma_ldpar<<" dW = "<<dW_ldpar<<" Delta = "<<Delta_ldpar<<" T = "<<T_ldpar<<" E0 = "<<E0_ldpar<<" Ex = "<<Ex_ldpar<<std::endl;
395 }
396
397 }
398
399
400
401
402