| 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
43 #include "G4NuDEXRandom.hh"
44 #include "G4NuDEXLevelDensity.hh"
45
46
47
48 G4NuDEXLevelDensity::G4NuDEXLevelDensity(G4int
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(__L
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; Del
63 }
64
65
66 G4int G4NuDEXLevelDensity::ReadLDParameters(co
67
68 char fname[100];
69 if(LDType==1 || LDType==3){ // Back-Shifted-
70 snprintf(fname,100,"%s/LevelDensities/leve
71 }
72 else{ // Constant Temperature
73 snprintf(fname,100,"%s/LevelDensities/leve
74 }
75 G4double EL=0,EU=0;
76
77 std::ifstream in(fname);
78 if(!in.good()){
79 std::cout<<" ######## Error opening file "
80 NuDEXException(__FILE__,std::to_string(__L
81 }
82 G4int aZ,aA;
83 char word[200];
84 in.ignore(10000,'\n');
85
86 //std::cout<<" LDType="<<LDType<<" "<<fname
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>>
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>>
97 Ed=(EL+EU)/2.;
98 }
99 else{
100 NuDEXException(__FILE__,std::to_string(__LIN
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(defaultinput
115 }
116 if(inputfname!=0){
117 SearchLDParametersInInputFile(inputfname);
118 }
119
120 if(!HasData){ //no data
121 G4int check=CalculateLDParameters_BSFG(dir
122 if(check==0){
123 HasData=true;
124 if(LDType==2){
125 LDType=1;
126 std::cout<<" ##### WARNING: level density mo
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::CalculateLDParamete
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_p
150 if((Z_Int%2)==0 && ((A_Int-Z_Int)%2)==0){n_p
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,
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/shellc
162 std::ifstream in(fname);
163 if(!in.good()){
164 std::cout<<" ######## Error opening file "
165 NuDEXException(__FILE__,std::to_string(__L
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::SearchLDParametersI
194
195 if(inputfname!=0){
196 std::ifstream in2(inputfname);
197 if(!in2.good()){
198 std::cout<<" ############## Error openin
199 NuDEXException(__FILE__,std::to_string(_
200 }
201 std::string word_tmp;
202 while(in2>>word_tmp){
203 if(word_tmp[0]=='#'){in2.ignore(10000,'\
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>>De
209 }
210 else if(LDType==2){
211 in2>>dW_ldpar>>gamma_ldpar>>ainf_ldpar>>De
212 }
213 else if(LDType==3){
214 in2>>ainf_ldpar>>Delta_ldpar;
215 }
216 else{
217 std::cout<<" ############## Error: Unknown
218 NuDEXException(__FILE__,std::to_string(__L
219 }
220 if(!in2.good()){
221 std::cout<<" ############## Error reading
222 NuDEXException(__FILE__,std::to_string(__L
223 }
224 HasData=true;
225 break;
226 }
227 }
228 in2.close();
229 }
230
231 return 0;
232 }
233
234 void G4NuDEXLevelDensity::PrintParametersInInp
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<<" "<<ai
241 }
242 else if(LDType==2){
243 out<<dW_ldpar<<" "<<gamma_ldpar<<" "<<ai
244 }
245 else if(LDType==3){
246 out<<ainf_ldpar<<" "<<Delta_ldpar<<std::e
247 }
248 out.precision(oldprc);
249 out<<std::endl;
250
251 }
252
253
254 G4double G4NuDEXLevelDensity::GetNucleusTemper
255
256 if(!HasData){
257 NuDEXException(__FILE__,std::to_string(__L
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/Uva
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(
278
279 if(!HasData){
280 NuDEXException(__FILE__,std::to_string(__L
281 }
282
283 //If A_Int even/odd --> spinx2 (spin_val*2)
284 if(((A_Int+(G4int)(spin*2+0.01))%2)!=0 && (T
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 (
293 if(LDType==3){
294 G4double sig2=0.0888*std::pow(A_mass,2./3.
295 G4double sig=std::sqrt(sig2);
296 G4double rho=0.05893*std::exp(2.*std::sqrt
297 G4double xj2=(spin+0.5)*(spin+0.5);
298 G4double fj=(2.*spin+1.)/2./sig2*std::exp(
299 return 0.5*fj*rho;
300 }
301 //------------------------------------------
302
303
304 //------------------------------------------
305 //statistical factor from eq. 39 of RIPL-3 m
306 G4double Uval_Sn=Sn-Delta_ldpar;
307 G4double a_ldpar=ainf_ldpar*(1.+dW_ldpar/Uva
308 G4double a_ldpar_Sn=ainf_ldpar*(1.+dW_ldpar/
309 G4double sigma2_f=0.01389*std::pow(A_mass,5.
310 G4double sigma2_f_Sn=0.01389*std::pow(A_mass
311 G4double sigma2_d=(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)*(si
319 }
320 else{
321 sigma2=sigma2_f;
322 }
323 G4double statfactor=1./2.*(2.*spin+1.)/(2.*s
324 if(TotalLevelDensity==true){
325 statfactor=1;
326 }
327 //------------------------------------------
328
329 //CT + BSFG: Gilbert & Cameron, Can.J.Phys.
330 if(LDType==2 && ExcEnergy<Ex_ldpar){
331 G4double totalrho=std::exp((ExcEnergy-E0_l
332 return totalrho*statfactor;
333 }
334
335 //Else: BSFGM (LDType==1 or ExcEnergy>Ex_ldp
336 G4double rhotot_f=1./std::sqrt(2.*sigma2)/12
337 G4double rhotot_0=std::exp(1.)*a_ldpar/12./s
338 G4double totalrho=1./(1./rhotot_f+1./rhotot_
339
340 return totalrho*statfactor;
341
342 }
343
344
345 G4double G4NuDEXLevelDensity::EstimateInverse(
346
347 //We assume that rho is a monotonically incr
348
349 G4double tolerance=0.001; //the result will
350
351 G4double xmin=0;
352 G4double xmax=1;
353 while(GetLevelDensity(xmax,spin,parity)<LevD
354 xmax*=2;
355 }
356
357 while(xmin/xmax<1-tolerance){
358 G4double x0=(xmin+xmax)/2.;
359 if(GetLevelDensity(x0,spin,parity)<LevDen_
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(G4doub
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:
388
389 out<<" Level density type: "<<LDType<<std::e
390 if(LDType==1){ // Back-Shifted-Fermi-Gas mod
391 out<<" ainf = "<<ainf_ldpar<<" gamma = "<
392 }
393 else{
394 out<<" ainf = "<<ainf_ldpar<<" gamma = "<
395 }
396
397 }
398
399
400
401
402