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 "G4NuDEXStatisticalNucleus.hh" 44 #include "G4NuDEXLevelDensity.hh" 45 #include "G4NuDEXInternalConversion.hh" 46 #include "G4NuDEXPSF.hh" 47 48 49 50 G4NuDEXStatisticalNucleus::G4NuDEXStatisticalN 51 52 //The default values for these flags are in 53 //Can be changed with G4NuDEXStatisticalNucl 54 LevelDensityType=-1; 55 PSFflag=-1; 56 maxspinx2=-1; 57 MinLevelsPerBand=-1; 58 BandWidth=0; 59 MaxExcEnergy=0; 60 BROpt=-1; 61 SampleGammaWidths=-1; 62 63 //The default values for these flags are in 64 //Can be changed via G4NuDEXStatisticalNucle 65 ElectronConversionFlag=-1; // All EC 66 KnownLevelsFlag=-1; //Use all known levels 67 PrimaryGammasIntensityNormFactor=-1; 68 PrimaryGammasEcut=-1; //If primary gammas, d 69 Ecrit=-1; 70 71 hasBeenInitialized=false; 72 NBands=-1; 73 theLevels=0; 74 theKnownLevels=0; 75 NKnownLevels=0; NUnknownLevels=0; NLevels=0; 76 theRandom1=0; 77 theRandom2=0; 78 theRandom3=0; 79 theLD=0; 80 theICC=0; 81 thePSF=0; 82 TotalGammaRho=0; 83 theThermalCaptureLevelCumulBR=0; 84 TotalCumulBR=0; 85 86 Z_Int=Z; 87 A_Int=A; 88 89 //Random generators: 90 seed1=1234567; 91 seed2=1234567; 92 seed3=1234567; 93 theRandom1= new G4NuDEXRandom(seed1); 94 theRandom2= new G4NuDEXRandom(seed2); 95 theRandom3= new G4NuDEXRandom(seed3); 96 Rand1seedProvided=false; Rand2seedProvided=f 97 } 98 99 void G4NuDEXStatisticalNucleus::SetInitialPara 100 if(hasBeenInitialized){ 101 //std::cout<<" ############## Error: G4NuD 102 NuDEXException(__FILE__,std::to_string(__L 103 } 104 if(knownLevelsFlag>=0){KnownLevelsFlag=known 105 if(electronConversionFlag>=0){ElectronConver 106 if(primGamNormFactor>=0){PrimaryGammasIntens 107 if(primGamEcut>=0){PrimaryGammasEcut=primGam 108 if(ecrit>=0){Ecrit=ecrit;} 109 110 } 111 112 void G4NuDEXStatisticalNucleus::SetSomeInitalP 113 114 if(hasBeenInitialized){ 115 std::cout<<" ############## Error: G4NuDEX 116 NuDEXException(__FILE__,std::to_string(__L 117 } 118 119 if(LDtype>0){LevelDensityType=LDtype;} 120 if(PSFFlag>=0){PSFflag=PSFFlag;} 121 if(MaxSpin>0){maxspinx2=(G4int)(2.*MaxSpin+0 122 if(minlevelsperband>0){MinLevelsPerBand=minl 123 if(BandWidth_MeV!=0){BandWidth=BandWidth_MeV 124 if(maxExcEnergy!=0){MaxExcEnergy=maxExcEnerg 125 if(BrOption>0){BROpt=BrOption;} 126 if(sampleGammaWidths>=0){SampleGammaWidths=s 127 if(aseed1>0){seed1=aseed1; theRandom1->SetSe 128 if(aseed2>0){seed2=aseed2; theRandom2->SetSe 129 if(aseed3>0){seed3=aseed3; theRandom3->SetSe 130 131 } 132 133 134 135 G4NuDEXStatisticalNucleus::~G4NuDEXStatistical 136 137 if(theLevels!=0){delete [] theLevels;} 138 for(G4int i=0;i<KnownLevelsVectorSize;i++){ 139 if(theKnownLevels[i].Ndecays>0){ 140 delete [] theKnownLevels[i].decayFractio 141 delete [] theKnownLevels[i].decayMode; 142 } 143 if(theKnownLevels[i].NGammas>0){ 144 delete [] theKnownLevels[i].FinalLevelID 145 delete [] theKnownLevels[i].multipolarit 146 delete [] theKnownLevels[i].Eg; 147 delete [] theKnownLevels[i].cumulPtot; 148 delete [] theKnownLevels[i].Pg; 149 delete [] theKnownLevels[i].Pe; 150 delete [] theKnownLevels[i].Icc; 151 } 152 } 153 if(theKnownLevels!=0){delete [] theKnownLeve 154 if(theRandom1!=0){delete theRandom1;} 155 if(theRandom2!=0){delete theRandom2;} 156 if(theRandom3!=0){delete theRandom3;} 157 if(theLD!=0){delete theLD;} 158 if(theICC!=0){delete theICC;} 159 if(thePSF!=0){delete thePSF;} 160 if(TotalGammaRho!=0){delete [] TotalGammaRho 161 if(theThermalCaptureLevelCumulBR!=0){delete 162 if(TotalCumulBR!=0){ 163 for(G4int i=0;i<NLevels;i++){ 164 if(TotalCumulBR[i]!=0){delete [] TotalCu 165 } 166 delete [] TotalCumulBR; 167 } 168 } 169 170 171 G4int G4NuDEXStatisticalNucleus::Init(const ch 172 173 hasBeenInitialized=true; 174 //------------------------------------------ 175 //First, we read data from files: 176 G4int check=0; 177 char fname[1000],defaultinputfname[1000]; 178 theLibDir=std::string(dirname); 179 180 //Special (default) input file: 181 snprintf(defaultinputfname,1000,"%s/SpecialI 182 G4int HasDefaultInput=ReadSpecialInputFile(d 183 char* definputfn=0; 184 if(HasDefaultInput>0){definputfn=defaultinpu 185 186 //General statistical parameters: 187 snprintf(fname,1000,"%s/GeneralStatNuclParam 188 check=ReadGeneralStatNuclParameters(fname); 189 190 //Some default, if not initialized yet: 191 if(ElectronConversionFlag<0){ElectronConvers 192 if(KnownLevelsFlag<0){KnownLevelsFlag=1;} // 193 if(PrimaryGammasIntensityNormFactor<0){Prima 194 if(PrimaryGammasEcut<0){PrimaryGammasEcut=0; 195 if(Ecrit<0){ 196 snprintf(fname,1000,"%s/KnownLevels/levels 197 check=ReadEcrit(fname); if(check<0){return 198 } 199 200 201 //Level density: 202 theLD=new G4NuDEXLevelDensity(Z_Int,A_Int,Le 203 check=theLD->ReadLDParameters(dirname,inputf 204 LevelDensityType=theLD->GetLDType(); //becau 205 if(check<0){ 206 delete theLD; theLD=0; 207 Sn=-1; D0=-1; I0=-1000; 208 } 209 else{ 210 theLD->GetSnD0I0Vals(Sn,D0,I0); 211 } 212 213 //Known level sheme: 214 snprintf(fname,1000,"%s/KnownLevels/z%03d.da 215 check=ReadKnownLevels(fname); if(check<0){re 216 I0=TakeTargetNucleiI0(fname,check); if(check 217 218 if(MaxExcEnergy<=0){ 219 if(Sn>0){ 220 MaxExcEnergy=Sn-MaxExcEnergy; 221 } 222 else{ 223 MaxExcEnergy=1-MaxExcEnergy; 224 } 225 } 226 227 //If we don't have level density and the kno 228 if(theLD==0 && Ecrit<MaxExcEnergy){ 229 std::cout<<" ###### WARNING: No level dens 230 return -1; 231 } 232 //------------------------------------------ 233 234 //------------------------------------------ 235 //Init some variables: 236 E_unk_min=Ecrit; 237 E_unk_max=MaxExcEnergy; 238 239 NBands=0; 240 if(BandWidth>0){//then we have to create som 241 NBands=0; 242 while(E_unk_min+BandWidth*NBands<MaxExcEne 243 NBands++; 244 } 245 E_unk_max=E_unk_min+BandWidth*NBands; 246 } 247 248 Emin_bands=E_unk_min; 249 Emax_bands=E_unk_max; 250 //------------------------------------------ 251 252 253 //Make some checks: 254 MakeSomeParameterChecks01(); 255 256 //Level scheme: 257 //std::cout<<" creating level scheme ..."<<s 258 CreateLevelScheme(); 259 //std::cout<<" ............. done"<<std::end 260 261 if(KnownLevelsFlag==1){ 262 InsertHighEnergyKnownLevels(); 263 } 264 265 //We set the precursors for each level: 266 for(G4int i=0;i<NLevels;i++){ 267 theLevels[NLevels-1-i].seed=theRandom2->In 268 } 269 270 //Internal conversion: 271 theICC=new G4NuDEXInternalConversion(Z_Int); 272 snprintf(fname,1000,"%s/ICC_factors.dat",dir 273 theICC->Init(fname); 274 theICC->SetRandom4Seed(theRandom3->GetSeed() 275 276 //PSF: 277 thePSF=new G4NuDEXPSF(Z_Int,A_Int); 278 thePSF->Init(dirname,theLD,inputfname,definp 279 280 //We compute the missing BR in the known par 281 ComputeKnownLevelsMissingBR(); 282 283 //Init TotalGammaRho: 284 TotalGammaRho=new G4double[NLevels]; 285 for(G4int i=0;i<NLevels-1;i++){ 286 TotalGammaRho[i]=-1; 287 } 288 289 //Thermal capture level: 290 if(Sn>0 && NLevels>1){ 291 CreateThermalCaptureLevel(); 292 GenerateThermalCaptureLevelBR(dirname); 293 } 294 295 //Init TotalCumulBR, if BROpt==1,2 296 if(BROpt==1 || BROpt==2){ 297 TotalCumulBR=new G4double*[NLevels]; 298 for(G4int i=0;i<NLevels;i++){ 299 TotalCumulBR[i]=0; 300 } 301 } 302 303 return 0; 304 } 305 306 307 void G4NuDEXStatisticalNucleus::MakeSomeParame 308 309 if(LevelDensityType<1 || LevelDensityType>3) 310 std::cout<<" ############## Error, LevelDe 311 } 312 313 if(maxspinx2<=0){ 314 std::cout<<" ############## Error, MaxSpin 315 } 316 317 if(MaxExcEnergy<=0){ 318 std::cout<<" ############## Error, MaxExcE 319 } 320 321 if(BROpt<0 || BROpt>2){ 322 std::cout<<" ############## Error, BROpt c 323 } 324 325 if(SampleGammaWidths<0 || SampleGammaWidths> 326 std::cout<<" ############## Error, SampleG 327 } 328 329 if(KnownLevelsFlag!=0 && KnownLevelsFlag!=1) 330 std::cout<<" ############## Error, KnownLe 331 } 332 333 if(ElectronConversionFlag!=0 && ElectronConv 334 std::cout<<" ############## Error, Electro 335 } 336 337 if(PrimaryGammasIntensityNormFactor<=0){ 338 std::cout<<" ############## Error, Primary 339 } 340 341 if(PrimaryGammasEcut<0){ 342 std::cout<<" ############## Error, Primary 343 } 344 345 if(Ecrit<0){ 346 std::cout<<" ############## Error, Ecrit c 347 } 348 349 } 350 351 352 //If InitialLevel==-1 then we start from the t 353 //If ExcitationEnergy>0 then is the excitation 354 //If ExcitationEnergy<0 then is a capture reac 355 // return Npar (number of particles emitted). 356 G4int G4NuDEXStatisticalNucleus::GenerateCasca 357 358 pType.clear(); 359 pEnergy.clear(); 360 pTime.clear(); 361 362 if(ExcitationEnergy<0){ 363 ExcitationEnergy=Sn-(A_Int-1.)/(G4double)A 364 } 365 if(ExcitationEnergy<=0){ 366 return 0; 367 } 368 369 G4int Npar=0; 370 G4int f_level=0,multipol=0; 371 G4double alpha,E_trans,Exc_ene_i,Exc_ene_f; 372 G4double EmissionTime=0; //in seconds 373 G4int NTransition=0; 374 //G4double TotalCascadeEnergy1=0,TotalCascad 375 376 //Start: 377 G4int i_level=InitialLevel; 378 Exc_ene_i=ExcitationEnergy; 379 380 381 if(i_level==0){ //could happen 382 pType.push_back('g'); 383 pEnergy.push_back(Exc_ene_i); 384 pTime.push_back(0); 385 Npar++; 386 } 387 388 //Loop: 389 while(i_level!=0){ 390 391 NTransition++; 392 //---------------------------------------- 393 //Sample final level: 394 if(i_level==-1){ //thermal level 395 if(!theThermalCaptureLevelCumulBR){ 396 f_level=0; 397 std::cout<<" ############## NuDEX: WARNING, 398 } 399 else{ 400 //Sample final level: 401 G4double randnumber=theRandom3->Uniform(); 402 f_level=-1; 403 for(G4int i=0;i<NLevelsBelowThermalCaptureLe 404 if(theThermalCaptureLevelCumulBR[i]>randnu 405 multipol=GetMultipolarity(&theThermalCap 406 f_level=i; 407 break; 408 } 409 } 410 } 411 if(f_level<0){ 412 NuDEXException(__FILE__,std::to_string(__LIN 413 } 414 Exc_ene_f=theLevels[f_level].Energy; 415 } 416 else if(i_level>0){ 417 f_level=SampleFinalLevel(i_level,multipo 418 } 419 else{ 420 NuDEXException(__FILE__,std::to_string(_ 421 } 422 //---------------------------------------- 423 424 //Energy of the transition: 425 Exc_ene_f=theLevels[f_level].Energy; 426 427 //We sample the final energy if it is a ba 428 if(theLevels[f_level].Width!=0){ 429 Exc_ene_f+=theRandom3->Uniform(-theLevel 430 } 431 E_trans=Exc_ene_i-Exc_ene_f; 432 if(E_trans<=0){ 433 //std::cout<<"Exc_ene_i = "<<Exc_ene_i<< 434 //std::cout<<" ####### WARNING: E_trans 435 return -1; 436 } 437 //---------------------------------------- 438 //Emission time: 439 if(i_level<NKnownLevels && i_level>0){ 440 if(theKnownLevels[i_level].T12>0){ 441 EmissionTime+=theRandom3->Exp(theKnownLevels 442 } 443 } 444 //---------------------------------------- 445 446 //---------------------------------------- 447 //calculate electron conversion: 448 G4bool ele_conv=false; 449 if(ElectronConversionFlag>0){ 450 if(i_level<NKnownLevels && i_level>0){ / 451 ele_conv=theICC->SampleInternalConversion(E_ 452 } 453 else if(ElectronConversionFlag==2){ 454 ele_conv=theICC->SampleInternalConvers 455 } 456 } 457 //---------------------------------------- 458 //std::cout<<" ---- "<<Exc_ene_i<<" "<<Ex 459 //TotalCascadeEnergy1+=E_trans; 460 461 //---------------------------------------- 462 //Fill result: 463 if(ele_conv){ 464 for(G4int i=0;i<theICC->Ne;i++){ 465 pType.push_back('e'); 466 pEnergy.push_back(theICC->Eele[i]); 467 pTime.push_back(EmissionTime); 468 Npar++; 469 } 470 for(G4int i=0;i<theICC->Ng;i++){ 471 pType.push_back('g'); 472 pEnergy.push_back(theICC->Egam[i]); 473 pTime.push_back(EmissionTime); 474 Npar++; 475 } 476 } 477 else{ 478 pType.push_back('g'); 479 pEnergy.push_back(E_trans); 480 pTime.push_back(EmissionTime); 481 Npar++; 482 } 483 //---------------------------------------- 484 i_level=f_level; 485 Exc_ene_i=Exc_ene_f; 486 } 487 488 //for(G4int i=0;i<Npar;i++){TotalCascadeEner 489 //std::cout<<" Total energy: "<<TotalCascade 490 491 492 if(i_level!=0){ 493 NuDEXException(__FILE__,std::to_string(__L 494 } 495 return Npar; 496 } 497 498 499 500 501 G4int G4NuDEXStatisticalNucleus::SampleFinalLe 502 503 if(i_level<=0 || i_level>=NLevels){ 504 NuDEXException(__FILE__,std::to_string(__L 505 } 506 507 G4double randnumber=theRandom3->Uniform(); 508 509 G4int i_knownLevel=-1; 510 if(i_level<NKnownLevels){ //then is a known 511 i_knownLevel=i_level; 512 } 513 if(theLevels[i_level].KnownLevelID>0){ //the 514 if(theKnownLevels[theLevels[i_level].Known 515 i_knownLevel=theLevels[i_level].KnownLev 516 } 517 } 518 519 if(i_knownLevel>=0){//known part of the spec 520 theSampledLevel=-1; 521 for(G4int j=0;j<theKnownLevels[i_knownLeve 522 if(theKnownLevels[i_knownLevel].cumulPto 523 multipolarity=theKnownLevels[i_knownLevel].m 524 icc_fac=theKnownLevels[i_knownLevel].Icc[j]; 525 return theKnownLevels[i_knownLevel].FinalLev 526 } 527 } 528 std::cout<<randnumber<<" "<<i_knownLevel< 529 for(G4int j=0;j<theKnownLevels[i_knownLeve 530 std::cout<<theKnownLevels[i_knownLevel]. 531 } 532 NuDEXException(__FILE__,std::to_string(__L 533 } 534 else{ 535 icc_fac=-1; 536 //---------------------------------------- 537 //If BROpt==1 or 2, then we store the BR, 538 if(BROpt==1 || (BROpt==2 && nTransition==1 539 //maybe the TotalGammaRho[i_level] and B 540 if(TotalCumulBR[i_level]==0){ 541 TotalCumulBR[i_level]=new G4double[i_level]; 542 TotalGammaRho[i_level]=ComputeDecayIntensiti 543 } 544 for(G4int j=0;j<i_level;j++){ 545 if(TotalCumulBR[i_level][j]>randnumber){ 546 multipolarity=GetMultipolarity(&theLevels[ 547 return j; 548 } 549 } 550 NuDEXException(__FILE__,std::to_string(_ 551 } 552 //---------------------------------------- 553 554 //BROpt==0 555 //---------------------------------------- 556 // If not, maybe the TotalGammaRho[i_level 557 if(TotalGammaRho[i_level]<0){//not compute 558 TotalGammaRho[i_level]=ComputeDecayInten 559 } 560 theSampledLevel=-1; 561 ComputeDecayIntensities(i_level,0,randnumb 562 multipolarity=theSampledMultipolarity; 563 return theSampledLevel; 564 //---------------------------------------- 565 } 566 567 NuDEXException(__FILE__,std::to_string(__LIN 568 return 0; 569 } 570 571 void G4NuDEXStatisticalNucleus::ChangeLevelSpi 572 573 if(i_level==-1){ //change BR of thermal, ign 574 if(Sn>0 && NLevels>1){ 575 CreateThermalCaptureLevel(seed); 576 GenerateThermalCaptureLevelBR(theLibDir. 577 } 578 return; 579 } 580 581 if(i_level<0 || i_level>=NLevels){ 582 std::cout<<" i_level = "<<i_level<<" ----- 583 NuDEXException(__FILE__,std::to_string(__L 584 } 585 586 //Do not apply to known levels: 587 if(i_level<NKnownLevels || theLevels[i_level 588 std::cout<<" ####### WARNING: you are tryi 589 return; 590 //NuDEXException(__FILE__,std::to_string(_ 591 } 592 593 theLevels[i_level].spinx2=newspinx2; 594 theLevels[i_level].parity=newParity; 595 if(seed>0){ 596 theLevels[i_level].seed=seed; 597 } 598 else{ 599 theLevels[i_level].seed=theRandom2->Intege 600 } 601 if(nlevels>=0){ 602 theLevels[i_level].NLevels=nlevels; 603 } 604 if(width>=0){ 605 theLevels[i_level].Width=width; 606 } 607 608 if(TotalGammaRho[i_level]>=0){ //then we hav 609 G4double* br_vector=0; 610 if(TotalCumulBR!=0){ 611 br_vector=TotalCumulBR[i_level]; 612 } 613 TotalGammaRho[i_level]=ComputeDecayIntensi 614 } 615 616 } 617 618 619 //if randnumber<0, return the total TotalGamma 620 //if randnumber>0, it is assumed that TotalGam 621 // (in the TotalGammaRho[] array or i 622 // The result is stored in theSampledLevel 623 G4double G4NuDEXStatisticalNucleus::ComputeDec 624 625 G4bool ComputeAlsoBR=false; 626 if(cumulativeBR!=0){ComputeAlsoBR=true;} 627 //------------------------------------------ 628 //Some checks: 629 if(i_level>=NLevels || i_level<0){ 630 std::cout<<" ############ Error: i = "<<i_ 631 NuDEXException(__FILE__,std::to_string(__L 632 } 633 if(randnumber>0){ 634 ComputeAlsoBR=false; 635 if(TotGR<=0){ 636 TotGR=TotalGammaRho[i_level]; 637 } 638 if(TotGR<=0){ 639 NuDEXException(__FILE__,std::to_string(_ 640 } 641 } 642 //------------------------------------------ 643 644 theRandom2->SetSeed(theLevels[i_level].seed) 645 G4double thisTotalGammaRho=0; 646 for(G4int j=0;j<i_level;j++){ 647 //If "solape" then zero: 648 if(theLevels[i_level].Energy-theLevels[i_l 649 thisTotalGammaRho+=0; //not cecessary, b 650 if(ComputeAlsoBR){ 651 cumulativeBR[j]=0; 652 } 653 } 654 else{ 655 G4double Eg=theLevels[i_level].Energy-th 656 657 //-------------------------------------- 658 //Check which are allowed transitions: 659 G4bool E1allowed=true,M1allowed=true,E2a 660 G4int Lmin=std::abs(theLevels[i_level].s 661 G4int Lmax=(theLevels[i_level].spinx2+th 662 if(theLevels[i_level].parity==theLevels[ 663 E1allowed=false; 664 } 665 else{ 666 M1allowed=false; E2allowed=false; 667 } 668 if(Lmin>1 || Lmax<1){ 669 E1allowed=false; M1allowed=false; 670 } 671 if(Lmin>2 || Lmax<2){ 672 E2allowed=false; 673 } 674 if(AllowE1){E1allowed=true;} 675 //-------------------------------------- 676 677 G4double GammaRho=0,Sumrand2; 678 theSampledMultipolarity=-50; 679 G4int RealNTransitions=theLevels[i_level 680 681 G4double rand; 682 G4double MaxNSamplesForChi2=1000; 683 684 if(E1allowed){ 685 Sumrand2=RealNTransitions; 686 if(SampleGammaWidths==1){ //Porter-Thomas fl 687 Sumrand2=0; 688 if(RealNTransitions>MaxNSamplesForChi2){ 689 Sumrand2=RealNTransitions*theRandom2->Ga 690 } 691 else{ 692 for(G4int ntr=0;ntr<RealNTransitions;ntr 693 rand=theRandom2->Gaus(0,1); 694 Sumrand2+=rand*rand; 695 } 696 } 697 } 698 GammaRho+=Sumrand2*Eg*Eg*Eg*thePSF->GetE1(Eg 699 if(randnumber>=0){ 700 if(thisTotalGammaRho+GammaRho>=TotGR*randn 701 theSampledMultipolarity=1; 702 } 703 } 704 } 705 if(M1allowed){ 706 Sumrand2=RealNTransitions; 707 if(SampleGammaWidths==1){ //Porter-Thomas fl 708 Sumrand2=0; 709 if(RealNTransitions>MaxNSamplesForChi2){ 710 Sumrand2=RealNTransitions*theRandom2->Ga 711 } 712 else{ 713 for(G4int ntr=0;ntr<RealNTransitions;ntr 714 rand=theRandom2->Gaus(0,1); 715 Sumrand2+=rand*rand; 716 } 717 } 718 } 719 GammaRho+=Sumrand2*Eg*Eg*Eg*thePSF->GetM1(Eg 720 if(randnumber>=0){ 721 if(thisTotalGammaRho+GammaRho>=TotGR*randn 722 theSampledMultipolarity=-1; 723 } 724 } 725 } 726 if(E2allowed){ 727 Sumrand2=RealNTransitions; 728 if(SampleGammaWidths==1){ //Porter-Thomas fl 729 Sumrand2=0; 730 if(RealNTransitions>MaxNSamplesForChi2){ 731 Sumrand2=RealNTransitions*theRandom2->Ga 732 } 733 else{ 734 for(G4int ntr=0;ntr<RealNTransitions;ntr 735 rand=theRandom2->Gaus(0,1); 736 Sumrand2+=rand*rand; 737 } 738 } 739 } 740 GammaRho+=Sumrand2*Eg*Eg*Eg*Eg*Eg*thePSF->Ge 741 if(randnumber>=0){ 742 if(thisTotalGammaRho+GammaRho>=TotGR*randn 743 theSampledMultipolarity=2; 744 } 745 } 746 } 747 748 /* 749 if(i_level==NLevels-1 && j<10){ 750 std::cout<<j<<" "<<GammaRho<<" "<<E1allow 751 } 752 */ 753 754 thisTotalGammaRho+=GammaRho; 755 if(ComputeAlsoBR){ 756 cumulativeBR[j]=GammaRho; 757 } 758 } 759 760 if(randnumber>=0 && thisTotalGammaRho>=Tot 761 if(theSampledMultipolarity==-50){ 762 NuDEXException(__FILE__,std::to_string(__LIN 763 } 764 theSampledLevel=j; 765 return -1; 766 } 767 } 768 769 //If there are no allowed transitions: 770 if(randnumber>=0 && thisTotalGammaRho==0){ / 771 return ComputeDecayIntensities(i_level,0,r 772 } 773 774 if(randnumber>=0){ //then we should not be h 775 NuDEXException(__FILE__,std::to_string(__L 776 } 777 778 if(thisTotalGammaRho>0){ 779 if(ComputeAlsoBR){ 780 for(G4int j=0;j<i_level;j++){ 781 cumulativeBR[j]/=thisTotalGammaRho; 782 } 783 for(G4int j=1;j<i_level;j++){ 784 cumulativeBR[j]+=cumulativeBR[j-1]; 785 } 786 if(std::fabs(cumulativeBR[i_level-1]-1)> 787 std::cout<<" ############### Warning: cumul 788 } 789 } 790 } 791 else{ 792 //std::cout<<" ############### WARNING: to 793 //NuDEXException(__FILE__,std::to_string(_ 794 thisTotalGammaRho=ComputeDecayIntensities( 795 } 796 797 return thisTotalGammaRho; 798 } 799 800 801 //retrieves the "lowest" allowed multipolarity 802 G4int G4NuDEXStatisticalNucleus::GetMultipolar 803 804 if(theInitialLevel->spinx2+theFinalLevel->sp 805 return 0; 806 } 807 G4int Lmin=std::abs(theInitialLevel->spinx2- 808 if(Lmin==0){Lmin=1;} 809 if(Lmin%2==0){ 810 if(theInitialLevel->parity==theFinalLevel- 811 return +Lmin; 812 } 813 else{ 814 return -Lmin; 815 } 816 } 817 else{ 818 if(theInitialLevel->parity==theFinalLevel- 819 return -Lmin; 820 } 821 else{ 822 return +Lmin; 823 } 824 } 825 826 return 0; 827 } 828 829 830 //Retrieves the closest level of the given spi 831 //if spinx2<0, then retrieves the closest leve 832 G4int G4NuDEXStatisticalNucleus::GetClosestLev 833 834 //std::cout<<" XXX finding closest level of 835 836 //------------------------------------------ 837 // We try to go closer to the solution, othe 838 G4int i_down=0,i_up=NLevels-1; 839 G4int i_close_down=0,i_close_up=NLevels-1; 840 G4int i_close=0; 841 while(i_close_up-i_close_down>10){ 842 i_close=(i_close_up+i_close_down)/2; 843 if(theLevels[i_close].Energy>Energy){ 844 i_close_up=i_close; 845 } 846 else{ 847 i_close_down=i_close; 848 } 849 } 850 851 for(G4int i=i_close_up;i<NLevels;i++){ 852 i_up=i; 853 if((theLevels[i].spinx2==spinx2 && theLeve 854 break; 855 } 856 } 857 for(G4int i=i_close_down;i>=0;i--){ 858 i_down=i; 859 if((theLevels[i].spinx2==spinx2 && theLeve 860 break; 861 } 862 } 863 //------------------------------------------ 864 865 G4double MinEnergyDistance=-1,EnergyDistance 866 G4int result=-1; 867 for(G4int i=i_down;i<=i_up;i++){ 868 EnergyDistance=std::fabs(theLevels[i].Ener 869 if((theLevels[i].spinx2==spinx2 && theLeve 870 if(EnergyDistance<MinEnergyDistance || M 871 MinEnergyDistance=EnergyDistance; 872 result=i; 873 } 874 } 875 } 876 //std::cout<<" XXX found --> "<<result<<std: 877 878 879 return result; 880 } 881 882 883 Level* G4NuDEXStatisticalNucleus::GetLevel(G4i 884 885 if(i_level>=0 && i_level<NLevels){ 886 return &theLevels[i_level]; 887 } 888 if(i_level==-1){ 889 return &theThermalCaptureLevel; 890 } 891 892 std::cout<<" ############ WARNING: for ZA="< 893 894 return 0; 895 } 896 897 G4double G4NuDEXStatisticalNucleus::GetLevelEn 898 899 if(i_level>=0 && i_level<NLevels){ 900 return theLevels[i_level].Energy; 901 } 902 903 return -1; 904 } 905 906 907 void G4NuDEXStatisticalNucleus::ComputeKnownLe 908 909 910 for(G4int i=1;i<NKnownLevels;i++){ 911 if(theKnownLevels[i].NGammas!=0){ 912 for(G4int j=0;j<theKnownLevels[i].NGamma 913 theKnownLevels[i].multipolarity[j]=GetMultip 914 } 915 } 916 if(theKnownLevels[i].NGammas==0){ 917 G4double* tmp=new G4double[i]; 918 ComputeDecayIntensities(i,tmp); 919 for(G4int j=1;j<i;j++){ 920 if(tmp[j]!=tmp[j-1]){theKnownLevels[i].NGamm 921 } 922 if(tmp[0]!=0){theKnownLevels[i].NGammas+ 923 if(theKnownLevels[i].NGammas<=0){ 924 NuDEXException(__FILE__,std::to_string(__LIN 925 } 926 927 theKnownLevels[i].FinalLevelID=new G4int 928 theKnownLevels[i].multipolarity=new G4in 929 theKnownLevels[i].Eg=new G4double[theKno 930 theKnownLevels[i].cumulPtot=new G4double 931 theKnownLevels[i].Pg=new G4double[theKno 932 theKnownLevels[i].Pe=new G4double[theKno 933 theKnownLevels[i].Icc=new G4double[theKn 934 G4int cont=0; 935 if(tmp[0]!=0){ 936 theKnownLevels[i].FinalLevelID[cont]=0; 937 theKnownLevels[i].Eg[cont]=theKnownLevels[i] 938 theKnownLevels[i].cumulPtot[cont]=tmp[0]; 939 theKnownLevels[i].multipolarity[cont]=GetMul 940 cont++; 941 } 942 for(G4int j=1;j<i;j++){ 943 if(tmp[j]!=tmp[j-1]){ 944 theKnownLevels[i].FinalLevelID[cont]=j; 945 theKnownLevels[i].Eg[cont]=theKnownLevels[ 946 theKnownLevels[i].cumulPtot[cont]=tmp[j]; 947 theKnownLevels[i].multipolarity[cont]=GetM 948 cont++; 949 } 950 } 951 delete [] tmp; 952 for(G4int j=0;j<theKnownLevels[i].NGamma 953 theKnownLevels[i].Icc[j]=0; 954 if(ElectronConversionFlag==2){ 955 if(theICC){ 956 theKnownLevels[i].Icc[j]=theICC->GetICC( 957 } 958 } 959 } 960 G4double alpha=theKnownLevels[i].Icc[0]; 961 theKnownLevels[i].Pg[0]=theKnownLevels[i 962 theKnownLevels[i].Pe[0]=theKnownLevels[i 963 for(G4int j=1;j<theKnownLevels[i].NGamma 964 alpha=theKnownLevels[i].Icc[j]; 965 theKnownLevels[i].Pg[j]=(theKnownLevels[i].c 966 theKnownLevels[i].Pe[j]=(theKnownLevels[i].c 967 } 968 } 969 } 970 971 972 } 973 974 975 976 977 void G4NuDEXStatisticalNucleus::CreateLevelSch 978 979 //The known levels have been read already 980 NLevels=-1; 981 Level* theUnkonwnLevels=0; 982 if(E_unk_min>=E_unk_max){//Then we know all 983 NUnknownLevels=0; //will be updated to 1 w 984 } 985 else{ 986 //======================================== 987 //Unknown levels: 988 G4int maxarraysize=EstimateNumberOfLevelsT 989 do{ 990 maxarraysize*=2; 991 if(theUnkonwnLevels!=0){delete [] theUnk 992 //std::cout<<" Max array size of "<<maxa 993 theUnkonwnLevels=new Level[maxarraysize] 994 NUnknownLevels=GenerateAllUnknownLevels( 995 }while(NUnknownLevels<0); 996 //======================================== 997 } 998 999 //========================================== 1000 //We define the final level scheme: 1001 NLevels=NKnownLevels+NUnknownLevels; 1002 //std::cout<<" There are "<<NLevels<<" leve 1003 theLevels=new Level[NLevels]; 1004 for(G4int i=0;i<NKnownLevels;i++){ 1005 CopyLevel(&theKnownLevels[i],&theLevels[i 1006 } 1007 for(G4int i=0;i<NUnknownLevels;i++){ 1008 CopyLevel(&theUnkonwnLevels[i],&theLevels 1009 } 1010 delete [] theUnkonwnLevels; 1011 //========================================= 1012 1013 //Final check: 1014 G4int TotalNIndividualLevels=1; 1015 for(G4int i=1;i<NLevels;i++){ 1016 TotalNIndividualLevels+=theLevels[i].NLev 1017 if(theLevels[i-1].Energy>theLevels[i].Ene 1018 std::cout<<" ########### Error creating 1019 NuDEXException(__FILE__,std::to_string( 1020 } 1021 } 1022 1023 std::cout<<" NuDEX: Generated statistical n 1024 1025 } 1026 1027 1028 void G4NuDEXStatisticalNucleus::CreateThermal 1029 1030 1031 G4int capturespinx2=((std::fabs(I0)+0.5)*2+ 1032 G4bool capturepar=true; if(I0<0){capturepar 1033 theThermalCaptureLevel.Energy=Sn; 1034 theThermalCaptureLevel.spinx2=capturespinx2 1035 theThermalCaptureLevel.parity=capturepar; 1036 if(seed>0){ 1037 theThermalCaptureLevel.seed=seed; 1038 } 1039 else{ 1040 theThermalCaptureLevel.seed=theRandom2->I 1041 } 1042 theThermalCaptureLevel.KnownLevelID=-1; 1043 theThermalCaptureLevel.NLevels=1; 1044 theThermalCaptureLevel.Width=0; 1045 1046 NLevelsBelowThermalCaptureLevel=0; 1047 for(G4int i=0;i<NLevels;i++){ 1048 if(theLevels[i].Energy<theThermalCaptureL 1049 NLevelsBelowThermalCaptureLevel++; 1050 } 1051 } 1052 NLevelsBelowThermalCaptureLevel--; // we ex 1053 1054 if(NLevelsBelowThermalCaptureLevel<=0){ 1055 NLevelsBelowThermalCaptureLevel=1; //we c 1056 } 1057 1058 } 1059 1060 G4int G4NuDEXStatisticalNucleus::GenerateLeve 1061 1062 //If A_Int even/odd --> spinx2 (spin_val*2) 1063 if(((A_Int+spinx2)%2)!=0){ 1064 return 0; 1065 } 1066 1067 //Get the level density: 1068 G4double meanNLevels=theLD->Integrate(Emin, 1069 1070 //Sample total number of levels: ???????? 1071 G4int thisNLevels=0; 1072 if(meanNLevels>0){ 1073 thisNLevels=(G4int)theRandom1->Poisson(me 1074 } 1075 1076 if(thisNLevels>=MaxNLevelsToFill){ 1077 std::cout<<" Warning: not enough space to 1078 return -1; 1079 } 1080 1081 //Distribute the levels in the energy inter 1082 for(G4int i=0;i<thisNLevels;i++){ 1083 someLevels[i].Energy=theRandom1->Uniform( 1084 someLevels[i].spinx2=spinx2; 1085 someLevels[i].parity=parity; 1086 someLevels[i].seed=0; 1087 someLevels[i].KnownLevelID=-1; 1088 someLevels[i].NLevels=1; 1089 someLevels[i].Width=0; 1090 } 1091 1092 return thisNLevels; 1093 } 1094 1095 G4int G4NuDEXStatisticalNucleus::GenerateLeve 1096 1097 G4int TotalNLevels=0; 1098 G4int NIntervals=1000; 1099 1100 // When the LD changes significantly betwee 1101 // So we will sample the levels according t 1102 // rho(E+<D>)/rho(E) is small, at higher en 1103 // The difference between "big" and "small" 1104 G4double WignerRatioThreshold=2; 1105 G4double LevDenThreshold=1; //if there is l 1106 1107 for(G4int i=0;i<NIntervals;i++){ 1108 G4double emin=Emin+(Emax-Emin)*i/(G4doubl 1109 G4double emax=Emin+(Emax-Emin)*(i+1)/(G4d 1110 G4double meanene=(emin+emax)/2.; 1111 G4double LevDen1=theLD->GetLevelDensity(m 1112 if(LevDen1>LevDenThreshold){ 1113 G4double LevDen2=theLD->GetLevelDensity 1114 if(LevDen2/LevDen1<WignerRatioThreshold 1115 //std::cout<<" Wigner way to generate level 1116 G4int newExtraLevels=GenerateWignerLevels(e 1117 if(newExtraLevels<0){return -1;} 1118 TotalNLevels+=newExtraLevels; 1119 break; 1120 } 1121 } 1122 //then use Poisson: 1123 G4int newExtraLevels=GenerateLevelsInSmal 1124 if(newExtraLevels<0){return -1;} 1125 TotalNLevels+=newExtraLevels; 1126 } 1127 1128 return TotalNLevels; 1129 } 1130 1131 1132 //Wigner law: p(x)=pi/2*rho*x*exp(-pi/4*rho*r 1133 G4int G4NuDEXStatisticalNucleus::GenerateWign 1134 1135 //If A_Int even/odd --> spinx2 (spin_val*2) 1136 if(((A_Int+spinx2)%2)!=0){ 1137 return 0; 1138 } 1139 1140 G4int TotalNLevels=0; 1141 1142 G4double previousELevel=Emin,nextELevel; 1143 while(previousELevel<Emax){ 1144 G4double LevDen=theLD->GetLevelDensity(pr 1145 G4double arandgamma=theRandom1->Uniform() 1146 G4double DeltaEMultipliedByLevDen=std::sq 1147 nextELevel=previousELevel+DeltaEMultiplie 1148 if(nextELevel<Emax){ 1149 someLevels[TotalNLevels].Energy=nextELe 1150 someLevels[TotalNLevels].spinx2=spinx2; 1151 someLevels[TotalNLevels].parity=parity; 1152 someLevels[TotalNLevels].seed=0; 1153 someLevels[TotalNLevels].KnownLevelID=- 1154 someLevels[TotalNLevels].NLevels=1; 1155 someLevels[TotalNLevels].Width=0; 1156 TotalNLevels++; 1157 if(TotalNLevels>=MaxNLevelsToFill){ 1158 std::cout<<" Warning: not enough space to f 1159 //NuDEXException(__FILE__,std::to_string(__ 1160 return -1; 1161 } 1162 } 1163 previousELevel=nextELevel; 1164 } 1165 1166 return TotalNLevels; 1167 1168 } 1169 1170 1171 //We genereate the levels directly in bands, 1172 G4int G4NuDEXStatisticalNucleus::GenerateBand 1173 1174 //If A_Int even/odd --> spinx2 (spin_val*2) 1175 if(((A_Int+spinx2)%2)!=0){ 1176 return 0; 1177 } 1178 1179 G4double Emin=Emin_bands; 1180 G4double Emax=Emax_bands; 1181 G4int TotalNLevels=0; 1182 1183 if(bandmax>NBands-1){ 1184 NuDEXException(__FILE__,std::to_string(__ 1185 } 1186 1187 for(G4int i=bandmin;i<=bandmax;i++){ 1188 G4double emin=Emin+(Emax-Emin)*i/(G4doubl 1189 G4double emax=Emin+(Emax-Emin)*(i+1.)/(G4 1190 G4double AverageNumberOfLevels=theLD->Int 1191 G4int NumberOfLevelsInThisBand=0; 1192 if(AverageNumberOfLevels>0){ 1193 NumberOfLevelsInThisBand=(G4int)theRand 1194 } 1195 if(NumberOfLevelsInThisBand>0){ 1196 someLevels[TotalNLevels].Energy=(emax+e 1197 someLevels[TotalNLevels].spinx2=spinx2; 1198 someLevels[TotalNLevels].parity=parity; 1199 someLevels[TotalNLevels].seed=0; 1200 someLevels[TotalNLevels].KnownLevelID=- 1201 someLevels[TotalNLevels].NLevels=Number 1202 someLevels[TotalNLevels].Width=emax-emi 1203 TotalNLevels++; 1204 if(TotalNLevels>=MaxNLevelsToFill){ 1205 std::cout<<" Warning: not enough space to f 1206 //NuDEXException(__FILE__,std::to_string(__ 1207 return -1; 1208 } 1209 } 1210 } 1211 1212 return TotalNLevels; 1213 } 1214 1215 1216 1217 G4int G4NuDEXStatisticalNucleus::GenerateAllU 1218 1219 G4int TotalNLevels=0,NLev; 1220 if(E_unk_min>=E_unk_max){return 0;} 1221 1222 for(G4int spinx2=0;spinx2<=maxspinx2;spinx2 1223 for(G4int ipar=0;ipar<2;ipar++){ 1224 //If A_Int even/odd --> spinx2 (spin_va 1225 if(((A_Int+spinx2)%2)==0){ 1226 //----------------------------------------- 1227 G4bool parity=true; 1228 if(ipar==1){parity=false;} 1229 1230 //We create random levels between E_unk_min 1231 //We will create the levels one by one at l 1232 //The limit between the two ranges will be 1233 G4double Emin=E_unk_min; 1234 G4double Emax=E_unk_max; 1235 G4double E_lim_onebyone=2.*E_unk_max; 1236 G4int i_Band_E_lim_onebyone=NBands+1; // ba 1237 1238 //----------------------------------------- 1239 //Calculate E_lim_onebyone: 1240 #ifndef GENERATEEXPLICITLYALLLEVELSCHEME 1241 if(NBands>0){ 1242 if(MinLevelsPerBand<=0){ // All the level 1243 E_lim_onebyone=0; 1244 i_Band_E_lim_onebyone=0; 1245 } 1246 else{ 1247 G4double bandwidth=(Emax_bands-Emin_ban 1248 G4double rho_lim_onebyone=3.*(MinLevels 1249 E_lim_onebyone=theLD->EstimateInverse(r 1250 } 1251 } 1252 if(E_unk_max-Emax_bands>0.001){ //then E_un 1253 E_lim_onebyone=2.*E_unk_max; 1254 i_Band_E_lim_onebyone=NBands+1; 1255 } 1256 1257 // E_lim_onebyone should be in a limit betw 1258 if(E_lim_onebyone>E_unk_min && E_lim_onebyo 1259 for(G4int i=0;i<NBands;i++){ 1260 G4double elow_band=Emin_bands+(Emax_ban 1261 if(elow_band>E_lim_onebyone){ 1262 E_lim_onebyone=elow_band; 1263 i_Band_E_lim_onebyone=i; 1264 break; 1265 } 1266 } 1267 } 1268 #endif 1269 //----------------------------------------- 1270 1271 1272 if(E_lim_onebyone>E_unk_min){ //then we hav 1273 if(E_lim_onebyone<Emax){ 1274 Emax=E_lim_onebyone; 1275 } 1276 NLev=GenerateLevelsInBigRange(Emin,Emax,s 1277 if(NLev<0){return -1;} 1278 if(NBands>0 && NLev>0){ 1279 NLev=CreateBandsFromLevels(NLev,&(someL 1280 } 1281 TotalNLevels+=NLev; 1282 } 1283 1284 if(i_Band_E_lim_onebyone<NBands){ //then we 1285 NLev=GenerateBandLevels(i_Band_E_lim_oneb 1286 if(NLev<0){return -1;} 1287 TotalNLevels+=NLev; 1288 } 1289 //----------------------------------------- 1290 } 1291 } 1292 } 1293 1294 1295 //Order levels by energy: 1296 qsort(someLevels,TotalNLevels,sizeof(Level) 1297 1298 return TotalNLevels; 1299 } 1300 1301 //Junta varios niveles en uno solo, creando d 1302 //Entiende que no hay otros spines ni paridad 1303 //Devuelve el numero de niveles actualizado. 1304 G4int G4NuDEXStatisticalNucleus::CreateBandsF 1305 1306 G4double Emin=Emin_bands; 1307 G4double Emax=Emax_bands; 1308 1309 Level* theBandLevels=new Level[NBands]; 1310 for(G4int i=0;i<NBands;i++){ 1311 G4double emin=Emin+(Emax-Emin)*i/(G4doubl 1312 G4double emax=Emin+(Emax-Emin)*(i+1.)/(G4 1313 theBandLevels[i].Energy=(emax+emin)/2.; 1314 theBandLevels[i].spinx2=spinx2; 1315 theBandLevels[i].parity=parity; 1316 theBandLevels[i].seed=0; 1317 theBandLevels[i].KnownLevelID=-1; 1318 theBandLevels[i].NLevels=0; 1319 theBandLevels[i].Width=emax-emin; 1320 G4int NLevelsInThisBand=0; 1321 for(G4int j=0;j<thisNLevels;j++){ 1322 if(someLevels[j].spinx2!=spinx2 || some 1323 NuDEXException(__FILE__,std::to_string(__LI 1324 } 1325 if(someLevels[j].Energy>=emin && someLe 1326 NLevelsInThisBand+=someLevels[j].NLevels; 1327 } 1328 } 1329 if(NLevelsInThisBand>=MinLevelsPerBand){ 1330 for(G4int j=0;j<thisNLevels;j++){ 1331 if(someLevels[j].Energy>=emin && someLevels 1332 theBandLevels[i].NLevels+=someLevels[j].N 1333 someLevels[j].Energy=-1; 1334 } 1335 } 1336 } 1337 } 1338 1339 G4int FinalNBands=NBands; 1340 1341 //Eliminate bands with cero levels: 1342 for(G4int i=0;i<FinalNBands;i++){ 1343 if(theBandLevels[i].NLevels==0){ 1344 if(i!=FinalNBands-1){ 1345 CopyLevel(&theBandLevels[FinalNBands-1],&th 1346 } 1347 i--; 1348 FinalNBands--; 1349 } 1350 } 1351 1352 //Replace levels with bands and update numb 1353 G4int NbandsCopied=0; 1354 for(G4int i=0;i<thisNLevels;i++){ 1355 if(someLevels[i].Energy<0){ 1356 if(NbandsCopied<FinalNBands){ //this le 1357 CopyLevel(&theBandLevels[NbandsCopied],&som 1358 NbandsCopied++; 1359 } 1360 else{ //there is no band to replace. Co 1361 CopyLevel(&someLevels[thisNLevels-1],&someL 1362 i--; 1363 thisNLevels--; 1364 } 1365 } 1366 } 1367 1368 //Simple check: 1369 if(NbandsCopied!=FinalNBands){ 1370 NuDEXException(__FILE__,std::to_string(__ 1371 } 1372 delete [] theBandLevels; 1373 return thisNLevels; 1374 } 1375 1376 1377 //Esto lo que hace es intentar reemplazar niv 1378 G4int G4NuDEXStatisticalNucleus::InsertHighEn 1379 1380 1381 1382 G4bool* HasBeenInserted=new G4bool[KnownLev 1383 for(G4int i=0;i<KnownLevelsVectorSize;i++){ 1384 HasBeenInserted[i]=false; 1385 } 1386 1387 for(G4int kk=0;kk<2;kk++){ //loop two times 1388 for(G4int k=1;k<5;k++){ //loop in the dis 1389 G4double MaxEnergyDistance=0.1*k; //The 1390 for(G4int i=NKnownLevels;i<KnownLevelsV 1391 if(theKnownLevels[i].Energy>Sn){break;} 1392 if(HasBeenInserted[i]==false && (theKnownLe 1393 //--------------------------------------- 1394 G4double MinEnergyDistance=-1,EnergyDista 1395 G4int thespinx2=theKnownLevels[i].spinx2; 1396 G4bool thepar=theKnownLevels[i].parity; 1397 G4int unknownLevelID=-1; 1398 for(G4int j=NKnownLevels;j<NLevels-1;j++) 1399 if(theLevels[j].spinx2==thespinx2 && th 1400 EnergyDistance=std::fabs(theKnownLeve 1401 if((EnergyDistance<MinEnergyDistance 1402 MinEnergyDistance=EnergyDistance; 1403 unknownLevelID=j; 1404 } 1405 //else if(EnergyDistance>MinEnergyDis 1406 //break; 1407 //} 1408 } 1409 } 1410 if(unknownLevelID>0 && theLevels[unknownL 1411 //std::cout<<" Copy Level "<<i<<" with 1412 CopyLevel(&theKnownLevels[i],&theLevels 1413 theLevels[unknownLevelID].KnownLevelID= 1414 HasBeenInserted[i]=true; 1415 } 1416 //--------------------------------------- 1417 } 1418 } 1419 } 1420 } 1421 delete [] HasBeenInserted; 1422 1423 //We re-order the levels: 1424 qsort(theLevels,NLevels,sizeof(Level), Comp 1425 1426 1427 1428 //----------------------------------------- 1429 //we cannot go to from an inserted level wi 1430 //so, if it is the case, we change the fina 1431 for(G4int i=NKnownLevels;i<NLevels;i++){ 1432 if(theLevels[i].KnownLevelID>0){ 1433 G4int knownID=theLevels[i].KnownLevelID 1434 for(G4int j=0;j<theKnownLevels[knownID] 1435 if(theKnownLevels[knownID].FinalLevelID[j]> 1436 //--------------------------------------- 1437 G4int i_finalknownlevel=theKnownLevels[kn 1438 G4double MinEnergyDistance=-1; 1439 G4int i_statlevel=-1; 1440 for(G4int k=0;k<i;k++){ 1441 G4double EnergyDistance=std::fabs(theKn 1442 if(EnergyDistance<MinEnergyDistance || 1443 MinEnergyDistance=EnergyDistance; 1444 i_statlevel=k; 1445 } 1446 } 1447 if(MinEnergyDistance<0){ 1448 NuDEXException(__FILE__,std::to_string( 1449 } 1450 //std::cout<<" Final known level "<<i_fin 1451 if(theLevels[i_statlevel].KnownLevelID>=0 1452 theKnownLevels[knownID].FinalLevelID[j] 1453 } 1454 else{ //is a real statistical level. We c 1455 theKnownLevels[knownID].FinalLevelID[j] 1456 theKnownLevels[knownID].multipolarity[j 1457 theKnownLevels[knownID].Eg[j]=theLevels 1458 theKnownLevels[knownID].Pg[j]=theKnownL 1459 theKnownLevels[knownID].Pe[j]=0; 1460 theKnownLevels[knownID].Icc[j]=0; //set 1461 } 1462 //--------------------------------------- 1463 } 1464 } 1465 } 1466 } 1467 //----------------------------------------- 1468 1469 return 0; 1470 } 1471 1472 1473 1474 G4int G4NuDEXStatisticalNucleus::EstimateNumb 1475 1476 if(E_unk_min>=E_unk_max){return 0;} 1477 1478 #ifndef GENERATEEXPLICITLYALLLEVELSCHEME 1479 if(BandWidth>0){ 1480 return maxspinx2*NBands*MinLevelsPerBand; 1481 } 1482 #endif 1483 1484 G4double Emin=E_unk_min; 1485 G4double Emax=E_unk_max; 1486 G4double TotalNLevels=0; 1487 G4double TotalNLevelsInsideBands=0; 1488 G4double MaxNLevelsWithSameSpinAndParity=0; 1489 G4double emin,emax,meanEnergy,LevDen,meanNL 1490 G4int NIntervals=1000; 1491 for(G4int spinx2=0;spinx2<=maxspinx2;spinx2 1492 G4double TotalNLevesWithSameSpin=0; 1493 for(G4int i=0;i<NIntervals;i++){ 1494 emin=Emin+(Emax-Emin)*i/(G4double)NInte 1495 emax=Emin+(Emax-Emin)*(i+1)/(G4double)N 1496 meanEnergy=(emax+emin)/2.; 1497 1498 //Positive parity: 1499 LevDen=theLD->GetLevelDensity(meanEnerg 1500 meanNLevels=LevDen*(emax-emin); 1501 TotalNLevels+=meanNLevels; 1502 TotalNLevesWithSameSpin+=meanNLevels; 1503 if(NBands>0 && meanEnergy>Emin_bands && 1504 TotalNLevelsInsideBands+=meanNLevels; 1505 } 1506 1507 1508 //Negative parity: 1509 LevDen=theLD->GetLevelDensity(meanEnerg 1510 meanNLevels=LevDen*(emax-emin); 1511 TotalNLevels+=meanNLevels; 1512 TotalNLevesWithSameSpin+=meanNLevels; 1513 if(NBands>0 && meanEnergy>Emin_bands && 1514 TotalNLevelsInsideBands+=meanNLevels; 1515 } 1516 1517 } 1518 if(MaxNLevelsWithSameSpinAndParity<TotalN 1519 } 1520 1521 MaxNLevelsWithSameSpinAndParity/=2.; 1522 1523 if(TotalNLevelsInsideBands>0){ 1524 //then there are some levels inside bands 1525 G4double TotalNLevelsOutsideBands=TotalNL 1526 G4double MaxNLevelsNeededForBands=NBands* 1527 TotalNLevels=MaxNLevelsWithSameSpinAndPar 1528 } 1529 1530 return (G4int)TotalNLevels; 1531 } 1532 1533 1534 G4double G4NuDEXStatisticalNucleus::TakeTarge 1535 1536 std::ifstream in(fname); 1537 if(!in.good()){ 1538 std::cout<<" ######## Error opening file 1539 NuDEXException(__FILE__,std::to_string(__ 1540 } 1541 check=0; 1542 1543 char buffer[1000]; 1544 G4int aZ,aA; 1545 while(in.get(buffer,6)){ 1546 in.get(buffer,6); aA=atoi(buffer); 1547 in.get(buffer,6); aZ=atoi(buffer); 1548 if(aZ==Z_Int && aA==A_Int-1){ 1549 break; 1550 } 1551 in.ignore(10000,'\n'); 1552 } 1553 if(!in.good()){ 1554 in.close(); 1555 check=-1; 1556 //NuDEXException(__FILE__,std::to_string( 1557 } 1558 in.ignore(10000,'\n'); 1559 G4double spin,par; 1560 in.get(buffer,16); 1561 in.get(buffer,6); spin=std::fabs(atof(buf 1562 in.get(buffer,4); par=atof(buffer); 1563 1564 in.close(); 1565 1566 if(par<0){return -spin;} 1567 1568 return spin; 1569 } 1570 1571 G4double G4NuDEXStatisticalNucleus::ReadKnown 1572 1573 1574 std::ifstream in(fname); 1575 if(!in.good()){ 1576 std::cout<<" ######## Error opening file 1577 NuDEXException(__FILE__,std::to_string(__ 1578 } 1579 1580 char buffer[1000]; 1581 G4int aZ,aA; 1582 while(in.get(buffer,6)){ 1583 in.get(buffer,6); aA=atoi(buffer); 1584 in.get(buffer,6); aZ=atoi(buffer); 1585 if(aZ==Z_Int && aA==A_Int){ 1586 in.get(buffer,6); KnownLevelsVectorSize 1587 in.get(buffer,16); 1588 in.get(buffer,13); G4double newSn=atof( 1589 if(Sn>0 && std::fabs(Sn-newSn)>0.01){ 1590 std::cout<<" ######## WARNING: Sn value fro 1591 } 1592 else if(Sn<0){ 1593 Sn=newSn; 1594 } 1595 break; 1596 } 1597 in.ignore(10000,'\n'); 1598 } 1599 1600 if(!in.good()){ 1601 in.close(); return -1; 1602 } 1603 1604 in.ignore(10000,'\n'); 1605 1606 NKnownLevels=0; 1607 theKnownLevels=new KnownLevel[KnownLevelsVe 1608 for(G4int i=0;i<KnownLevelsVectorSize;i++){ 1609 G4double spin,par; 1610 for(G4int i=0;i<KnownLevelsVectorSize;i++){ 1611 in.get(buffer,4); theKnownLevels[i].id= 1612 in.get(buffer,2); 1613 in.get(buffer,11); theKnownLevels[i].Ene 1614 in.get(buffer,2); 1615 in.get(buffer,6); spin=atof(buffer); 1616 in.get(buffer,4); par=atof(buffer); 1617 if((spin<0 || par==0) && theKnownLevels[i 1618 std::cout<<" ######## WARNING: Spin and 1619 if(spin<0){ 1620 spin=0; 1621 if(i>1){ //Random spin, same as one of the 1622 G4int i_sampleSpin=theRandom1->Integer(i- 1623 spin=theKnownLevels[i_sampleSpin].spinx2/ 1624 } 1625 } 1626 if(par==0){ 1627 par=1; 1628 if(theRandom1->Uniform(-1,1)<0){par=-1;} 1629 } 1630 } 1631 in.get(buffer,2); 1632 in.get(buffer,11); theKnownLevels[i].T12 1633 in.get(buffer,4); theKnownLevels[i].NGa 1634 if(theKnownLevels[i].NGammas>0){ 1635 if(spin<0){ 1636 spin=0; 1637 if(i>1){ //Random spin, same as one of the 1638 G4int i_sampleSpin=theRandom1->Integer(i- 1639 spin=theKnownLevels[i_sampleSpin].spinx2/ 1640 } 1641 } 1642 if(par==0){ 1643 par=1; 1644 if(theRandom1->Uniform(-1,1)<0){par=-1;} 1645 } 1646 } 1647 theKnownLevels[i].spinx2=(G4int)(spin*2+0 1648 if(par>0){theKnownLevels[i].parity=true;} 1649 1650 //--------------------------------- 1651 //decay modes: 1652 in.get(buffer,27); //dummy 1653 in.get(buffer,4); 1654 G4int decays=theKnownLevels[i].Ndecays=at 1655 if(decays>0){ 1656 theKnownLevels[i].decayFraction=new G4d 1657 theKnownLevels[i].decayMode=new std::st 1658 } 1659 for(G4int j=0;j<decays;j++){ 1660 in.get(buffer,5); 1661 in.get(buffer,11); 1662 theKnownLevels[i].decayFraction[j]=atof 1663 in.get(buffer,2); 1664 in.get(buffer,8); 1665 theKnownLevels[i].decayMode[j]=std::str 1666 } 1667 //---------------------------------- 1668 1669 in.ignore(10000,'\n'); 1670 1671 if(theKnownLevels[i].NGammas>0){ 1672 theKnownLevels[i].FinalLevelID=new G4in 1673 theKnownLevels[i].multipolarity=new G4i 1674 theKnownLevels[i].Eg=new G4double[theKn 1675 theKnownLevels[i].cumulPtot=new G4doubl 1676 theKnownLevels[i].Pg=new G4double[theKn 1677 theKnownLevels[i].Pe=new G4double[theKn 1678 theKnownLevels[i].Icc=new G4double[theK 1679 } 1680 1681 for(G4int j=0;j<theKnownLevels[i].NGammas 1682 in.get(buffer,40); 1683 1684 in.get(buffer,5); theKnownLevels[i].Fi 1685 theKnownLevels[i].multipolarity[j]=0; 1686 in.get(buffer,2); 1687 in.get(buffer,11); theKnownLevels[i].E 1688 in.get(buffer,2); 1689 in.get(buffer,11); theKnownLevels[i].P 1690 in.get(buffer,2); 1691 in.get(buffer,11); theKnownLevels[i].P 1692 in.get(buffer,2); 1693 in.get(buffer,11); theKnownLevels[i].I 1694 theKnownLevels[i].cumulPtot[j]=theKnown 1695 in.ignore(10000,'\n'); 1696 if(theKnownLevels[i].FinalLevelID[j]>=i 1697 std::cout<<" ######## Error reading file "< 1698 NuDEXException(__FILE__,std::to_string(__LI 1699 } 1700 } 1701 1702 if(theKnownLevels[i].id!=i || !in.good()) 1703 std::cout<<" ######## Error reading fil 1704 std::cout<<" Level "<<i<<" has id = "<< 1705 NuDEXException(__FILE__,std::to_string( 1706 } 1707 1708 //--------------------------------------- 1709 //normalize, and put cumulPtot as cumulat 1710 G4double totalEmissionProb=0; 1711 for(G4int j=0;j<theKnownLevels[i].NGammas 1712 totalEmissionProb+=theKnownLevels[i].cu 1713 } 1714 //------------------------------------ 1715 if(totalEmissionProb==0){//sometimes all 1716 for(G4int j=0;j<theKnownLevels[i].NGamm 1717 theKnownLevels[i].cumulPtot[j]=1./theKnownL 1718 theKnownLevels[i].Pg[j]=theKnownLevels[i].c 1719 } 1720 totalEmissionProb=1; 1721 } 1722 //------------------------------------ 1723 for(G4int j=0;j<theKnownLevels[i].NGammas 1724 theKnownLevels[i].cumulPtot[j]/=totalEm 1725 theKnownLevels[i].Pg[j]/=totalEmissionP 1726 theKnownLevels[i].Pe[j]=theKnownLevels[ 1727 } 1728 for(G4int j=1;j<theKnownLevels[i].NGammas 1729 theKnownLevels[i].cumulPtot[j]+=theKnow 1730 } 1731 //--------------------------------------- 1732 1733 if(theKnownLevels[i].Energy<=Ecrit*1.0001 1734 NKnownLevels++; 1735 } 1736 /* 1737 else{ 1738 break; 1739 } 1740 */ 1741 } 1742 1743 if(!in.good()){ 1744 in.close(); return -1; 1745 } 1746 1747 in.close(); 1748 1749 return 0; 1750 } 1751 1752 1753 1754 G4double G4NuDEXStatisticalNucleus::ReadEcrit 1755 1756 std::ifstream in(fname); 1757 if(!in.good()){ 1758 std::cout<<" ######## Error opening file 1759 NuDEXException(__FILE__,std::to_string(__ 1760 } 1761 1762 G4int aZ,aA; 1763 char word[200]; 1764 Ecrit=-1; 1765 for(G4int i=0;i<4;i++){in.ignore(10000,'\n' 1766 while(in>>aZ>>aA){ 1767 if(aZ==Z_Int && aA==A_Int){ 1768 in>>word>>word>>word>>word>>word>>word> 1769 } 1770 in.ignore(10000,'\n'); 1771 } 1772 in.close(); 1773 1774 return Ecrit; 1775 } 1776 1777 G4int G4NuDEXStatisticalNucleus::ReadSpecialI 1778 1779 std::string word; 1780 std::ifstream in(fname); 1781 if(!in.good()){ 1782 in.close(); 1783 return -1; 1784 } 1785 G4double MaxSpin; 1786 while(in>>word){ 1787 if(word.c_str()[0]=='#'){in.ignore(10000, 1788 if(word==std::string("END")){break;} 1789 //now we will take values only if they ha 1790 else if(word==std::string("LEVELDENSITYTY 1791 else if(word==std::string("MAXSPIN")){if( 1792 else if(word==std::string("MINLEVELSPERBA 1793 else if(word==std::string("BANDWIDTH_MEV" 1794 else if(word==std::string("MAXEXCENERGY_M 1795 else if(word==std::string("ECRIT_MEV")){i 1796 else if(word==std::string("KNOWNLEVELSFLA 1797 1798 else if(word==std::string("PSF_FLAG")){if 1799 else if(word==std::string("BROPTION")){if 1800 else if(word==std::string("SAMPLEGAMMAWID 1801 1802 else if(word==std::string("ELECTRONCONVER 1803 else if(word==std::string("PRIMARYTHCAPGA 1804 else if(word==std::string("PRIMARYGAMMASE 1805 } 1806 in.close(); 1807 return 1; 1808 } 1809 1810 G4int G4NuDEXStatisticalNucleus::ReadGeneralS 1811 1812 std::ifstream in(fname); 1813 if(!in.good()){ 1814 std::cout<<" ######## Error opening file 1815 NuDEXException(__FILE__,std::to_string(__ 1816 } 1817 char line[1000]; 1818 in.getline(line,1000); 1819 in.getline(line,1000); 1820 G4int tmpZ,tmpA,tmpLDtype,tmpPSFflag,tmpMax 1821 G4double tmpBandWidth,tmpMaxExcEnergy; 1822 unsigned int tmpseed1,tmpseed2,tmpseed3; 1823 G4int finalLDtype=0,finalPSFflag=0,finalMax 1824 G4double finalBandWidth=0,finalMaxExcEnergy 1825 unsigned int finalseed1=0,finalseed2=0,fina 1826 G4bool NuclDataHasBeenRead=false; 1827 G4bool DefaulDataHasBeenRead=false; 1828 while(in>>tmpZ>>tmpA>>tmpLDtype>>tmpPSFflag 1829 G4bool TakeData=false; 1830 if(tmpZ==Z_Int && tmpA==A_Int){ //then th 1831 NuclDataHasBeenRead=true; 1832 TakeData=true; 1833 } 1834 else if(tmpZ==0 && tmpA==0 && NuclDataHas 1835 DefaulDataHasBeenRead=true; 1836 TakeData=true; 1837 } 1838 if(TakeData){ 1839 finalLDtype=tmpLDtype; finalPSFflag=tmp 1840 finalBandWidth=tmpBandWidth; finalMaxEx 1841 finalseed1=tmpseed1; finalseed2=tmpseed 1842 } 1843 } 1844 in.close(); 1845 1846 if(NuclDataHasBeenRead==false && DefaulData 1847 std::cout<<" ######## Error reading "<<fn 1848 NuDEXException(__FILE__,std::to_string(__ 1849 } 1850 1851 // Replace data only if it has not been set 1852 if(LevelDensityType<0){LevelDensityType=fin 1853 if(PSFflag<0){PSFflag=finalPSFflag;} 1854 if(maxspinx2<0){maxspinx2=(G4int)(2.*finalM 1855 if(MinLevelsPerBand<0){MinLevelsPerBand=fin 1856 if(BandWidth==0){BandWidth=finalBandWidth;} 1857 if(MaxExcEnergy==0){MaxExcEnergy=finalMaxEx 1858 if(BROpt<0){BROpt=finalBrOption;} 1859 if(SampleGammaWidths<0){SampleGammaWidths=f 1860 if(Rand1seedProvided==false){seed1=finalsee 1861 if(Rand2seedProvided==false){seed2=finalsee 1862 if(Rand3seedProvided==false){seed3=finalsee 1863 1864 //----------------------------------------- 1865 //Now some checks: 1866 if(maxspinx2<1){ 1867 std::cout<<" ######## Error: maximum spin 1868 } 1869 if(MinLevelsPerBand<=0 && BandWidth<0){ 1870 std::cout<<" ######## Error: MinLevelsPer 1871 } 1872 if(BROpt!=0 && BROpt!=1 && BROpt!=2){ 1873 std::cout<<" ######## Error: BROpt for ge 1874 } 1875 if(SampleGammaWidths!=0 && SampleGammaWidth 1876 std::cout<<" ######## Error: SampleGammaW 1877 } 1878 //----------------------------------------- 1879 1880 1881 return 0; 1882 } 1883 1884 1885 void G4NuDEXStatisticalNucleus::GenerateTherm 1886 1887 char fname[1000]; 1888 snprintf(fname,1000,"%s/PrimaryCaptureGamma 1889 1890 G4int aZA,ng=0; 1891 char word[200]; 1892 G4double ELevel; 1893 G4double ThEg[1000],ThI[1000]; 1894 //G4double TotalThI=0; 1895 1896 //We read the gamma intensities from the fi 1897 std::ifstream in(fname); 1898 while(in>>word){ 1899 if(word[0]=='Z' && word[1]=='A'){ 1900 in>>aZA; 1901 if(aZA==Z_Int*1000+A_Int){ 1902 in.ignore(10000,'\n'); 1903 in>>ELevel>>ng; 1904 in.ignore(10000,'\n'); 1905 ELevel*=1.e-3; // keV to MeV 1906 //Check if ELevel is very close to Sn: 1907 if(ELevel/Sn>1.001 || ELevel/Sn<0.999){ 1908 std::cout<<" ########## WARNING: ELevel = 1909 } 1910 for(G4int i=0;i<ng;i++){ 1911 in>>ThEg[i]>>ThI[i]; 1912 ThEg[i]/=1.e3; // keV to MeV 1913 ThI[i]/=100.; // percent to no-percent 1914 ThI[i]*=PrimaryGammasIntensityNormFactor; 1915 //TotalThI+=ThI[i]; 1916 } 1917 break; 1918 } 1919 } 1920 in.ignore(10000,'\n'); 1921 } 1922 in.close(); 1923 1924 if(theThermalCaptureLevelCumulBR){delete [] 1925 theThermalCaptureLevelCumulBR=new G4double[ 1926 for(G4int i=0;i<NLevelsBelowThermalCaptureL 1927 theThermalCaptureLevelCumulBR[i]=0; 1928 } 1929 1930 //----------------------------------------- 1931 //We calculate which transitrions go to an 1932 G4double totalThGInt=0,ENDSFLevelEnergy=0,M 1933 G4int i_closest=0,i_closest_known=0; 1934 G4double MaxAllowedLevelDistance=0.010; //1 1935 G4bool ComputePrimaryGammasEcut=false; 1936 if(PrimaryGammasEcut==0){ComputePrimaryGamm 1937 1938 //We take only those intensities going to o 1939 for(G4int i=0;i<ng;i++){ 1940 ENDSFLevelEnergy=ELevel-ThEg[i]; 1941 if(ComputePrimaryGammasEcut && PrimaryGam 1942 PrimaryGammasEcut=ENDSFLevelEnergy; 1943 } 1944 i_closest=0; 1945 MinlevelDist=1.e20; 1946 i_closest_known=0; 1947 MinlevelDist_known=1.e20; 1948 for(G4int j=0;j<NLevelsBelowThermalCaptur 1949 LevDist=std::fabs(ENDSFLevelEnergy-theL 1950 if(theLevels[j].KnownLevelID>=0){ //the 1951 if(LevDist<MinlevelDist_known){ 1952 MinlevelDist_known=LevDist; 1953 i_closest_known=j; 1954 } 1955 } 1956 if(LevDist<MinlevelDist){ 1957 MinlevelDist=LevDist; 1958 i_closest=j; 1959 } 1960 } 1961 if(MinlevelDist_known<MaxAllowedLevelDist 1962 theThermalCaptureLevelCumulBR[i_closest 1963 totalThGInt+=theThermalCaptureLevelCumu 1964 } 1965 else if(MinlevelDist<MaxAllowedLevelDista 1966 theThermalCaptureLevelCumulBR[i_closest 1967 totalThGInt+=theThermalCaptureLevelCumu 1968 } 1969 } 1970 //if(TotalThI>0){std::cout<<" NuDEX: Primar 1971 //else{std::cout<<" Primary thermal gammas 1972 std::cout<<" NuDEX: Primary thermal gammas 1973 //----------------------------------------- 1974 1975 //----------------------------------------- 1976 //If we don't have all the info, we take th 1977 if(totalThGInt<0.95){ 1978 G4double TotalNeededIntensity=1.-totalThG 1979 G4double oldInt,TotalOldIntensity=0; 1980 //------------------------- 1981 //We put the capture level replacing one 1982 Level tmpLevel; 1983 CopyLevel(&theLevels[NLevelsBelowThermalC 1984 CopyLevel(&theThermalCaptureLevel,&theLev 1985 G4double* CumulBR_th_v2=new G4double[NLev 1986 ComputeDecayIntensities(NLevelsBelowTherm 1987 CopyLevel(&tmpLevel,&theLevels[NLevelsBel 1988 //------------------------- 1989 1990 for(G4int i=0;i<NLevelsBelowThermalCaptur 1991 if(i==0){oldInt=CumulBR_th_v2[i];}else{ 1992 if(theThermalCaptureLevelCumulBR[i]==0 1993 } 1994 1995 if(TotalOldIntensity>0){ 1996 for(G4int i=0;i<NLevelsBelowThermalCapt 1997 if(theThermalCaptureLevelCumulBR[i]==0 && t 1998 if(i==0){oldInt=CumulBR_th_v2[i];}else{ol 1999 theThermalCaptureLevelCumulBR[i]=oldInt*T 2000 } 2001 } 2002 } 2003 delete [] CumulBR_th_v2; 2004 } 2005 //----------------------------------------- 2006 2007 //theThermalCaptureLevelCumulBR still not c 2008 2009 for(G4int i=1;i<NLevelsBelowThermalCaptureL 2010 theThermalCaptureLevelCumulBR[i]+=theTher 2011 } 2012 for(G4int i=0;i<NLevelsBelowThermalCaptureL 2013 theThermalCaptureLevelCumulBR[i]/=theTherm 2014 } 2015 2016 } 2017 2018 //Cambiamos las intensidades de los "primary 2019 void G4NuDEXStatisticalNucleus::ChangeThermal 2020 2021 if(!theThermalCaptureLevelCumulBR){return;} 2022 G4int level_id=GetClosestLevel(LevelEnergy, 2023 if(level_id<0 || level_id>=NLevelsBelowTher 2024 std::cout<<" ############## WARNING in "< 2025 std::cout<<" ---> "<<level_id<<" "<<Lev 2026 } 2027 2028 for(G4int i=NLevelsBelowThermalCaptureLevel 2029 theThermalCaptureLevelCumulBR[i]-=theTher 2030 } 2031 G4double OldIntensity=theThermalCaptureLeve 2032 theThermalCaptureLevelCumulBR[level_id]=abs 2033 2034 for(G4int i=1;i<NLevelsBelowThermalCaptureL 2035 theThermalCaptureLevelCumulBR[i]+=theTher 2036 } 2037 for(G4int i=0;i<NLevelsBelowThermalCaptureL 2038 theThermalCaptureLevelCumulBR[i]/=theTherm 2039 } 2040 if(level_id==0){ 2041 std::cout<<" Thermal primary gammas to le 2042 } 2043 else{ 2044 std::cout<<" Thermal primary gammas to le 2045 } 2046 } 2047 2048 void G4NuDEXStatisticalNucleus::PrintParamete 2049 2050 out<<" #################################### 2051 out<<" GENERAL_PARS"<<std::endl; 2052 out<<" Z = "<<Z_Int<<" A = "<<A_Int<<std:: 2053 out<<" Sn = "<<Sn<<" I0(ZA-1) = "<<I0<<std 2054 if(theLD!=0){theLD->PrintParameters(out);} 2055 else{out<<" No level density"<<std::endl;} 2056 out<<" PSFflag = "<<PSFflag<<std::endl; 2057 out<<" Ecrit = "<<Ecrit<<std::endl; 2058 out<<" E_unknown_min = "<<E_unk_min<<" E_u 2059 out<<" maxspin = "<<maxspinx2/2.<<std::endl 2060 out<<" MaxExcEnergy = "<<MaxExcEnergy<<std: 2061 out<<" NBands = "<<NBands<<" MinLevelsPerB 2062 out<<" Emin_bands = "<<Emin_bands<<" Emax_ 2063 out<<" NLevels = "<<NLevels<<" NKnownLeve 2064 out<<" BROpt = "<<BROpt<<" SampleGammaWid 2065 out<<" PrimaryGammasIntensityNormFactor = " 2066 out<<" KnownLevelsFlag = "<<KnownLevelsFlag 2067 out<<" ElectronConversionFlag = "<<Electron 2068 out<<" #################################### 2069 2070 } 2071 2072 void G4NuDEXStatisticalNucleus::PrintKnownLev 2073 2074 out<<" #################################### 2075 out<<" KNOWN_LEVEL_SCHEME "<<std::endl; 2076 out<<" NKnownLevels = "<<NKnownLevels<<std: 2077 char buffer[1000]; 2078 2079 //for(G4int i=0;i<NKnownLevels;i++){ 2080 for(G4int i=0;i<KnownLevelsVectorSize;i++){ 2081 snprintf(buffer,1000,"%3d %10.4g %5g %2d 2082 out<<buffer; 2083 for(G4int j=0;j<theKnownLevels[i].Ndecays 2084 snprintf(buffer,1000," %10.4g %7s",t 2085 out<<buffer; 2086 } 2087 out<<std::endl; 2088 for(G4int j=0;j<theKnownLevels[i].NGammas 2089 snprintf(buffer,1000," 2090 out<<buffer<<std::endl; 2091 } 2092 } 2093 out<<" #################################### 2094 2095 } 2096 2097 void G4NuDEXStatisticalNucleus::PrintKnownLev 2098 2099 out<<" #################################### 2100 out<<" KNOWN_LEVES_DEGEN "<<std::endl; 2101 out<<" NKnownLevels = "<<NKnownLevels<<std: 2102 char buffer[1000]; 2103 2104 for(G4int i=0;i<NKnownLevels;i++){ 2105 G4double MaxIntens=-100; 2106 G4double GammaEnergy; 2107 for(G4int j=0;j<theKnownLevels[i].NGammas 2108 if(theKnownLevels[i].Pg[j]>MaxIntens){M 2109 } 2110 for(G4int j=0;j<theKnownLevels[i].NGammas 2111 //snprintf(buffer,1000,"%10.3f %9.3f %9 2112 GammaEnergy=theKnownLevels[i].Energy-th 2113 snprintf(buffer,1000,"%10.3f %9.3f %9.3 2114 out<<buffer<<std::endl; 2115 } 2116 } 2117 out<<" #################################### 2118 2119 } 2120 2121 void G4NuDEXStatisticalNucleus::PrintLevelDen 2122 2123 if(theLD==0){return;} 2124 2125 G4double Emin=0; 2126 G4double Emax=E_unk_max; 2127 G4int np=100; 2128 2129 out<<" #################################### 2130 out<<" LEVELDENSITY"<<std::endl; 2131 G4double ene,exp=0; 2132 G4double *ld=new G4double[maxspinx2+2]; 2133 G4bool *WriteThisSpin=new G4bool[maxspinx2+ 2134 2135 for(G4int spx2=0;spx2<=maxspinx2;spx2++){ 2136 WriteThisSpin[spx2]=true; 2137 if(((A_Int+spx2)%2)!=0){ 2138 WriteThisSpin[spx2]=false; 2139 } 2140 } 2141 2142 out<<np<<" "<<Emin<<" "<<Emax<<" "<<Ecri 2143 out<<"ENE EXP TOT SUM(J)"; 2144 for(G4int spx2=0;spx2<=maxspinx2;spx2++){ 2145 if(WriteThisSpin[spx2]){out<<" J="<<spx 2146 } 2147 out<<std::endl; 2148 2149 for(G4int i=0;i<np;i++){ 2150 ene=Emin+(Emax-Emin)*i/(G4double)(np-1); 2151 exp=0; 2152 for(G4int j=0;j<NLevels;j++){if(theLevels 2153 out<<ene<<" "<<exp<<" "; 2154 ld[maxspinx2+1]=0; 2155 for(G4int spx2=0;spx2<=maxspinx2;spx2++){ 2156 ld[spx2]=2*theLD->GetLevelDensity(ene,s 2157 ld[maxspinx2+1]+=ld[spx2]; 2158 } 2159 //out<<ld[maxspinx2+1]; 2160 out<<theLD->GetLevelDensity(ene,0,true,tr 2161 for(G4int spx2=0;spx2<=maxspinx2;spx2++){ 2162 if(WriteThisSpin[spx2]){out<<" "<<ld[ 2163 } 2164 out<<std::endl; 2165 } 2166 out<<" #################################### 2167 2168 delete [] ld; 2169 delete [] WriteThisSpin; 2170 } 2171 2172 void G4NuDEXStatisticalNucleus::PrintLevelSch 2173 2174 std::ofstream out(fname); 2175 char buffer[1000]; 2176 for(G4int i=0;i<NLevels;i++){ 2177 if(theLevels[i].Energy>Ecrit && (MaxLevel 2178 snprintf(buffer,1000,"%13.5f %17.8f %17 2179 out<<buffer<<std::endl; 2180 } 2181 } 2182 out.close(); 2183 2184 } 2185 2186 void G4NuDEXStatisticalNucleus::PrintLevelSch 2187 out<<" #################################### 2188 out<<" LEVELSCHEME"<<std::endl; 2189 for(G4int i=0;i<NLevels;i++){ 2190 out<<i<<" "<<theLevels[i].Energy<<" "<< 2191 } 2192 out<<" #################################### 2193 } 2194 2195 void G4NuDEXStatisticalNucleus::PrintThermalP 2196 2197 out<<" #################################### 2198 out<<" THERMAL PRIMARY TRANSITIONS"<<std::e 2199 out<<" "<<NLevelsBelowThermalCaptureLevel<< 2200 if(theThermalCaptureLevelCumulBR!=0){ 2201 out<<" "<<0<<" "<<theLevels[0].Energy<<" 2202 for(G4int i=1;i<NLevelsBelowThermalCaptur 2203 out<<" "<<i<<" "<<theLevels[i].Energy< 2204 } 2205 } 2206 out<<" #################################### 2207 2208 G4double ThresholdIntensity=0.01; 2209 out<<" #################################### 2210 out<<" STRONGEST THERMAL PRIMARY TRANSITION 2211 out<<" "<<NLevelsBelowThermalCaptureLevel<< 2212 if(theThermalCaptureLevelCumulBR!=0){ 2213 if(theThermalCaptureLevelCumulBR[0]>Thres 2214 for(G4int i=1;i<NLevelsBelowThermalCaptur 2215 if(theThermalCaptureLevelCumulBR[i]-the 2216 } 2217 } 2218 out<<" #################################### 2219 } 2220 2221 void G4NuDEXStatisticalNucleus::PrintTotalCum 2222 2223 if(TotalCumulBR[i_level]!=0){ 2224 out<<" ################################## 2225 out<<" CUMULBR FROM LEVEL "<<i_level<<" w 2226 for(G4int i=0;i<i_level;i++){ 2227 out<<theLevels[i].Energy<<" "<<theLeve 2228 } 2229 out<<" ################################## 2230 } 2231 2232 } 2233 2234 void G4NuDEXStatisticalNucleus::PrintBR(G4int 2235 2236 if(TotalCumulBR[i_level]!=0){ 2237 out<<" ################################## 2238 out<<" BR FROM LEVEL "<<i_level<<" with E 2239 for(G4int i=0;i<i_level;i++){ 2240 if(theLevels[i].Energy<MaxExcEneToPrint 2241 if(i==0){ 2242 out<<theLevels[i].Energy<<" "<<theLevels 2243 } 2244 else{ 2245 out<<theLevels[i].Energy<<" "<<theLevels 2246 } 2247 } 2248 } 2249 out<<" ################################## 2250 } 2251 2252 2253 } 2254 2255 2256 void G4NuDEXStatisticalNucleus::PrintPSF(std: 2257 2258 thePSF->PrintPSFParameters(out); 2259 2260 G4int NVals=400; 2261 G4int nEnePSF=(G4int)Sn+1; //number of exci 2262 G4double EnePSF[200]; 2263 G4double Emin=0; 2264 G4double Emax=10; 2265 G4double xval,e1,m1,e2; 2266 2267 out<<" #################################### 2268 out<<" PSF"<<std::endl; 2269 out<<" "<<NVals<<" "<<Emin<<" "<<Emax<<" 2270 EnePSF[0]=Sn; 2271 for(G4int i=1;i<nEnePSF;i++){ 2272 EnePSF[i]=i; 2273 } 2274 for(G4int i=0;i<nEnePSF;i++){ 2275 out<<" "<<EnePSF[i]; 2276 } 2277 out<<std::endl; 2278 char word[1000]; 2279 out<<" E E1 M1 E2 2280 for(G4int i=0;i<nEnePSF;i++){ 2281 for(G4int j=0;j<NVals;j++){ 2282 xval=Emin+(Emax-Emin)*j/(NVals-1.); 2283 if(xval==0){xval=1.e-6;} 2284 e1=thePSF->GetE1(xval,EnePSF[i]); 2285 m1=thePSF->GetM1(xval,EnePSF[i]); 2286 e2=thePSF->GetE2(xval,EnePSF[i]); 2287 snprintf(word,1000," %10.4E %10.4E %10. 2288 out<<word<<std::endl; 2289 } 2290 } 2291 out<<" #################################### 2292 2293 } 2294 2295 2296 void G4NuDEXStatisticalNucleus::PrintICC(std: 2297 2298 theICC->PrintICC(out); 2299 2300 } 2301 2302 void G4NuDEXStatisticalNucleus::PrintAll(std: 2303 2304 PrintParameters(out); 2305 PrintKnownLevels(out); 2306 PrintLevelDensity(out); 2307 PrintLevelScheme(out); 2308 PrintThermalPrimaryTransitions(out); 2309 PrintPSF(out); 2310 PrintICC(out); 2311 2312 } 2313 2314 2315 2316 void G4NuDEXStatisticalNucleus::PrintInput01( 2317 2318 out<<"LEVELDENSITYTYPE "<<LevelDensityType< 2319 out<<"MAXSPIN "<<maxspinx2/2.<<std::endl; 2320 out<<"MINLEVELSPERBAND "<<MinLevelsPerBand< 2321 out<<"BANDWIDTH_MEV "<<BandWidth<<std::endl 2322 out<<"MAXEXCENERGY_MEV "<<MaxExcEnergy<<std 2323 out<<"ECRIT_MEV "<<Ecrit<<std::endl; 2324 out<<"KNOWNLEVELSFLAG "<<KnownLevelsFlag<<s 2325 out<<std::endl; 2326 out<<"PSF_FLAG "<<PSFflag<<std::endl; 2327 out<<"BROPTION "<<BROpt<<std::endl; 2328 out<<"SAMPLEGAMMAWIDTHS "<<SampleGammaWidth 2329 out<<std::endl; 2330 out<<"SEED1 "<<seed1<<std::endl; 2331 out<<"SEED2 "<<seed2<<std::endl; 2332 out<<"SEED3 "<<seed3<<std::endl; 2333 out<<std::endl; 2334 out<<"ELECTRONCONVERSIONFLAG "<<ElectronCon 2335 out<<"PRIMARYTHCAPGAMNORM "<<PrimaryGammasI 2336 out<<"PRIMARYGAMMASECUT "<<PrimaryGammasEcu 2337 out<<std::endl; 2338 theLD->PrintParametersInInputFileFormat(out 2339 thePSF->PrintPSFParametersInInputFileFormat 2340 out<<std::endl; 2341 out<<"END"<<std::endl; 2342 2343 } 2344 2345 G4int ComparisonLevels(const void* va, const 2346 { 2347 Level* a, *b; 2348 a = (Level*) va; 2349 b = (Level*) vb; 2350 if( a->Energy == b->Energy ) return 0; 2351 return( ( a->Energy ) > ( b->Energy ) ) ? 1: 2352 } 2353 2354 void CopyLevel(Level* a,Level* b){ 2355 b->Energy=a->Energy; 2356 b->spinx2=a->spinx2; 2357 b->parity=a->parity; 2358 b->seed=a->seed; 2359 b->KnownLevelID=a->KnownLevelID; 2360 b->NLevels=a->NLevels; 2361 b->Width=a->Width; 2362 } 2363 2364 2365 void CopyLevel(KnownLevel* a,Level* b){ 2366 b->Energy=a->Energy; 2367 b->spinx2=a->spinx2; 2368 b->parity=a->parity; 2369 b->seed=0; 2370 b->KnownLevelID=-1; 2371 b->NLevels=1; 2372 b->Width=0; 2373 } 2374 2375 2376 2377 2378