| 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