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 "G4NuDEXInternalConversion.hh" 43 44 45 46 //If alpha>0, use that value 47 G4bool G4NuDEXInternalConversion::SampleIntern 48 49 if(theZ<MINZINTABLES){ //then we have no inf 50 if(alpha<0){ 51 Ne=0; 52 Ng=0; 53 return false; 54 } 55 else{ 56 G4double rand=theRandom4->Uniform(0,alph 57 if(rand<alpha){ //then electron conversi 58 Ne=1; 59 Ng=0; 60 Eele[0]=Ene; //which is not correct, but we 61 return true; 62 } 63 return false; 64 } 65 } 66 67 68 Ne=0; 69 Ng=0; 70 71 if(multipolarity==0){ //maybe it is better t 72 //return true; 73 if(alpha<=0){ 74 return false; 75 } 76 //NuDEXException(__FILE__,std::to_string(_ 77 } 78 79 G4bool usegivenalpha=true; 80 if(NShells==0 || std::abs(multipolarity)>ICC 81 if(alpha<0){ 82 usegivenalpha=false; 83 alpha=GetICC(Ene,multipolarity); 84 } 85 86 G4double rand=theRandom4->Uniform(0,alpha+1) 87 if(rand<alpha){ //then electron conversion 88 if(!CalculateProducts){return true;} 89 //Select the orbital: 90 if(usegivenalpha){rand=rand*GetICC(Ene,mul 91 G4double cumul=0; 92 for(G4int i=1;i<NShells;i++){ 93 cumul+=GetICC(Ene,multipolarity,i); 94 //std::cout<<Ene<<" "<<multipolarity<<" 95 if(cumul>=rand || multipolarity==0){ //t 96 Ne=1; 97 Eele[0]=Ene-BindingEnergy[i]; 98 FillElectronHole(i); //now there is a hole t 99 if(Eele[0]<0){ 100 std::cout<<" For Z = "<<theZ<<" and orbita 101 std::cout<<" Given alpha is "<<alpha<<" (" 102 for(G4int j=1;j<=NShells;j++){ 103 std::cout<<j<<" "<<GetICC(Ene,multipola 104 } 105 Eele[0]=0; 106 } 107 return true; 108 } 109 } 110 std::cout<<" ############ Warning in "<<__ 111 std::cout<<" Given alpha is "<<alpha<<" an 112 for(G4int i=1;i<=NShells;i++){ 113 std::cout<<i<<" "<<GetICC(Ene,multipola 114 } 115 Ne=1; 116 Eele[0]=Ene-BindingEnergy[NShells-1]; 117 return true; 118 } 119 120 return false; 121 } 122 123 124 void G4NuDEXInternalConversion::FillElectronHo 125 126 //A very simplified version of the process ( 127 128 G4double fluoyield=0; 129 if(i_shell==1){ //K-shell 130 //Hubbell et al. (1994) formula for the fl 131 G4double C0=0.0370,C1=0.03112,C2=5.44e-5,C 132 G4double w_fac=std::pow(C0+C1*theZ+C2*theZ 133 fluoyield=w_fac/(1.+w_fac); 134 } 135 else if(i_shell>=2 && i_shell<=4){ //L-shell 136 //Hubbell et al. (1994) formula for the fl 137 if(theZ>=3 && theZ<=36){ 138 fluoyield=1.939e-8*std::pow(theZ,3.8874) 139 } 140 else if(theZ>36){ 141 G4double C0=0.17765,C1=0.00298937,C2=8.9 142 G4double w_fac=std::pow(C0+C1*theZ+C2*th 143 fluoyield=w_fac/(1.+w_fac); 144 } 145 } 146 147 148 G4double rand=theRandom4->Uniform(0,1); 149 if(rand<fluoyield){ //gamma emission 150 Egam[Ng]=BindingEnergy[i_shell]; 151 Ng++; 152 } 153 else{ //electron emission 154 Eele[Ne]=BindingEnergy[i_shell]; 155 Ne++; 156 } 157 158 159 } 160 161 162 163 164 165 166 167 //If i_shell<0 --> the total alpha 168 G4double G4NuDEXInternalConversion::GetICC(G4d 169 170 if(theZ<MINZINTABLES){ //then we have no inf 171 return 0; 172 } 173 174 if(NShells==0 || std::abs(multipolarity)>ICC 175 176 //----------------------------------------- 177 //Total: 178 //The following line does not work, due to i 179 //if(i_shell<0){i_shell=NShells;} 180 181 if(i_shell<0){ 182 G4double result=0; 183 for(G4int i=1;i<NShells;i++){ 184 result+=GetICC(Ene,multipolarity,i); 185 } 186 return result; 187 } 188 //----------------------------------------- 189 190 191 if(Ene<BindingEnergy[i_shell]){return 0;} 192 193 if(np[i_shell]==0){ 194 std::cout<<" shell "<<i_shell<<" has not b 195 NuDEXException(__FILE__,std::to_string(__L 196 } 197 198 if(i_shell==NShells && Ene<Eg[i_shell][0]){ 199 G4double total=0; 200 for(G4int i=1;i<NShells;i++){total+=GetICC 201 return total; 202 } 203 204 if(multipolarity>0){ 205 return Interpolate(Ene,np[i_shell],Eg[i_sh 206 } 207 else if(multipolarity<0){ 208 return Interpolate(Ene,np[i_shell],Eg[i_sh 209 } 210 return 0; 211 } 212 213 214 G4NuDEXInternalConversion::G4NuDEXInternalConv 215 theZ=Z; 216 NShells=0; 217 Ne=Ng=0; 218 for(G4int i=0;i<ICC_MAXNSHELLS;i++){ 219 Eg[i]=0; np[i]=0; BindingEnergy[i]=0; 220 for(G4int j=0;j<ICC_NMULTIP;j++){ 221 Icc_E[j][i]=0; Icc_M[j][i]=0; 222 } 223 } 224 theRandom4= new G4NuDEXRandom(1234567); 225 } 226 227 228 G4NuDEXInternalConversion::~G4NuDEXInternalCon 229 for(G4int i=0;i<ICC_MAXNSHELLS;i++){ 230 if(Eg[i]!=0){delete [] Eg[i];} 231 for(G4int j=0;j<ICC_NMULTIP;j++){ 232 if(Icc_E[j][i]!=0){delete [] Icc_E[j][i] 233 if(Icc_M[j][i]!=0){delete [] Icc_M[j][i] 234 } 235 } 236 delete theRandom4; 237 } 238 239 void G4NuDEXInternalConversion::PrintICC(std:: 240 241 char word[1000]; 242 out<<" ##################################### 243 out<<" ICC"<<std::endl; 244 out<<" Z = "<<theZ<<std::endl; 245 out<<" NShells = "<<NShells<<std::endl; 246 out<<" ------------------------------------- 247 out<<" Total calculated from the sum of the 248 out<<" E_g E1 E2 249 for(G4int j=0;j<np[NShells];j++){ 250 snprintf(word,1000,"%10.4g",Eg[NShells][j] 251 for(G4int k=0;k<ICC_NMULTIP;k++){ 252 snprintf(word,1000," %10.4g",Icc_E[k][N 253 } 254 for(G4int k=0;k<ICC_NMULTIP;k++){ 255 snprintf(word,1000," %10.4g",Icc_M[k][N 256 } 257 out<<std::endl; 258 } 259 out<<" ------------------------------------- 260 for(G4int i=0;i<NShells;i++){ 261 out<<" ----------------------------------- 262 out<<" Binding energy = "<<BindingEnergy[i 263 out<<" E_g E1 E2 264 for(G4int j=0;j<np[i];j++){ 265 snprintf(word,1000,"%10.4g",Eg[i][j]); o 266 for(G4int k=0;k<ICC_NMULTIP;k++){ 267 snprintf(word,1000," %10.4g",Icc_E[k][i][j] 268 } 269 for(G4int k=0;k<ICC_NMULTIP;k++){ 270 snprintf(word,1000," %10.4g",Icc_M[k][i][j] 271 } 272 out<<std::endl; 273 } 274 out<<" ----------------------------------- 275 } 276 out<<" ##################################### 277 } 278 279 void G4NuDEXInternalConversion::Init(const cha 280 281 if(theZ<MINZINTABLES){ //then we have no inf 282 return; 283 } 284 285 if(NShells!=0){ //Init only once 286 NuDEXException(__FILE__,std::to_string(__L 287 } 288 289 std::ifstream in(fname); 290 if(!in.good()){ 291 std::cout<<" ################ Error openin 292 NuDEXException(__FILE__,std::to_string(__L 293 } 294 std::string word; 295 NShells=1; 296 while(in>>word){ 297 if(word.c_str()[0]=='Z' && word.c_str()[1] 298 if(std::atoi(&(word.c_str()[2]))==theZ){ 299 in>>word>>word; 300 G4int orbindex; 301 if(word==std::string("Total")){ 302 in.ignore(1000,'\n'); 303 in.ignore(1000,'\n'); 304 orbindex=0; 305 } 306 else{ 307 orbindex=NShells; 308 in>>word>>word>>BindingEnergy[NShells]>>wo 309 BindingEnergy[NShells]*=1.e-6; // all in M 310 in.ignore(1000,'\n'); 311 in.ignore(1000,'\n'); 312 } 313 //------------------------------------------ 314 size_t sz,sz2; 315 G4int np_tmp=0; 316 G4double Eg_tmp[1000],Icc_E_tmp[ICC_NMULTIP] 317 while(getline(in,word)){ 318 if(word.size()<100){ 319 np[orbindex]=np_tmp; 320 Eg[orbindex]=new G4double[np_tmp]; 321 for(G4int j=0;j<np_tmp;j++){ 322 Eg[orbindex][j]=Eg_tmp[j]; 323 } 324 for(G4int i=0;i<ICC_NMULTIP;i++){ 325 Icc_E[i][orbindex]=new G4double[np_tmp 326 Icc_M[i][orbindex]=new G4double[np_tmp 327 } 328 for(G4int i=0;i<ICC_NMULTIP;i++){ 329 for(G4int j=0;j<np_tmp;j++){ 330 Icc_E[i][orbindex][j]=Icc_E_tmp[i][j]; 331 Icc_M[i][orbindex][j]=Icc_M_tmp[i][j]; 332 } 333 } 334 if(orbindex!=0){NShells++;} 335 break; 336 } 337 else{ 338 sz=0; 339 Eg_tmp[np_tmp]=std::stof(word,&sz2); 340 Eg_tmp[np_tmp]*=1.e-3; //all to MeV 341 sz+=sz2; 342 for(G4int i=0;i<ICC_NMULTIP;i++){ 343 Icc_E_tmp[i][np_tmp]=std::stof(word.su 344 } 345 for(G4int i=0;i<ICC_NMULTIP;i++){ 346 Icc_M_tmp[i][np_tmp]=std::stof(word.su 347 } 348 if((G4int)(std::stof(word.substr(sz),&sz 349 NuDEXException(__FILE__,std::to_string 350 } 351 sz+=sz2; 352 sz2=word.find_first_not_of(' ',sz); 353 if(np_tmp==0){OrbitalName[orbindex]=word 354 np_tmp++; 355 } 356 } 357 if(orbindex==0){break;} 358 //------------------------------------------ 359 } 360 } 361 } 362 if(!in.good()){ 363 NuDEXException(__FILE__,std::to_string(__L 364 } 365 in.close(); 366 367 MakeTotal(); 368 } 369 370 371 372 // Total Icc goes to nShell=NShells 373 void G4NuDEXInternalConversion::MakeTotal(){ 374 375 if(np[0]==0 || Eg[0]==0){ 376 NuDEXException(__FILE__,std::to_string(__L 377 } 378 379 //We evaluate it in the same frame as the to 380 BindingEnergy[NShells]=0; 381 np[NShells]=np[0]; 382 Eg[NShells]=new G4double[np[NShells]]; 383 for(G4int i=0;i<ICC_NMULTIP;i++){ 384 Icc_E[i][NShells]=new G4double[np[NShells] 385 Icc_M[i][NShells]=new G4double[np[NShells] 386 } 387 for(G4int k=0;k<np[NShells];k++){ 388 for(G4int j=0;j<ICC_NMULTIP;j++){ 389 Icc_E[j][NShells][k]=0; 390 Icc_M[j][NShells][k]=0; 391 } 392 } 393 394 395 for(G4int k=0;k<np[NShells];k++){ 396 Eg[NShells][k]=Eg[0][k]; 397 for(G4int i=1;i<NShells;i++){ 398 for(G4int j=0;j<ICC_NMULTIP;j++){ 399 Icc_E[j][NShells][k]+=GetICC(Eg[NShells][k], 400 Icc_M[j][NShells][k]+=GetICC(Eg[NShells][k], 401 } 402 } 403 } 404 405 } 406 407 //if val>x[npmax] then ---> return 0 408 G4double G4NuDEXInternalConversion::Interpolat 409 410 G4int i_interp=-1; 411 for(G4int i=1;i<npoints;i++){ 412 if(x[i]>=val){i_interp=i-1; break;} 413 } 414 if(i_interp<0){return 0;} 415 416 417 /* 418 //y=a0+a1*x 419 G4double a1=(y[i_interp+1]-y[i_interp])/(x[i 420 G4double a0=y[i_interp]-a1*x[i_interp]; 421 422 return (a0+a1*val); 423 */ 424 425 //It is better a log-log interpolation: 426 if(y[i_interp+1]<=0 || y[i_interp]<=0 || x[i 427 //y=a0+a1*x 428 G4double a1=(y[i_interp+1]-y[i_interp])/(x 429 G4double a0=y[i_interp]-a1*x[i_interp]; 430 431 return (a0+a1*val); 432 } 433 434 //log(y)=a0+a1*log(x) 435 G4double a1=std::log(y[i_interp+1]/y[i_inter 436 G4double a0=std::log(y[i_interp])-a1*std::lo 437 438 G4double result=std::exp(a0+a1*std::log(val) 439 return result; 440 } 441