| 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 "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 and 6 MeV
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 from the inputfname instead of the dirname
74 G4int G4NuDEXPSF::Init(const char* dirname,G4NuDEXLevelDensity* aLD,const char* inputfname,const char* defaultinputfname,G4int PSFflag){
75
76 theLD=aLD;
77
78 //Three options: very detailed model, if not --> gdr-parameters&errors-exp-MLO.dat (RIPL-3), if not --> theorethical values
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(defaultinputfname);
92 if(IsDone){return 0;}
93 }
94
95 //Detailed model
96 snprintf(fname,500,"%s/PSF/PSF_param.dat",dirname);
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_E1_v01.dat",dirname);
103 IsDone=TakePSFFromIAEA01(fname);
104 if(IsDone){return 0;}
105 }
106 else if(PSFflag!=1){
107 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
108 }
109
110 //RIPL-MLO values:
111 snprintf(fname,500,"%s/PSF/gdr-parameters&errors-exp-MLO.dat",dirname);
112 IsDone=TakePSFFromRIPL01(fname);
113 if(IsDone){return 0;}
114
115 //RIPL-Theorethical values:
116 snprintf(fname,500,"%s/PSF/gdr-parameters-theor.dat",dirname);
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; //SLO-old (RIPL-2)
125 //G4double a=27.47,b=22.063,c=0.0277,d=1.222;//SLO (RIPL-3)
126 G4double a=28.69,b=21.731,c=0.0285,d=1.267;//MLO (RIPL-3)
127
128 E_E1[nR_E1]=a*std::pow(A_Int,-1./3.)+b*std::pow(A_Int,-1./6.);
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_Int/(G4double)A_Int/G_E1[nR_E1];
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 MeV
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]/std::pow(A_Int,1./3)/G_E2[nR_E2];
158 PSFType_E2[nR_E2]=0;
159 nR_E2++;
160
161 }
162
163 G4bool G4NuDEXPSF::TakePSFFromRIPL02(const char* fname){
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.267;//MLO
180 G4double E_E1_0=a*std::pow(A_Int,-1./3.)+b*std::pow(A_Int,-1./6.);
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_Int)*Z_Int/(G4double)A_Int/G_E1_0;
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 char* fname){
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 resonance
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 char* fname){
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_E1[nR_E1]>>dum>>dum>>s_E1[nR_E1];
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("Er2")){
251 in>>dum>>E_E1[nR_E1]>>dum>>dum>>G_E1[nR_E1]>>dum>>dum>>s_E1[nR_E1]>>dum>>beta;
252 PSFType_E1[nR_E1]=11;
253 nR_E1++;
254
255 }
256 else{
257 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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-12840-1)
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_Int,9./10.);
281 nR_M1++;
282 //upbend:
283 PSFType_M1[nR_M1]=21;
284 E_M1[nR_M1]=0.4035*std::exp(-6.0*std::fabs(beta));
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]/std::pow(A_Int,1./3)/G_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 char* fname){
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")){break;}
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]>>p3_E1[i];}
322 if(PSFType_E1[i]==40 || PSFType_E1[i]==41){ //only one pointwise function is allowed
323 if(x_E1!=0){NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");}
324 in>>np_E1;
325 x_E1=new G4double[np_E1]; y_E1=new G4double[np_E1];
326 for(G4int j=0;j<np_E1;j++){in>>x_E1[j]>>y_E1[j];}
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]>>p3_M1[i];}
337 if(PSFType_M1[i]==40 || PSFType_M1[i]==41){//only one pointwise function is allowed
338 if(x_M1!=0){NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");}
339 in>>np_M1;
340 x_M1=new G4double[np_M1]; y_M1=new G4double[np_M1];
341 for(G4int j=0;j<np_M1;j++){in>>x_M1[j]>>y_M1[j];}
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]>>p3_E2[i];}
352 if(PSFType_E2[i]==40 || PSFType_E2[i]==41){//only one pointwise function is allowed
353 if(x_E2!=0){NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");}
354 in>>np_E2;
355 x_E2=new G4double[np_E2]; y_E2=new G4double[np_E2];
356 for(G4int j=0;j<np_E2;j++){in>>x_E2[j]>>y_E2[j];}
357 in>>E2_normFac;
358 }
359 }
360 break;
361 }
362 }
363
364 Renormalize(); // if XX_normFac>0 --> renormalization of the PSF
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)/npIntegral;
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<<" "<<ScaleFactor_M1<<std::endl; getchar();
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(const char* fname){
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]>>p3_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]>>p3_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]>>p3_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,G4double ExcitationEnergy){
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],s_E1[i]);
469 }
470 else if(PSFType_E1[i]==1){
471 result+=8.674E-8*EGLO(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy);
472 }
473 else if(PSFType_E1[i]==2){
474 result+=8.674E-8*SMLO(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy);
475 }
476 else if(PSFType_E1[i]==3){
477 result+=8.674E-8*GLO(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy);
478 }
479 else if(PSFType_E1[i]==4){
480 result+=8.674E-8*MGLO(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy);
481 }
482 else if(PSFType_E1[i]==5){
483 result+=8.674E-8*KMF(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy);
484 }
485 else if(PSFType_E1[i]==6){
486 result+=8.674E-8*GH(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy);
487 }
488 else if(PSFType_E1[i]==7){
489 result+=8.674E-8*MEGLO(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy,p1_E1[i],p1_E1[i]);
490 }
491 else if(PSFType_E1[i]==8){
492 result+=8.674E-8*MEGLO(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy,p1_E1[i],p2_E1[i]);
493 }
494 else if(PSFType_E1[i]==9){
495 result+=8.674E-8*MEGLO(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy,p1_E1[i],p1_E1[i],p2_E1[i]);
496 }
497 else if(PSFType_E1[i]==10){
498 result+=8.674E-8*MEGLO(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy,p1_E1[i],p2_E1[i],p3_E1[i]);
499 }
500 else if(PSFType_E1[i]==11){
501 result+=8.674E-8*SMLO_v2(Eg,E_E1[i],G_E1[i],s_E1[i],ExcitationEnergy);
502 }
503 else if(PSFType_E1[i]==20){
504 result+=8.674E-8*Gauss(Eg,E_E1[i],G_E1[i],s_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_E1);
511 }
512 else if(PSFType_E1[i]==41){
513 result+=std::pow(10.,EvaluateFunction(Eg,np_E1,x_E1,y_E1));
514 }
515 else{
516 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
517 }
518 }
519
520 if(result!=result){ // nan
521 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
522 }
523
524 return result*ScaleFactor_E1;
525 }
526
527 G4double G4NuDEXPSF::GetM1(G4double Eg,G4double ExcitationEnergy){
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],s_M1[i]);
533 }
534 else if(PSFType_M1[i]==1){
535 result+=8.674E-8*EGLO(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy);
536 }
537 else if(PSFType_M1[i]==2){
538 result+=8.674E-8*SMLO(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy);
539 }
540 else if(PSFType_M1[i]==3){
541 result+=8.674E-8*GLO(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy);
542 }
543 else if(PSFType_M1[i]==4){
544 result+=8.674E-8*MGLO(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy);
545 }
546 else if(PSFType_M1[i]==5){
547 result+=8.674E-8*KMF(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy);
548 }
549 else if(PSFType_M1[i]==6){
550 result+=8.674E-8*GH(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy);
551 }
552 else if(PSFType_M1[i]==7){
553 result+=8.674E-8*MEGLO(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy,p1_M1[i],p1_M1[i]);
554 }
555 else if(PSFType_M1[i]==8){
556 result+=8.674E-8*MEGLO(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy,p1_M1[i],p2_M1[i]);
557 }
558 else if(PSFType_M1[i]==9){
559 result+=8.674E-8*MEGLO(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy,p1_M1[i],p1_M1[i],p2_M1[i]);
560 }
561 else if(PSFType_M1[i]==10){
562 result+=8.674E-8*MEGLO(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy,p1_M1[i],p2_M1[i],p3_M1[i]);
563 }
564 else if(PSFType_M1[i]==11){
565 result+=8.674E-8*SMLO_v2(Eg,E_M1[i],G_M1[i],s_M1[i],ExcitationEnergy);
566 }
567 else if(PSFType_M1[i]==20){
568 result+=8.674E-8*Gauss(Eg,E_M1[i],G_M1[i],s_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_M1);
575 }
576 else if(PSFType_M1[i]==41){
577 result+=std::pow(10.,EvaluateFunction(Eg,np_M1,x_M1,y_M1));
578 }
579 else{
580 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
581 }
582 }
583
584 if(result!=result){ // nan
585 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
586 }
587
588 return result*ScaleFactor_M1;
589 }
590
591 G4double G4NuDEXPSF::GetE2(G4double Eg,G4double ExcitationEnergy){
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_E2[i]);
597 }
598 else if(PSFType_E2[i]==1){
599 result+=5.22E-8*EGLO(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy);
600 }
601 else if(PSFType_E2[i]==2){
602 result+=5.22E-8*SMLO(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy);
603 }
604 else if(PSFType_E2[i]==3){
605 result+=5.22E-8*GLO(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy);
606 }
607 else if(PSFType_E2[i]==4){
608 result+=5.22E-8*MGLO(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy);
609 }
610 else if(PSFType_E2[i]==5){
611 result+=5.22E-8*KMF(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy);
612 }
613 else if(PSFType_E2[i]==6){
614 result+=5.22E-8*GH(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy);
615 }
616 else if(PSFType_E2[i]==7){
617 result+=5.22E-8*MEGLO(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy,p1_E2[i],p1_E2[i]);
618 }
619 else if(PSFType_E2[i]==8){
620 result+=5.22E-8*MEGLO(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy,p1_E2[i],p2_E2[i]);
621 }
622 else if(PSFType_E2[i]==9){
623 result+=5.22E-8*MEGLO(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy,p1_E2[i],p1_E2[i],p2_E2[i]);
624 }
625 else if(PSFType_E2[i]==10){
626 result+=5.22E-8*MEGLO(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy,p1_E2[i],p2_E2[i],p3_E2[i]);
627 }
628 else if(PSFType_E2[i]==11){
629 result+=5.22E-8*SMLO_v2(Eg,E_E2[i],G_E2[i],s_E2[i],ExcitationEnergy);
630 }
631 else if(PSFType_E2[i]==20){
632 result+=5.22E-8*Gauss(Eg,E_E2[i],G_E2[i],s_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_E2);
639 }
640 else if(PSFType_E2[i]==41){
641 result+=std::pow(10.,EvaluateFunction(Eg,np_E2,x_E2,y_E2));
642 }
643 else{
644 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
645 }
646 }
647
648 if(result!=result){ // nan
649 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
650 }
651
652 return result*ScaleFactor_E2;
653 }
654
655
656 //**********************************************************************************************************
657 //**********************************************************************************************************
658 //**********************************************************************************************************
659
660 //Defined as in RIPL-3 when possible. Some of them come from other references:
661 //http://dx.doi.org/10.1103/PhysRevC.88.034317
662
663 G4double G4NuDEXPSF::SLO(G4double Eg,G4double Er,G4double Gr,G4double sr){
664
665 return sr*Gr*Eg*Gr/((Eg*Eg-Er*Er)*(Eg*Eg-Er*Er)+Eg*Eg*Gr*Gr);
666
667 }
668
669 //Kadmenskij-Markushev-Furman model (KMF) --> not well described in RIPL-3 manual, taken from another document
670 G4double G4NuDEXPSF::KMF(G4double Eg,G4double Er,G4double Gr,G4double sr,G4double ExcitationEnergy){
671
672 G4double Tf=0;
673 if(theLD!=0){
674 Tf=theLD->GetNucleusTemperature(ExcitationEnergy-Eg);
675 }
676 G4double Gc=Gr/Er/Er*(Eg*Eg+4*3.141592*3.141592*Tf*Tf);
677
678 if(Eg==Er){return 0;}
679
680 return 0.7*Er*Gr*sr*Gc/((Eg*Eg-Er*Er)*(Eg*Eg-Er*Er));
681 }
682
683
684 G4double G4NuDEXPSF::EGLO(G4double Eg,G4double Er,G4double Gr,G4double sr,G4double ExcitationEnergy){
685
686 G4double result=EGLO_GLO_MGLO(Eg,Er,Gr,sr,ExcitationEnergy,0);
687
688 return result;
689 }
690
691 G4double G4NuDEXPSF::GLO(G4double Eg,G4double Er,G4double Gr,G4double sr,G4double ExcitationEnergy){
692
693 G4double result=EGLO_GLO_MGLO(Eg,Er,Gr,sr,ExcitationEnergy,1);
694
695 return result;
696 }
697
698 G4double G4NuDEXPSF::MGLO(G4double Eg,G4double Er,G4double Gr,G4double sr,G4double ExcitationEnergy){
699
700 G4double result=EGLO_GLO_MGLO(Eg,Er,Gr,sr,ExcitationEnergy,2);
701
702 return result;
703 }
704
705 //Hybrid model
706 G4double G4NuDEXPSF::GH(G4double Eg,G4double Er,G4double Gr,G4double sr,G4double ExcitationEnergy){
707
708 G4double Tf=0;
709 if(theLD!=0){
710 Tf=theLD->GetNucleusTemperature(ExcitationEnergy-Eg);
711 }
712
713 G4double Gamma_h=0.63*Gr/Eg/Er*(Eg*Eg+4*3.141592*3.141592*Tf*Tf);
714
715 return sr*Gr*Eg*Gamma_h/((Eg*Eg-Er*Er)*(Eg*Eg-Er*Er)+Eg*Eg*Gr*Gamma_h);
716
717 }
718
719 G4double G4NuDEXPSF::SMLO(G4double Eg,G4double Er,G4double Gr,G4double sr,G4double ExcitationEnergy){
720
721 G4double Tf=0;
722 if(theLD!=0){
723 Tf=theLD->GetNucleusTemperature(ExcitationEnergy-Eg);
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)*(Eg*Eg-Er*Er)+Eg*Eg*Gk_Eg*Gk_Eg);
730 }
731
732
733 G4double G4NuDEXPSF::SMLO_v2(G4double Eg,G4double Er,G4double Gr,G4double sr,G4double ExcitationEnergy){
734
735 G4double Tf=0;
736 if(Eg<ExcitationEnergy){
737 Tf=std::sqrt((ExcitationEnergy-Eg)/(A_Int/10.));
738 }
739
740 G4double Lambda=1/(1.-std::exp(-Eg/Tf));
741 G4double sig_trk=60.*(A_Int-Z_Int)*Z_Int/(G4double)A_Int;
742 G4double Gk_Eg=Gr/Er*(Eg+4*3.141592*3.141592*Tf*Tf/Er);
743
744 return Lambda*sig_trk*2./3.141592*sr*Eg*Gk_Eg/((Eg*Eg-Er*Er)*(Eg*Eg-Er*Er)+Eg*Eg*Gk_Eg*Gk_Eg);
745 }
746
747
748 G4double G4NuDEXPSF::Gauss(G4double Eg,G4double Er,G4double sigma,G4double Area){
749
750 return Area*(1./(sigma*std::sqrt(2.*3.141592)))*std::exp(-0.5*std::pow((Eg-Er)/sigma,2.));
751
752 }
753
754 G4double G4NuDEXPSF::Expo(G4double Eg,G4double C,G4double eta){
755
756 return C*std::exp(-eta*Eg);
757
758 }
759
760
761 G4double G4NuDEXPSF::MEGLO(G4double Eg,G4double Er,G4double Gr,G4double sr,G4double ExcitationEnergy,G4double k_param1,G4double k_param2,G4double Temp){
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(ExcitationEnergy);
770 Tf=theLD->GetNucleusTemperature(ExcitationEnergy-Eg);
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); // in most of the references they use just one temperature
776
777 return sr*Gr*(Eg*Gk_Eg/((Eg*Eg-Er*Er)*(Eg*Eg-Er*Er)+Eg*Eg*Gk_Eg*Gk_Eg)+0.7*Gk_0/Er/Er/Er);
778
779 }
780
781
782
783
784 //Ti, Tf --> initial/final temperature of the nucleus
785 //Opt = 0,1,2 --> EGLO, GLO, MGLO
786 G4double G4NuDEXPSF::EGLO_GLO_MGLO(G4double Eg,G4double Er,G4double Gr,G4double sr,G4double ExcitationEnergy,G4int Opt){
787
788 G4double Ti=0,Tf=0;
789 if(theLD!=0){
790 Ti=theLD->GetNucleusTemperature(ExcitationEnergy);
791 Tf=theLD->GetNucleusTemperature(ExcitationEnergy-Eg);
792 }
793
794 //k_param could be modified according to experimental data.
795 //The following expression is just a general recomendation
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::exp(-0.18*(A_Int-148));
800 }
801 G4double result=0;
802 if(Opt==0){//EGLO
803 result=FlexibleGLOType(Eg,Er,Gr,sr,Tf,k_param,Ti,k_param);
804 }
805 else if(Opt==1){//GLO --> same as EGLO, but k_param=1
806 result=FlexibleGLOType(Eg,Er,Gr,sr,Tf,1,Ti,1);
807 }
808 else if(Opt==2){//MGLO --> same as EGLO, but k_param2=1
809 result=FlexibleGLOType(Eg,Er,Gr,sr,Tf,k_param,Ti,1);
810 }
811 else{
812 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
813 }
814
815 return result;
816 }
817
818
819 G4double G4NuDEXPSF::FlexibleGLOType(G4double Eg,G4double Er,G4double Gr,G4double sr,G4double Temp1,G4double k_param1,G4double /*Temp2*/,G4double k_param2){
820
821 G4double Gk_Eg=Gamma_k(Eg,Er,Gr,Temp1,k_param1);
822 //G4double Gk_0=Gamma_k(0,Er,Gr,Temp2,k_param2);
823 G4double Gk_0=Gamma_k(0,Er,Gr,Temp1,k_param2); // in most of the references they use just one temperature
824
825 return sr*Gr*(Eg*Gk_Eg/((Eg*Eg-Er*Er)*(Eg*Eg-Er*Er)+Eg*Eg*Gk_Eg*Gk_Eg)+0.7*Gk_0/Er/Er/Er);
826
827 }
828
829
830 G4double G4NuDEXPSF::Gamma_k(G4double Eg,G4double Er,G4double Gr,G4double Temp,G4double k_param){
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)/(Er-eps0_param);
836 }
837 G4double C_coll=Gr/Er/Er*Chi;
838 G4double Gamma_k=C_coll*(Eg*Eg+4*3.141592*3.141592*Temp*Temp);
839
840 return Gamma_k;
841 }
842
843
844
845
846 //**********************************************************************************************************
847 //**********************************************************************************************************
848 //**********************************************************************************************************
849
850
851
852 void G4NuDEXPSF::PrintPSFParameters(std::ostream &out){
853
854 out<<" ###################################################################################### "<<std::endl;
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]<<" "<<G_E1[i]<<" "<<s_E1[i]<<std::endl;
859 if(PSFType_E1[i]==7){out<<" "<<p1_E1[i]<<std::endl;}
860 if(PSFType_E1[i]==8){out<<" "<<p1_E1[i]<<" "<<p2_E1[i]<<std::endl;}
861 if(PSFType_E1[i]==9){out<<" "<<p1_E1[i]<<" "<<p2_E1[i]<<std::endl;}
862 if(PSFType_E1[i]==10){out<<" "<<p1_E1[i]<<" "<<p2_E1[i]<<" "<<p3_E1[i]<<std::endl;}
863 if(PSFType_E1[i]==40 || PSFType_E1[i]==41){out<<np_E1; for(G4int j=0;j<np_E1;j++){out<<" "<<x_E1[j]<<" "<<y_E1[j];} out<<std::endl;}
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]<<" "<<G_M1[i]<<" "<<s_M1[i]<<std::endl;
868 if(PSFType_M1[i]==7){out<<" "<<p1_M1[i]<<std::endl;}
869 if(PSFType_M1[i]==8){out<<" "<<p1_M1[i]<<" "<<p2_M1[i]<<std::endl;}
870 if(PSFType_M1[i]==9){out<<" "<<p1_M1[i]<<" "<<p2_M1[i]<<std::endl;}
871 if(PSFType_M1[i]==10){out<<" "<<p1_M1[i]<<" "<<p2_M1[i]<<" "<<p3_M1[i]<<std::endl;}
872 if(PSFType_M1[i]==40 || PSFType_M1[i]==41){out<<np_M1; for(G4int j=0;j<np_M1;j++){out<<" "<<x_M1[j]<<" "<<y_M1[j];} out<<std::endl;}
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]<<" "<<G_E2[i]<<" "<<s_E2[i]<<std::endl;
877 if(PSFType_E2[i]==7){out<<" "<<p1_E2[i]<<std::endl;}
878 if(PSFType_E2[i]==8){out<<" "<<p1_E2[i]<<" "<<p2_E2[i]<<std::endl;}
879 if(PSFType_E2[i]==9){out<<" "<<p1_E2[i]<<" "<<p2_E2[i]<<std::endl;}
880 if(PSFType_E2[i]==10){out<<" "<<p1_E2[i]<<" "<<p2_E2[i]<<" "<<p3_E2[i]<<std::endl;}
881 if(PSFType_E2[i]==40 || PSFType_E2[i]==41){out<<np_E2; for(G4int j=0;j<np_E2;j++){out<<" "<<x_E2[j]<<" "<<y_E2[j];} out<<std::endl;}
882 }
883 out<<" ###################################################################################### "<<std::endl;
884
885 }
886
887
888 void G4NuDEXPSF::PrintPSFParametersInInputFileFormat(std::ostream &out){
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]<<" "<<G_E1[i]<<" "<<s_E1[i];
895 if(PSFType_E1[i]==7){out<<" "<<p1_E1[i];}
896 if(PSFType_E1[i]==8){out<<" "<<p1_E1[i]<<" "<<p2_E1[i];}
897 if(PSFType_E1[i]==9){out<<" "<<p1_E1[i]<<" "<<p2_E1[i];}
898 if(PSFType_E1[i]==10){out<<" "<<p1_E1[i]<<" "<<p2_E1[i]<<" "<<p3_E1[i];}
899 if(PSFType_E1[i]==40 || PSFType_E1[i]==41){out<<np_E1; for(G4int j=0;j<np_E1;j++){out<<" "<<x_E1[j]<<" "<<y_E1[j];} }
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]<<" "<<G_M1[i]<<" "<<s_M1[i];
905 if(PSFType_M1[i]==7){out<<" "<<p1_M1[i];}
906 if(PSFType_M1[i]==8){out<<" "<<p1_M1[i]<<" "<<p2_M1[i];}
907 if(PSFType_M1[i]==9){out<<" "<<p1_M1[i]<<" "<<p2_M1[i];}
908 if(PSFType_M1[i]==10){out<<" "<<p1_M1[i]<<" "<<p2_M1[i]<<" "<<p3_M1[i];}
909 if(PSFType_M1[i]==40 || PSFType_M1[i]==41){out<<np_M1; for(G4int j=0;j<np_M1;j++){out<<" "<<x_M1[j]<<" "<<y_M1[j];}}
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]<<" "<<G_E2[i]<<" "<<s_E2[i];
915 if(PSFType_E2[i]==7){out<<" "<<p1_E2[i];}
916 if(PSFType_E2[i]==8){out<<" "<<p1_E2[i]<<" "<<p2_E2[i];}
917 if(PSFType_E2[i]==9){out<<" "<<p1_E2[i]<<" "<<p2_E2[i];}
918 if(PSFType_E2[i]==10){out<<" "<<p1_E2[i]<<" "<<p2_E2[i]<<" "<<p3_E2[i];}
919 if(PSFType_E2[i]==40 || PSFType_E2[i]==41){out<<np_E2; for(G4int j=0;j<np_E2;j++){out<<" "<<x_E2[j]<<" "<<y_E2[j];}}
920 out<<std::endl;
921 }
922 out.precision(oldprc);
923 }
924
925 G4double G4NuDEXPSF::EvaluateFunction(G4double xval,G4int np,G4double* x,G4double* y){
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_eval-1]);
940 b=y[i_eval]-m*x[i_eval];
941
942 return m*xval+b;
943 }
944
945