| 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 "G4NuDEXRandom.hh"
43 #include "G4NuDEXLevelDensity.hh"
44 #include "G4NuDEXPSF.hh"
45
46
47 G4NuDEXPSF::G4NuDEXPSF(G4int aZ,G4int aA){
48 Z_Int=aZ;
49 A_Int=aA;
50 nR_E1= nR_M1= nR_E2=0;
51 np_E1= np_M1= np_E2=0;
52 x_E1= y_E1=nullptr;
53 x_M1= y_M1=nullptr;
54 x_E2= y_E2=nullptr;
55 E1_normFac=-1; M1_normFac=-1; E2_normFac=-1;
56 NormEmin=0; NormEmax=6; //Integral between 0
57 ScaleFactor_E1=1;
58 ScaleFactor_M1=1;
59 ScaleFactor_E2=1;
60 theLD=nullptr;
61 }
62
63 G4NuDEXPSF::~G4NuDEXPSF(){
64 delete [] x_E1;
65 delete [] y_E1;
66 delete [] x_M1;
67 delete [] y_M1;
68 delete [] x_E2;
69 delete [] y_E2;
70 }
71
72
73 //If inputfname!=0 then we take the PSF data f
74 G4int G4NuDEXPSF::Init(const char* dirname,G4N
75
76 theLD=aLD;
77
78 //Three options: very detailed model, if not
79
80 char fname[500];
81 G4bool IsDone=false;
82
83 //input:
84 if(inputfname!=0){
85 IsDone=TakePSFFromInputFile(inputfname);
86 if(IsDone){return 0;}
87 }
88
89 //default input:
90 if(defaultinputfname!=0){
91 IsDone=TakePSFFromInputFile(defaultinputfn
92 if(IsDone){return 0;}
93 }
94
95 //Detailed model
96 snprintf(fname,500,"%s/PSF/PSF_param.dat",di
97 IsDone=TakePSFFromDetailedParFile(fname);
98 if(IsDone){return 0;}
99
100 //IAEA - 2019 values:
101 if(PSFflag==0){
102 snprintf(fname,500,"%s/PSF/CRP_IAEA_SMLO_E
103 IsDone=TakePSFFromIAEA01(fname);
104 if(IsDone){return 0;}
105 }
106 else if(PSFflag!=1){
107 NuDEXException(__FILE__,std::to_string(__L
108 }
109
110 //RIPL-MLO values:
111 snprintf(fname,500,"%s/PSF/gdr-parameters&er
112 IsDone=TakePSFFromRIPL01(fname);
113 if(IsDone){return 0;}
114
115 //RIPL-Theorethical values:
116 snprintf(fname,500,"%s/PSF/gdr-parameters-th
117 IsDone=TakePSFFromRIPL02(fname);
118 if(IsDone){return 0;}
119
120 //Theorethical values:
121 // E1 for spherical nucleus:
122 nR_E1=0;
123 PSFType_E1[nR_E1]=2;
124 //G4double a=31.2,b=20.6,c=0.026,d=1.05; //S
125 //G4double a=27.47,b=22.063,c=0.0277,d=1.222
126 G4double a=28.69,b=21.731,c=0.0285,d=1.267;/
127
128 E_E1[nR_E1]=a*std::pow(A_Int,-1./3.)+b*std::
129 G_E1[nR_E1]=c*std::pow(E_E1[nR_E1],1.9);
130 s_E1[nR_E1]=120/3.141592*d*(A_Int-Z_Int)*Z_I
131 nR_E1++;
132 GenerateM1AndE2FromE1();
133
134 return 0;
135 }
136
137 void G4NuDEXPSF::GenerateM1AndE2FromE1(){
138
139 //M1:
140 nR_M1=0;
141 E_M1[nR_M1]=41*std::pow(A_Int,-1./3.);
142 G_M1[nR_M1]=4;
143 s_M1[nR_M1]=1;
144 PSFType_M1[nR_M1]=0;
145 nR_M1++;
146
147 //f(E1)/f(M1) = 0.0588*A**0.878 at +-7 Me
148 G4double fE1=GetE1(7,7);
149 G4double fM1=GetM1(7,7);
150 s_M1[0]=fE1/0.0588/std::pow(A_Int,0.878)/fM1
151
152
153 //E2:
154 nR_E2=0;
155 E_E2[nR_E2]=63*std::pow(A_Int,-1./3.);
156 G_E2[nR_E2]=6.11-0.021*A_Int;
157 s_E2[nR_E2]=0.00014*Z_Int*Z_Int*E_E2[nR_E2]/
158 PSFType_E2[nR_E2]=0;
159 nR_E2++;
160
161 }
162
163 G4bool G4NuDEXPSF::TakePSFFromRIPL02(const cha
164
165 G4bool result=false;
166 G4int aA,aZ;
167 std::ifstream in(fname);
168 char dum[200];
169
170 for(G4int i=0;i<4;i++){in.ignore(10000,'\n')
171 while(in>>aZ>>aA){
172 if(aZ==Z_Int && aA==A_Int){
173 result=true;
174 in>>dum>>dum;
175 nR_E1=2;
176 in>>E_E1[0]>>G_E1[0]>>E_E1[1]>>G_E1[1];
177 PSFType_E1[0]=2; PSFType_E1[1]=2; //SMLO
178
179 G4double a=28.69,b=21.731,c=0.0285,d=1.2
180 G4double E_E1_0=a*std::pow(A_Int,-1./3.)
181 G4double G_E1_0=c*std::pow(E_E1_0,1.9);
182 G4double s_E1_0=120/3.141592*d*(A_Int-Z_
183 s_E1[0]=s_E1_0/3.;
184 s_E1[1]=2.*s_E1_0/3.;
185
186 break;
187 }
188 in.ignore(10000,'\n');
189 }
190 in.close();
191 if(result){GenerateM1AndE2FromE1();}
192
193 return result;
194 }
195
196
197 G4bool G4NuDEXPSF::TakePSFFromRIPL01(const cha
198
199 G4bool result=false;
200 G4int aA,aZ;
201 std::ifstream in(fname);
202 char dum[200];
203
204 for(G4int i=0;i<7;i++){in.ignore(10000,'\n')
205 while(in>>aZ>>aA){
206 if(aZ==Z_Int && aA==A_Int){
207 result=true;
208 in>>dum>>dum;
209 nR_E1=0;
210 in>>E_E1[nR_E1]>>s_E1[nR_E1]>>G_E1[nR_E1
211 PSFType_E1[nR_E1]=2; //SMLO
212 nR_E1++;
213 //sometimes there is a second resonance:
214 in>>E_E1[nR_E1]>>dum>>G_E1[nR_E1];
215 if(dum[0]!='-'){ //there is a second res
216 s_E1[nR_E1]=std::atof(dum);
217 PSFType_E1[nR_E1]=2;
218 nR_E1++;
219 }
220 break;
221 }
222 in.ignore(10000,'\n');
223 }
224 in.close();
225 if(result){GenerateM1AndE2FromE1();}
226
227 return result;
228 }
229
230
231 G4bool G4NuDEXPSF::TakePSFFromIAEA01(const cha
232
233 G4bool result=false;
234 G4int aA,aZ;
235 char dum[200];
236 G4double beta=0;
237 std::ifstream in(fname);
238 while(in>>aZ>>aA){
239 if(aZ==Z_Int && aA==A_Int){
240 result=true;
241 nR_E1=0;
242 in>>dum>>dum>>E_E1[nR_E1]>>dum>>dum>>G_E
243 PSFType_E1[nR_E1]=11;
244 nR_E1++;
245 in>>dum;
246 if(std::string(dum)==std::string("beta="
247 in>>beta;
248 break;
249 }
250 else if(std::string(dum)==std::string("E
251 in>>dum>>E_E1[nR_E1]>>dum>>dum>>G_E1[nR_E1]>
252 PSFType_E1[nR_E1]=11;
253 nR_E1++;
254
255 }
256 else{
257 NuDEXException(__FILE__,std::to_string(__LIN
258 }
259 break;
260 }
261 in.ignore(10000,'\n');
262 }
263 if(!result){
264 return result;
265 }
266
267 //------------------------------------------
268 //Now M1 (https://doi.org/10.1140/epja/i2019
269 nR_M1=0;
270 //Spin-flip:
271 PSFType_M1[nR_M1]=0;
272 E_M1[nR_M1]=18.0*std::pow(A_Int,-1./6.);
273 G_M1[nR_M1]=4;
274 s_M1[nR_M1]=0.03*std::pow(A_Int,5./6.);
275 nR_M1++;
276 //Scissors-mode:
277 PSFType_M1[nR_M1]=0;
278 E_M1[nR_M1]=5.0*std::pow(A_Int,-1./10.);
279 G_M1[nR_M1]=1.5;
280 s_M1[nR_M1]=0.02*std::fabs(beta)*std::pow(A_
281 nR_M1++;
282 //upbend:
283 PSFType_M1[nR_M1]=21;
284 E_M1[nR_M1]=0.4035*std::exp(-6.0*std::fabs(b
285 G_M1[nR_M1]=0.8;
286 s_M1[nR_M1]=0;
287 nR_M1++;
288 //------------------------------------------
289
290 //------------------------------------------
291 //E2 same as in the old RIPL recommendations
292 nR_E2=0;
293 E_E2[nR_E2]=63*std::pow(A_Int,-1./3.);
294 G_E2[nR_E2]=6.11-0.021*A_Int;
295 s_E2[nR_E2]=0.00014*Z_Int*Z_Int*E_E2[nR_E2]/
296 PSFType_E2[nR_E2]=0;
297 nR_E2++;
298 //------------------------------------------
299
300 return result;
301 }
302
303
304
305 G4bool G4NuDEXPSF::TakePSFFromInputFile(const
306
307 G4bool result=false;
308 char word[1000];
309 std::ifstream in(fname);
310 while(in>>word){
311 if(word[0]=='#'){in.ignore(10000,'\n');}
312 if(std::string(word)==std::string("END")){
313 if(std::string(word)==std::string("PSF")){
314 result=true;
315 in>>nR_E1;
316 for(G4int i=0;i<nR_E1;i++){
317 in>>PSFType_E1[i]>>E_E1[i]>>G_E1[i]>>s_E1[i]
318 if(PSFType_E1[i]==7){in>>p1_E1[i];}
319 if(PSFType_E1[i]==8){in>>p1_E1[i]>>p2_E1[i];
320 if(PSFType_E1[i]==9){in>>p1_E1[i]>>p2_E1[i];
321 if(PSFType_E1[i]==10){in>>p1_E1[i]>>p2_E1[i]
322 if(PSFType_E1[i]==40 || PSFType_E1[i]==41){
323 if(x_E1!=0){NuDEXException(__FILE__,std::t
324 in>>np_E1;
325 x_E1=new G4double[np_E1]; y_E1=new G4doubl
326 for(G4int j=0;j<np_E1;j++){in>>x_E1[j]>>y_
327 in>>E1_normFac;
328 }
329 }
330 in>>nR_M1;
331 for(G4int i=0;i<nR_M1;i++){
332 in>>PSFType_M1[i]>>E_M1[i]>>G_M1[i]>>s_M1[i]
333 if(PSFType_M1[i]==7){in>>p1_M1[i];}
334 if(PSFType_M1[i]==8){in>>p1_M1[i]>>p2_M1[i];
335 if(PSFType_M1[i]==9){in>>p1_M1[i]>>p2_M1[i];
336 if(PSFType_M1[i]==10){in>>p1_M1[i]>>p2_M1[i]
337 if(PSFType_M1[i]==40 || PSFType_M1[i]==41){/
338 if(x_M1!=0){NuDEXException(__FILE__,std::t
339 in>>np_M1;
340 x_M1=new G4double[np_M1]; y_M1=new G4doubl
341 for(G4int j=0;j<np_M1;j++){in>>x_M1[j]>>y_
342 in>>M1_normFac;
343 }
344 }
345 in>>nR_E2;
346 for(G4int i=0;i<nR_E2;i++){
347 in>>PSFType_E2[i]>>E_E2[i]>>G_E2[i]>>s_E2[i]
348 if(PSFType_E2[i]==7){in>>p1_E2[i];}
349 if(PSFType_E2[i]==8){in>>p1_E2[i]>>p2_E2[i];
350 if(PSFType_E2[i]==9){in>>p1_E2[i]>>p2_E2[i];
351 if(PSFType_E2[i]==10){in>>p1_E2[i]>>p2_E2[i]
352 if(PSFType_E2[i]==40 || PSFType_E2[i]==41){/
353 if(x_E2!=0){NuDEXException(__FILE__,std::t
354 in>>np_E2;
355 x_E2=new G4double[np_E2]; y_E2=new G4doubl
356 for(G4int j=0;j<np_E2;j++){in>>x_E2[j]>>y_
357 in>>E2_normFac;
358 }
359 }
360 break;
361 }
362 }
363
364 Renormalize(); // if XX_normFac>0 --> renorm
365
366 return result;
367 }
368
369
370 void G4NuDEXPSF::Renormalize(){
371
372 G4int npIntegral=1000;
373 G4double Integral=0,x_eval,y_eval;
374 G4double binWidth=(NormEmax-NormEmin)/npInte
375
376 //------------------------------------------
377 if(E1_normFac>0){
378 Integral=0;
379 for(G4int i=0;i<npIntegral;i++){
380 x_eval=NormEmin+binWidth*(i+0.5);
381 y_eval=GetE1(x_eval,NormEmax);
382 Integral+=y_eval;
383 }
384 Integral*=binWidth;
385 ScaleFactor_E1=E1_normFac/Integral;
386 }
387 //------------------------------------------
388 //------------------------------------------
389 if(M1_normFac>0){
390 Integral=0;
391 for(G4int i=0;i<npIntegral;i++){
392 x_eval=NormEmin+binWidth*(i+0.5);
393 y_eval=GetM1(x_eval,NormEmax);
394 Integral+=y_eval;
395 }
396 Integral*=binWidth;
397 ScaleFactor_M1=M1_normFac/Integral;
398 //std::cout<<M1_normFac<<" "<<Integral<<"
399 }
400 //------------------------------------------
401 //------------------------------------------
402 if(E2_normFac>0){
403 Integral=0;
404 for(G4int i=0;i<npIntegral;i++){
405 x_eval=NormEmin+binWidth*(i+0.5);
406 y_eval=GetE2(x_eval,NormEmax);
407 Integral+=y_eval;
408 }
409 Integral*=binWidth;
410 ScaleFactor_E2=E2_normFac/Integral;
411 }
412 //------------------------------------------
413
414 }
415
416
417
418 G4bool G4NuDEXPSF::TakePSFFromDetailedParFile(
419
420 G4bool result=false;
421 G4int aA,aZ;
422 std::ifstream in(fname);
423 while(in>>aZ>>aA){
424 if(aZ==Z_Int && aA==A_Int){
425 result=true;
426 in>>nR_E1;
427 for(G4int i=0;i<nR_E1;i++){
428 in>>PSFType_E1[i]>>E_E1[i]>>G_E1[i]>>s_E1[i]
429 if(PSFType_E1[i]==7){in>>p1_E1[i];}
430 if(PSFType_E1[i]==8){in>>p1_E1[i]>>p2_E1[i];
431 if(PSFType_E1[i]==9){in>>p1_E1[i]>>p2_E1[i];
432 if(PSFType_E1[i]==10){in>>p1_E1[i]>>p2_E1[i]
433 }
434 in>>nR_M1;
435 for(G4int i=0;i<nR_M1;i++){
436 in>>PSFType_M1[i]>>E_M1[i]>>G_M1[i]>>s_M1[i]
437 if(PSFType_M1[i]==7){in>>p1_M1[i];}
438 if(PSFType_M1[i]==8){in>>p1_M1[i]>>p2_M1[i];
439 if(PSFType_M1[i]==9){in>>p1_M1[i]>>p2_M1[i];
440 if(PSFType_M1[i]==10){in>>p1_M1[i]>>p2_M1[i]
441 }
442 in>>nR_E2;
443 for(G4int i=0;i<nR_E2;i++){
444 in>>PSFType_E2[i]>>E_E2[i]>>G_E2[i]>>s_E2[i]
445 if(PSFType_E2[i]==7){in>>p1_E2[i];}
446 if(PSFType_E2[i]==8){in>>p1_E2[i]>>p2_E2[i];
447 if(PSFType_E2[i]==9){in>>p1_E2[i]>>p2_E2[i];
448 if(PSFType_E2[i]==10){in>>p1_E2[i]>>p2_E2[i]
449 }
450 break;
451 }
452 in.ignore(10000,'\n');
453 }
454 in.close();
455
456 return result;
457 }
458
459
460
461
462
463 G4double G4NuDEXPSF::GetE1(G4double Eg,G4doubl
464
465 G4double result=0;
466 for(G4int i=0;i<nR_E1;i++){
467 if(PSFType_E1[i]==0){
468 result+=8.674E-8*SLO(Eg,E_E1[i],G_E1[i],
469 }
470 else if(PSFType_E1[i]==1){
471 result+=8.674E-8*EGLO(Eg,E_E1[i],G_E1[i]
472 }
473 else if(PSFType_E1[i]==2){
474 result+=8.674E-8*SMLO(Eg,E_E1[i],G_E1[i]
475 }
476 else if(PSFType_E1[i]==3){
477 result+=8.674E-8*GLO(Eg,E_E1[i],G_E1[i],
478 }
479 else if(PSFType_E1[i]==4){
480 result+=8.674E-8*MGLO(Eg,E_E1[i],G_E1[i]
481 }
482 else if(PSFType_E1[i]==5){
483 result+=8.674E-8*KMF(Eg,E_E1[i],G_E1[i],
484 }
485 else if(PSFType_E1[i]==6){
486 result+=8.674E-8*GH(Eg,E_E1[i],G_E1[i],s
487 }
488 else if(PSFType_E1[i]==7){
489 result+=8.674E-8*MEGLO(Eg,E_E1[i],G_E1[i
490 }
491 else if(PSFType_E1[i]==8){
492 result+=8.674E-8*MEGLO(Eg,E_E1[i],G_E1[i
493 }
494 else if(PSFType_E1[i]==9){
495 result+=8.674E-8*MEGLO(Eg,E_E1[i],G_E1[i
496 }
497 else if(PSFType_E1[i]==10){
498 result+=8.674E-8*MEGLO(Eg,E_E1[i],G_E1[i
499 }
500 else if(PSFType_E1[i]==11){
501 result+=8.674E-8*SMLO_v2(Eg,E_E1[i],G_E1
502 }
503 else if(PSFType_E1[i]==20){
504 result+=8.674E-8*Gauss(Eg,E_E1[i],G_E1[i
505 }
506 else if(PSFType_E1[i]==21){
507 result+=8.674E-8*Expo(Eg,E_E1[i],G_E1[i]
508 }
509 else if(PSFType_E1[i]==40){
510 result+=EvaluateFunction(Eg,np_E1,x_E1,y
511 }
512 else if(PSFType_E1[i]==41){
513 result+=std::pow(10.,EvaluateFunction(Eg
514 }
515 else{
516 NuDEXException(__FILE__,std::to_string(_
517 }
518 }
519
520 if(result!=result){ // nan
521 NuDEXException(__FILE__,std::to_string(__L
522 }
523
524 return result*ScaleFactor_E1;
525 }
526
527 G4double G4NuDEXPSF::GetM1(G4double Eg,G4doubl
528
529 G4double result=0;
530 for(G4int i=0;i<nR_M1;i++){
531 if(PSFType_M1[i]==0){
532 result+=8.674E-8*SLO(Eg,E_M1[i],G_M1[i],
533 }
534 else if(PSFType_M1[i]==1){
535 result+=8.674E-8*EGLO(Eg,E_M1[i],G_M1[i]
536 }
537 else if(PSFType_M1[i]==2){
538 result+=8.674E-8*SMLO(Eg,E_M1[i],G_M1[i]
539 }
540 else if(PSFType_M1[i]==3){
541 result+=8.674E-8*GLO(Eg,E_M1[i],G_M1[i],
542 }
543 else if(PSFType_M1[i]==4){
544 result+=8.674E-8*MGLO(Eg,E_M1[i],G_M1[i]
545 }
546 else if(PSFType_M1[i]==5){
547 result+=8.674E-8*KMF(Eg,E_M1[i],G_M1[i],
548 }
549 else if(PSFType_M1[i]==6){
550 result+=8.674E-8*GH(Eg,E_M1[i],G_M1[i],s
551 }
552 else if(PSFType_M1[i]==7){
553 result+=8.674E-8*MEGLO(Eg,E_M1[i],G_M1[i
554 }
555 else if(PSFType_M1[i]==8){
556 result+=8.674E-8*MEGLO(Eg,E_M1[i],G_M1[i
557 }
558 else if(PSFType_M1[i]==9){
559 result+=8.674E-8*MEGLO(Eg,E_M1[i],G_M1[i
560 }
561 else if(PSFType_M1[i]==10){
562 result+=8.674E-8*MEGLO(Eg,E_M1[i],G_M1[i
563 }
564 else if(PSFType_M1[i]==11){
565 result+=8.674E-8*SMLO_v2(Eg,E_M1[i],G_M1
566 }
567 else if(PSFType_M1[i]==20){
568 result+=8.674E-8*Gauss(Eg,E_M1[i],G_M1[i
569 }
570 else if(PSFType_M1[i]==21){
571 result+=8.674E-8*Expo(Eg,E_M1[i],G_M1[i]
572 }
573 else if(PSFType_M1[i]==40){
574 result+=EvaluateFunction(Eg,np_M1,x_M1,y
575 }
576 else if(PSFType_M1[i]==41){
577 result+=std::pow(10.,EvaluateFunction(Eg
578 }
579 else{
580 NuDEXException(__FILE__,std::to_string(_
581 }
582 }
583
584 if(result!=result){ // nan
585 NuDEXException(__FILE__,std::to_string(__L
586 }
587
588 return result*ScaleFactor_M1;
589 }
590
591 G4double G4NuDEXPSF::GetE2(G4double Eg,G4doubl
592
593 G4double result=0;
594 for(G4int i=0;i<nR_E2;i++){
595 if(PSFType_E2[i]==0){
596 result+=5.22E-8*SLO(Eg,E_E2[i],G_E2[i],s
597 }
598 else if(PSFType_E2[i]==1){
599 result+=5.22E-8*EGLO(Eg,E_E2[i],G_E2[i],
600 }
601 else if(PSFType_E2[i]==2){
602 result+=5.22E-8*SMLO(Eg,E_E2[i],G_E2[i],
603 }
604 else if(PSFType_E2[i]==3){
605 result+=5.22E-8*GLO(Eg,E_E2[i],G_E2[i],s
606 }
607 else if(PSFType_E2[i]==4){
608 result+=5.22E-8*MGLO(Eg,E_E2[i],G_E2[i],
609 }
610 else if(PSFType_E2[i]==5){
611 result+=5.22E-8*KMF(Eg,E_E2[i],G_E2[i],s
612 }
613 else if(PSFType_E2[i]==6){
614 result+=5.22E-8*GH(Eg,E_E2[i],G_E2[i],s_
615 }
616 else if(PSFType_E2[i]==7){
617 result+=5.22E-8*MEGLO(Eg,E_E2[i],G_E2[i]
618 }
619 else if(PSFType_E2[i]==8){
620 result+=5.22E-8*MEGLO(Eg,E_E2[i],G_E2[i]
621 }
622 else if(PSFType_E2[i]==9){
623 result+=5.22E-8*MEGLO(Eg,E_E2[i],G_E2[i]
624 }
625 else if(PSFType_E2[i]==10){
626 result+=5.22E-8*MEGLO(Eg,E_E2[i],G_E2[i]
627 }
628 else if(PSFType_E2[i]==11){
629 result+=5.22E-8*SMLO_v2(Eg,E_E2[i],G_E2[
630 }
631 else if(PSFType_E2[i]==20){
632 result+=5.22E-8*Gauss(Eg,E_E2[i],G_E2[i]
633 }
634 else if(PSFType_E2[i]==21){
635 result+=5.22E-8*Expo(Eg,E_E2[i],G_E2[i])
636 }
637 else if(PSFType_E2[i]==40){
638 result+=EvaluateFunction(Eg,np_E2,x_E2,y
639 }
640 else if(PSFType_E2[i]==41){
641 result+=std::pow(10.,EvaluateFunction(Eg
642 }
643 else{
644 NuDEXException(__FILE__,std::to_string(_
645 }
646 }
647
648 if(result!=result){ // nan
649 NuDEXException(__FILE__,std::to_string(__L
650 }
651
652 return result*ScaleFactor_E2;
653 }
654
655
656 //********************************************
657 //********************************************
658 //********************************************
659
660 //Defined as in RIPL-3 when possible. Some of
661 //http://dx.doi.org/10.1103/PhysRevC.88.034317
662
663 G4double G4NuDEXPSF::SLO(G4double Eg,G4double
664
665 return sr*Gr*Eg*Gr/((Eg*Eg-Er*Er)*(Eg*Eg-Er*
666
667 }
668
669 //Kadmenskij-Markushev-Furman model (KMF) -->
670 G4double G4NuDEXPSF::KMF(G4double Eg,G4double
671
672 G4double Tf=0;
673 if(theLD!=0){
674 Tf=theLD->GetNucleusTemperature(Excitation
675 }
676 G4double Gc=Gr/Er/Er*(Eg*Eg+4*3.141592*3.141
677
678 if(Eg==Er){return 0;}
679
680 return 0.7*Er*Gr*sr*Gc/((Eg*Eg-Er*Er)*(Eg*Eg
681 }
682
683
684 G4double G4NuDEXPSF::EGLO(G4double Eg,G4double
685
686 G4double result=EGLO_GLO_MGLO(Eg,Er,Gr,sr,Ex
687
688 return result;
689 }
690
691 G4double G4NuDEXPSF::GLO(G4double Eg,G4double
692
693 G4double result=EGLO_GLO_MGLO(Eg,Er,Gr,sr,Ex
694
695 return result;
696 }
697
698 G4double G4NuDEXPSF::MGLO(G4double Eg,G4double
699
700 G4double result=EGLO_GLO_MGLO(Eg,Er,Gr,sr,Ex
701
702 return result;
703 }
704
705 //Hybrid model
706 G4double G4NuDEXPSF::GH(G4double Eg,G4double E
707
708 G4double Tf=0;
709 if(theLD!=0){
710 Tf=theLD->GetNucleusTemperature(Excitation
711 }
712
713 G4double Gamma_h=0.63*Gr/Eg/Er*(Eg*Eg+4*3.14
714
715 return sr*Gr*Eg*Gamma_h/((Eg*Eg-Er*Er)*(Eg*E
716
717 }
718
719 G4double G4NuDEXPSF::SMLO(G4double Eg,G4double
720
721 G4double Tf=0;
722 if(theLD!=0){
723 Tf=theLD->GetNucleusTemperature(Excitation
724 }
725
726 G4double Lambda=1/(1.-std::exp(-Eg/Tf));
727 G4double Gk_Eg=Gr/Er*ExcitationEnergy;
728
729 return Lambda*sr*Gr*Eg*Gk_Eg/((Eg*Eg-Er*Er)*
730 }
731
732
733 G4double G4NuDEXPSF::SMLO_v2(G4double Eg,G4dou
734
735 G4double Tf=0;
736 if(Eg<ExcitationEnergy){
737 Tf=std::sqrt((ExcitationEnergy-Eg)/(A_Int/
738 }
739
740 G4double Lambda=1/(1.-std::exp(-Eg/Tf));
741 G4double sig_trk=60.*(A_Int-Z_Int)*Z_Int/(G4
742 G4double Gk_Eg=Gr/Er*(Eg+4*3.141592*3.141592
743
744 return Lambda*sig_trk*2./3.141592*sr*Eg*Gk_E
745 }
746
747
748 G4double G4NuDEXPSF::Gauss(G4double Eg,G4doubl
749
750 return Area*(1./(sigma*std::sqrt(2.*3.141592
751
752 }
753
754 G4double G4NuDEXPSF::Expo(G4double Eg,G4double
755
756 return C*std::exp(-eta*Eg);
757
758 }
759
760
761 G4double G4NuDEXPSF::MEGLO(G4double Eg,G4doubl
762
763 G4double /*Ti=0,*/Tf=0;
764 if(Temp>=0){
765 //Ti=Temp;
766 Tf=Temp;
767 }
768 else if(theLD!=0){
769 //Ti=theLD->GetNucleusTemperature(Excitati
770 Tf=theLD->GetNucleusTemperature(Excitation
771 }
772
773 G4double Gk_Eg=Gamma_k(Eg,Er,Gr,Tf,k_param1)
774 //G4double Gk_0=Gamma_k(0,Er,Gr,Ti,k_param2)
775 G4double Gk_0=Gamma_k(0,Er,Gr,Tf,k_param2);
776
777 return sr*Gr*(Eg*Gk_Eg/((Eg*Eg-Er*Er)*(Eg*Eg
778
779 }
780
781
782
783
784 //Ti, Tf --> initial/final temperature of the
785 //Opt = 0,1,2 --> EGLO, GLO, MGLO
786 G4double G4NuDEXPSF::EGLO_GLO_MGLO(G4double Eg
787
788 G4double Ti=0,Tf=0;
789 if(theLD!=0){
790 Ti=theLD->GetNucleusTemperature(Excitation
791 Tf=theLD->GetNucleusTemperature(Excitation
792 }
793
794 //k_param could be modified according to exp
795 //The following expression is just a general
796 //If k_param==1 --> GLO
797 G4double k_param=1;
798 if(A_Int>=148){
799 k_param=1+0.09*(A_Int-148)*(A_Int-148)*std
800 }
801 G4double result=0;
802 if(Opt==0){//EGLO
803 result=FlexibleGLOType(Eg,Er,Gr,sr,Tf,k_pa
804 }
805 else if(Opt==1){//GLO --> same as EGLO, but
806 result=FlexibleGLOType(Eg,Er,Gr,sr,Tf,1,Ti
807 }
808 else if(Opt==2){//MGLO --> same as EGLO, but
809 result=FlexibleGLOType(Eg,Er,Gr,sr,Tf,k_pa
810 }
811 else{
812 NuDEXException(__FILE__,std::to_string(__L
813 }
814
815 return result;
816 }
817
818
819 G4double G4NuDEXPSF::FlexibleGLOType(G4double
820
821 G4double Gk_Eg=Gamma_k(Eg,Er,Gr,Temp1,k_para
822 //G4double Gk_0=Gamma_k(0,Er,Gr,Temp2,k_para
823 G4double Gk_0=Gamma_k(0,Er,Gr,Temp1,k_param2
824
825 return sr*Gr*(Eg*Gk_Eg/((Eg*Eg-Er*Er)*(Eg*Eg
826
827 }
828
829
830 G4double G4NuDEXPSF::Gamma_k(G4double Eg,G4dou
831
832 G4double eps0_param=4.5;
833 G4double Chi=1;
834 if(Er>eps0_param){
835 Chi=k_param+(1-k_param)*(Eg-eps0_param)/(E
836 }
837 G4double C_coll=Gr/Er/Er*Chi;
838 G4double Gamma_k=C_coll*(Eg*Eg+4*3.141592*3.
839
840 return Gamma_k;
841 }
842
843
844
845
846 //********************************************
847 //********************************************
848 //********************************************
849
850
851
852 void G4NuDEXPSF::PrintPSFParameters(std::ostre
853
854 out<<" #####################################
855 out<<" PSF_PARAMS"<<std::endl;
856 out<<" E1: nRes = "<<nR_E1<<std::endl;
857 for(G4int i=0;i<nR_E1;i++){
858 out<<" "<<PSFType_E1[i]<<" "<<E_E1[i]<<
859 if(PSFType_E1[i]==7){out<<"
860 if(PSFType_E1[i]==8){out<<"
861 if(PSFType_E1[i]==9){out<<"
862 if(PSFType_E1[i]==10){out<<"
863 if(PSFType_E1[i]==40 || PSFType_E1[i]==41)
864 }
865 out<<" M1: nRes = "<<nR_M1<<std::endl;
866 for(G4int i=0;i<nR_M1;i++){
867 out<<" "<<PSFType_M1[i]<<" "<<E_M1[i]<<
868 if(PSFType_M1[i]==7){out<<"
869 if(PSFType_M1[i]==8){out<<"
870 if(PSFType_M1[i]==9){out<<"
871 if(PSFType_M1[i]==10){out<<"
872 if(PSFType_M1[i]==40 || PSFType_M1[i]==41)
873 }
874 out<<" E2: nRes = "<<nR_E2<<std::endl;
875 for(G4int i=0;i<nR_E2;i++){
876 out<<" "<<PSFType_E2[i]<<" "<<E_E2[i]<<
877 if(PSFType_E2[i]==7){out<<"
878 if(PSFType_E2[i]==8){out<<"
879 if(PSFType_E2[i]==9){out<<"
880 if(PSFType_E2[i]==10){out<<"
881 if(PSFType_E2[i]==40 || PSFType_E2[i]==41)
882 }
883 out<<" #####################################
884
885 }
886
887
888 void G4NuDEXPSF::PrintPSFParametersInInputFile
889
890 out<<" PSF"<<std::endl;
891 G4long oldprc = out.precision(15);
892 out<<nR_E1<<std::endl;
893 for(G4int i=0;i<nR_E1;i++){
894 out<<" "<<PSFType_E1[i]<<" "<<E_E1[i]<<
895 if(PSFType_E1[i]==7){out<<" "<<p1_E1[i];
896 if(PSFType_E1[i]==8){out<<" "<<p1_E1[i]<
897 if(PSFType_E1[i]==9){out<<" "<<p1_E1[i]<
898 if(PSFType_E1[i]==10){out<<" "<<p1_E1[i]<
899 if(PSFType_E1[i]==40 || PSFType_E1[i]==41)
900 out<<std::endl;
901 }
902 out<<nR_M1<<std::endl;
903 for(G4int i=0;i<nR_M1;i++){
904 out<<" "<<PSFType_M1[i]<<" "<<E_M1[i]<<
905 if(PSFType_M1[i]==7){out<<"
906 if(PSFType_M1[i]==8){out<<"
907 if(PSFType_M1[i]==9){out<<"
908 if(PSFType_M1[i]==10){out<<"
909 if(PSFType_M1[i]==40 || PSFType_M1[i]==41)
910 out<<std::endl;
911 }
912 out<<nR_E2<<std::endl;
913 for(G4int i=0;i<nR_E2;i++){
914 out<<" "<<PSFType_E2[i]<<" "<<E_E2[i]<<
915 if(PSFType_E2[i]==7){out<<"
916 if(PSFType_E2[i]==8){out<<"
917 if(PSFType_E2[i]==9){out<<"
918 if(PSFType_E2[i]==10){out<<"
919 if(PSFType_E2[i]==40 || PSFType_E2[i]==41)
920 out<<std::endl;
921 }
922 out.precision(oldprc);
923 }
924
925 G4double G4NuDEXPSF::EvaluateFunction(G4double
926
927 if(xval<x[0]){return y[0];}
928 if(xval>x[np-1]){return y[np-1];}
929
930 G4double m,b;
931 G4int i_eval=np-1;
932 for(G4int i=1;i<np;i++){
933 if(x[i]>=xval){
934 i_eval=i;
935 break;
936 }
937 }
938
939 m=(y[i_eval]-y[i_eval-1])/(x[i_eval]-x[i_eva
940 b=y[i_eval]-m*x[i_eval];
941
942 return m*xval+b;
943 }
944
945