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