| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // Author: E.Mendoza
30 //
31 // Creation date: May 2024
32 //
33 // Modifications:
34 //
35 // -------------------------------------------------------------------
36 //
37 // NuDEX code (https://doi.org/10.1016/j.nima.2022.167894)
38 //
39
40
41
42
43 #include "G4NuDEXStatisticalNucleus.hh"
44 #include "G4NuDEXLevelDensity.hh"
45 #include "G4NuDEXInternalConversion.hh"
46 #include "G4NuDEXPSF.hh"
47
48
49
50 G4NuDEXStatisticalNucleus::G4NuDEXStatisticalNucleus(G4int Z,G4int A){
51
52 //The default values for these flags are in "GeneralStatNuclParameters.dat"
53 //Can be changed with G4NuDEXStatisticalNucleus::SetSomeInitalParameters(...)
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 G4NuDEXStatisticalNucleus::Init(...)
64 //Can be changed via G4NuDEXStatisticalNucleus::SetInitialParameters02(...):
65 ElectronConversionFlag=-1; // All EC
66 KnownLevelsFlag=-1; //Use all known levels
67 PrimaryGammasIntensityNormFactor=-1;
68 PrimaryGammasEcut=-1; //If primary gammas, do not create new primary gammas going to the "Primary gammas" region.
69 Ecrit=-1;
70
71 hasBeenInitialized=false;
72 NBands=-1;
73 theLevels=0;
74 theKnownLevels=0;
75 NKnownLevels=0; NUnknownLevels=0; NLevels=0; KnownLevelsVectorSize=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=false; Rand3seedProvided=false;
97 }
98
99 void G4NuDEXStatisticalNucleus::SetInitialParameters02(G4int knownLevelsFlag,G4int electronConversionFlag,G4double primGamNormFactor,G4double primGamEcut,G4double ecrit){
100 if(hasBeenInitialized){
101 //std::cout<<" ############## Error: G4NuDEXStatisticalNucleus::SetInitialParameters02 cannot be used after initializing the nucleus ##############"<<std::endl;
102 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
103 }
104 if(knownLevelsFlag>=0){KnownLevelsFlag=knownLevelsFlag;}
105 if(electronConversionFlag>=0){ElectronConversionFlag=electronConversionFlag;}
106 if(primGamNormFactor>=0){PrimaryGammasIntensityNormFactor=primGamNormFactor;}
107 if(primGamEcut>=0){PrimaryGammasEcut=primGamEcut;}
108 if(ecrit>=0){Ecrit=ecrit;}
109
110 }
111
112 void G4NuDEXStatisticalNucleus::SetSomeInitalParameters(G4int LDtype,G4int PSFFlag,G4double MaxSpin,G4int minlevelsperband,G4double BandWidth_MeV,G4double maxExcEnergy,G4int BrOption,G4int sampleGammaWidths,unsigned int aseed1,unsigned int aseed2,unsigned int aseed3){
113
114 if(hasBeenInitialized){
115 std::cout<<" ############## Error: G4NuDEXStatisticalNucleus::SetSomeInitalParameters cannot be used after initializing the nucleus ##############"<<std::endl;
116 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
117 }
118
119 if(LDtype>0){LevelDensityType=LDtype;}
120 if(PSFFlag>=0){PSFflag=PSFFlag;}
121 if(MaxSpin>0){maxspinx2=(G4int)(2.*MaxSpin+0.01);}
122 if(minlevelsperband>0){MinLevelsPerBand=minlevelsperband;}
123 if(BandWidth_MeV!=0){BandWidth=BandWidth_MeV;}
124 if(maxExcEnergy!=0){MaxExcEnergy=maxExcEnergy;}
125 if(BrOption>0){BROpt=BrOption;}
126 if(sampleGammaWidths>=0){SampleGammaWidths=sampleGammaWidths;}
127 if(aseed1>0){seed1=aseed1; theRandom1->SetSeed(seed1); Rand1seedProvided=true;}
128 if(aseed2>0){seed2=aseed2; theRandom2->SetSeed(seed2); Rand2seedProvided=true;}
129 if(aseed3>0){seed3=aseed3; theRandom3->SetSeed(seed3); Rand3seedProvided=true;}
130
131 }
132
133
134
135 G4NuDEXStatisticalNucleus::~G4NuDEXStatisticalNucleus(){
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].decayFraction;
141 delete [] theKnownLevels[i].decayMode;
142 }
143 if(theKnownLevels[i].NGammas>0){
144 delete [] theKnownLevels[i].FinalLevelID;
145 delete [] theKnownLevels[i].multipolarity;
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 [] theKnownLevels;}
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 [] theThermalCaptureLevelCumulBR;}
162 if(TotalCumulBR!=0){
163 for(G4int i=0;i<NLevels;i++){
164 if(TotalCumulBR[i]!=0){delete [] TotalCumulBR[i];}
165 }
166 delete [] TotalCumulBR;
167 }
168 }
169
170
171 G4int G4NuDEXStatisticalNucleus::Init(const char* dirname,const char* inputfname){
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/SpecialInputs/ZA_%d.dat",dirname,Z_Int*1000+A_Int);
182 G4int HasDefaultInput=ReadSpecialInputFile(defaultinputfname);
183 char* definputfn=0;
184 if(HasDefaultInput>0){definputfn=defaultinputfname;}
185
186 //General statistical parameters:
187 snprintf(fname,1000,"%s/GeneralStatNuclParameters.dat",dirname);
188 check=ReadGeneralStatNuclParameters(fname); if(check<0){return -1;}
189
190 //Some default, if not initialized yet:
191 if(ElectronConversionFlag<0){ElectronConversionFlag=2;} // All EC
192 if(KnownLevelsFlag<0){KnownLevelsFlag=1;} //Use all known levels
193 if(PrimaryGammasIntensityNormFactor<0){PrimaryGammasIntensityNormFactor=1;}
194 if(PrimaryGammasEcut<0){PrimaryGammasEcut=0;}
195 if(Ecrit<0){
196 snprintf(fname,1000,"%s/KnownLevels/levels-param.data",dirname);
197 check=ReadEcrit(fname); if(check<0){return -1;}
198 }
199
200
201 //Level density:
202 theLD=new G4NuDEXLevelDensity(Z_Int,A_Int,LevelDensityType);
203 check=theLD->ReadLDParameters(dirname,inputfname,definputfn); //if(check<0){return -1;}
204 LevelDensityType=theLD->GetLDType(); //because it can be changed by inputfname or due to lack of data
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.dat",dirname,Z_Int);
215 check=ReadKnownLevels(fname); if(check<0){return -1;} //here we get/crosscheck Sn
216 I0=TakeTargetNucleiI0(fname,check); if(check<0){return -1;} //if no I0 --> out
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 known level scheme is not complete, then we can do nothing ...
228 if(theLD==0 && Ecrit<MaxExcEnergy){
229 std::cout<<" ###### WARNING: No level density and level scheme not complete for ZA="<<1000*Z_Int+A_Int<<" --> Ecrit="<<Ecrit<<" MeV and MaxExcEnergy = "<<MaxExcEnergy<<" MeV ######"<<std::endl;
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 some bands
241 NBands=0;
242 while(E_unk_min+BandWidth*NBands<MaxExcEnergy){
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 ..."<<std::endl;
258 CreateLevelScheme();
259 //std::cout<<" ............. done"<<std::endl;
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->Integer(4294967295)+1;
268 }
269
270 //Internal conversion:
271 theICC=new G4NuDEXInternalConversion(Z_Int);
272 snprintf(fname,1000,"%s/ICC_factors.dat",dirname);
273 theICC->Init(fname);
274 theICC->SetRandom4Seed(theRandom3->GetSeed()); //same seed as for generating the cascades
275
276 //PSF:
277 thePSF=new G4NuDEXPSF(Z_Int,A_Int);
278 thePSF->Init(dirname,theLD,inputfname,definputfn,PSFflag);
279
280 //We compute the missing BR in the known part of the level scheme:
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::MakeSomeParameterChecks01(){
308
309 if(LevelDensityType<1 || LevelDensityType>3){
310 std::cout<<" ############## Error, LevelDensityType cannot be set to: "<<LevelDensityType<<" ##############"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
311 }
312
313 if(maxspinx2<=0){
314 std::cout<<" ############## Error, MaxSpin cannot be set to: "<<maxspinx2/2.<<" ##############"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
315 }
316
317 if(MaxExcEnergy<=0){
318 std::cout<<" ############## Error, MaxExcEnergy cannot be set to: "<<MaxExcEnergy<<" ##############"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
319 }
320
321 if(BROpt<0 || BROpt>2){
322 std::cout<<" ############## Error, BROpt cannot be set to: "<<BROpt<<" ##############"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
323 }
324
325 if(SampleGammaWidths<0 || SampleGammaWidths>1){
326 std::cout<<" ############## Error, SampleGammaWidths cannot be set to: "<<SampleGammaWidths<<" ##############"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
327 }
328
329 if(KnownLevelsFlag!=0 && KnownLevelsFlag!=1){
330 std::cout<<" ############## Error, KnownLevelsFlag cannot be set to: "<<KnownLevelsFlag<<" ##############"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
331 }
332
333 if(ElectronConversionFlag!=0 && ElectronConversionFlag!=1 && ElectronConversionFlag!=2){
334 std::cout<<" ############## Error, ElectronConversionFlag cannot be set to: "<<ElectronConversionFlag<<" ##############"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
335 }
336
337 if(PrimaryGammasIntensityNormFactor<=0){
338 std::cout<<" ############## Error, PrimaryGammasIntensityNormFactor cannot be set to: "<<PrimaryGammasIntensityNormFactor<<" ##############"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
339 }
340
341 if(PrimaryGammasEcut<0){
342 std::cout<<" ############## Error, PrimaryGammasEcut cannot be set to: "<<PrimaryGammasEcut<<" ##############"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
343 }
344
345 if(Ecrit<0){
346 std::cout<<" ############## Error, Ecrit cannot be set to: "<<Ecrit<<" ##############"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
347 }
348
349 }
350
351
352 //If InitialLevel==-1 then we start from the thermal capture level
353 //If ExcitationEnergy>0 then is the excitation energy of the nucleus
354 //If ExcitationEnergy<0 then is a capture reaction of a neutron with energy -ExcitationEnergy
355 // return Npar (number of particles emitted). If something goes wrong, returns negative value (for example negative energy transition, which could happen).
356 G4int G4NuDEXStatisticalNucleus::GenerateCascade(G4int InitialLevel,G4double ExcitationEnergy,std::vector<char>& pType,std::vector<double>& pEnergy,std::vector<double>& pTime){
357
358 pType.clear();
359 pEnergy.clear();
360 pTime.clear();
361
362 if(ExcitationEnergy<0){
363 ExcitationEnergy=Sn-(A_Int-1.)/(G4double)A_Int*ExcitationEnergy;
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; //icc factor, energy of the transition, initial/final excitation energy
372 G4double EmissionTime=0; //in seconds
373 G4int NTransition=0;
374 //G4double TotalCascadeEnergy1=0,TotalCascadeEnergy2=0;
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, there are no thermal capture for ZA="<<A_Int+1000*Z_Int-1<<" , with Sn = "<<Sn<<" ##############"<<std::endl;
398 }
399 else{
400 //Sample final level:
401 G4double randnumber=theRandom3->Uniform();
402 f_level=-1;
403 for(G4int i=0;i<NLevelsBelowThermalCaptureLevel;i++){
404 if(theThermalCaptureLevelCumulBR[i]>randnumber){
405 multipol=GetMultipolarity(&theThermalCaptureLevel,&theLevels[i]);
406 f_level=i;
407 break;
408 }
409 }
410 }
411 if(f_level<0){
412 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
413 }
414 Exc_ene_f=theLevels[f_level].Energy;
415 }
416 else if(i_level>0){
417 f_level=SampleFinalLevel(i_level,multipol,alpha,NTransition);
418 }
419 else{
420 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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 band of levels:
428 if(theLevels[f_level].Width!=0){
429 Exc_ene_f+=theRandom3->Uniform(-theLevels[f_level].Width,+theLevels[f_level].Width);
430 }
431 E_trans=Exc_ene_i-Exc_ene_f;
432 if(E_trans<=0){
433 //std::cout<<"Exc_ene_i = "<<Exc_ene_i<<" Exc_ene_f = "<<Exc_ene_f<<std::endl;
434 //std::cout<<" ####### WARNING: E_trans = "<<E_trans<<" for i="<<i_level<<" with E = "<<theLevels[std::max(i_level,0)].Energy<<" to f="<<f_level<<" with E = "<<theLevels[f_level].Energy<<" ########"<<std::endl;
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[i_level].T12/std::log(2));
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){ //ElectronConversionFlag=1,2
451 ele_conv=theICC->SampleInternalConversion(E_trans,multipol,alpha); //use the alpha value from the know level value
452 }
453 else if(ElectronConversionFlag==2){
454 ele_conv=theICC->SampleInternalConversion(E_trans,multipol); //calculate alpha (icc factor)
455 }
456 }
457 //------------------------------------------------------------
458 //std::cout<<" ---- "<<Exc_ene_i<<" "<<Exc_ene_f<<" "<<E_trans<<" "<<multipol<<" "<<ele_conv<<std::endl;
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++){TotalCascadeEnergy2+=pEnergy[i];}
489 //std::cout<<" Total energy: "<<TotalCascadeEnergy1<<" "<<TotalCascadeEnergy2<<std::endl; getchar();
490
491
492 if(i_level!=0){
493 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
494 }
495 return Npar;
496 }
497
498
499
500
501 G4int G4NuDEXStatisticalNucleus::SampleFinalLevel(G4int i_level,G4int& multipolarity,G4double &icc_fac,G4int nTransition){
502
503 if(i_level<=0 || i_level>=NLevels){
504 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
505 }
506
507 G4double randnumber=theRandom3->Uniform();
508
509 G4int i_knownLevel=-1;
510 if(i_level<NKnownLevels){ //then is a known level
511 i_knownLevel=i_level;
512 }
513 if(theLevels[i_level].KnownLevelID>0){ //then is in the unknown part, but we use it as a known level
514 if(theKnownLevels[theLevels[i_level].KnownLevelID].NGammas>0){
515 i_knownLevel=theLevels[i_level].KnownLevelID;
516 }
517 }
518
519 if(i_knownLevel>=0){//known part of the spectrum
520 theSampledLevel=-1;
521 for(G4int j=0;j<theKnownLevels[i_knownLevel].NGammas;j++){
522 if(theKnownLevels[i_knownLevel].cumulPtot[j]>randnumber){
523 multipolarity=theKnownLevels[i_knownLevel].multipolarity[j];
524 icc_fac=theKnownLevels[i_knownLevel].Icc[j];
525 return theKnownLevels[i_knownLevel].FinalLevelID[j];
526 }
527 }
528 std::cout<<randnumber<<" "<<i_knownLevel<<" "<<theKnownLevels[i_knownLevel].NGammas<<std::endl;
529 for(G4int j=0;j<theKnownLevels[i_knownLevel].NGammas;j++){
530 std::cout<<theKnownLevels[i_knownLevel].cumulPtot[j]<<std::endl;
531 }
532 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
533 }
534 else{
535 icc_fac=-1;
536 //------------------------------------------------------------------------------
537 //If BROpt==1 or 2, then we store the BR, if not computed, or calculate the final level from it
538 if(BROpt==1 || (BROpt==2 && nTransition==1)){
539 //maybe the TotalGammaRho[i_level] and BR have not been computed yet:
540 if(TotalCumulBR[i_level]==0){
541 TotalCumulBR[i_level]=new G4double[i_level];
542 TotalGammaRho[i_level]=ComputeDecayIntensities(i_level,TotalCumulBR[i_level]);
543 }
544 for(G4int j=0;j<i_level;j++){
545 if(TotalCumulBR[i_level][j]>randnumber){
546 multipolarity=GetMultipolarity(&theLevels[i_level],&theLevels[j]);
547 return j;
548 }
549 }
550 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
551 }
552 //------------------------------------------------------------------------------
553
554 //BROpt==0
555 //------------------------------------------------------------------------------
556 // If not, maybe the TotalGammaRho[i_level] has not been computed yet:
557 if(TotalGammaRho[i_level]<0){//not computed, we compute it:
558 TotalGammaRho[i_level]=ComputeDecayIntensities(i_level);
559 }
560 theSampledLevel=-1;
561 ComputeDecayIntensities(i_level,0,randnumber); // here we compute the final level
562 multipolarity=theSampledMultipolarity;
563 return theSampledLevel;
564 //------------------------------------------------------------------------------
565 }
566
567 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
568 return 0;
569 }
570
571 void G4NuDEXStatisticalNucleus::ChangeLevelSpinParityAndBR(G4int i_level,G4int newspinx2,G4bool newParity,G4int nlevels,G4double width,unsigned int seed){
572
573 if(i_level==-1){ //change BR of thermal, ignore arguments
574 if(Sn>0 && NLevels>1){
575 CreateThermalCaptureLevel(seed);
576 GenerateThermalCaptureLevelBR(theLibDir.c_str());
577 }
578 return;
579 }
580
581 if(i_level<0 || i_level>=NLevels){
582 std::cout<<" i_level = "<<i_level<<" ------ NLevels = "<<NLevels<<std::endl;
583 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
584 }
585
586 //Do not apply to known levels:
587 if(i_level<NKnownLevels || theLevels[i_level].KnownLevelID>0){
588 std::cout<<" ####### WARNING: you are trying to change the BR, spin, parity, etc. of a known level --> nothing is done ############"<<std::endl;
589 return;
590 //NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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->Integer(4294967295)+1;
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 have to change TotalGammaRho[i_level]
609 G4double* br_vector=0;
610 if(TotalCumulBR!=0){
611 br_vector=TotalCumulBR[i_level];
612 }
613 TotalGammaRho[i_level]=ComputeDecayIntensities(i_level,br_vector);
614 }
615
616 }
617
618
619 //if randnumber<0, return the total TotalGammaRho, and if cumulativeBR!=0, the corresponding cumulativeBR vector is calculated
620 //if randnumber>0, it is assumed that TotalGammaRho has been already computed
621 // (in the TotalGammaRho[] array or in the TotGR argument) and is used to sample the transition
622 // The result is stored in theSampledLevel and theSampledMultipolarity variables
623 G4double G4NuDEXStatisticalNucleus::ComputeDecayIntensities(G4int i_level,G4double* cumulativeBR,G4double randnumber,G4double TotGR,G4bool AllowE1){
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_level<<" out of range. NLevels = "<<NLevels<<std::endl;
631 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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(__LINE__).c_str(),"##### Error in NuDEX #####");
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_level].Width<=theLevels[j].Energy+theLevels[j].Width){
649 thisTotalGammaRho+=0; //not cecessary, but for understanding ...
650 if(ComputeAlsoBR){
651 cumulativeBR[j]=0;
652 }
653 }
654 else{
655 G4double Eg=theLevels[i_level].Energy-theLevels[j].Energy;
656
657 //------------------------------------------------------------------
658 //Check which are allowed transitions:
659 G4bool E1allowed=true,M1allowed=true,E2allowed=true;
660 G4int Lmin=std::abs(theLevels[i_level].spinx2-theLevels[j].spinx2)/2;
661 G4int Lmax=(theLevels[i_level].spinx2+theLevels[j].spinx2)/2;
662 if(theLevels[i_level].parity==theLevels[j].parity){
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].NLevels*theLevels[j].NLevels;
680
681 G4double rand;
682 G4double MaxNSamplesForChi2=1000;
683
684 if(E1allowed){
685 Sumrand2=RealNTransitions;
686 if(SampleGammaWidths==1){ //Porter-Thomas fluctuations
687 Sumrand2=0;
688 if(RealNTransitions>MaxNSamplesForChi2){
689 Sumrand2=RealNTransitions*theRandom2->Gaus(1,std::sqrt(2./RealNTransitions));
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,theLevels[i_level].Energy);
699 if(randnumber>=0){
700 if(thisTotalGammaRho+GammaRho>=TotGR*randnumber){
701 theSampledMultipolarity=1;
702 }
703 }
704 }
705 if(M1allowed){
706 Sumrand2=RealNTransitions;
707 if(SampleGammaWidths==1){ //Porter-Thomas fluctuations
708 Sumrand2=0;
709 if(RealNTransitions>MaxNSamplesForChi2){
710 Sumrand2=RealNTransitions*theRandom2->Gaus(1,std::sqrt(2./RealNTransitions));
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,theLevels[i_level].Energy);
720 if(randnumber>=0){
721 if(thisTotalGammaRho+GammaRho>=TotGR*randnumber){
722 theSampledMultipolarity=-1;
723 }
724 }
725 }
726 if(E2allowed){
727 Sumrand2=RealNTransitions;
728 if(SampleGammaWidths==1){ //Porter-Thomas fluctuations
729 Sumrand2=0;
730 if(RealNTransitions>MaxNSamplesForChi2){
731 Sumrand2=RealNTransitions*theRandom2->Gaus(1,std::sqrt(2./RealNTransitions));
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->GetE2(Eg,theLevels[i_level].Energy);
741 if(randnumber>=0){
742 if(thisTotalGammaRho+GammaRho>=TotGR*randnumber && theSampledMultipolarity<-10){
743 theSampledMultipolarity=2;
744 }
745 }
746 }
747
748 /*
749 if(i_level==NLevels-1 && j<10){
750 std::cout<<j<<" "<<GammaRho<<" "<<E1allowed<<" "<<M1allowed<<" "<<E2allowed<<" "<<GammaRho/Eg/Eg/Eg/thePSF->GetE1(Eg,theLevels[i_level].Energy)<<std::endl;
751 }
752 */
753
754 thisTotalGammaRho+=GammaRho;
755 if(ComputeAlsoBR){
756 cumulativeBR[j]=GammaRho;
757 }
758 }
759
760 if(randnumber>=0 && thisTotalGammaRho>=TotGR*randnumber){
761 if(theSampledMultipolarity==-50){
762 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
763 }
764 theSampledLevel=j;
765 return -1;
766 }
767 }
768
769 //If there are no allowed transitions:
770 if(randnumber>=0 && thisTotalGammaRho==0){ //if randnumber>0 then TotalGammaRho[i_lev] has been already computed allowing E1 transitions
771 return ComputeDecayIntensities(i_level,0,randnumber,TotGR,true);
772 }
773
774 if(randnumber>=0){ //then we should not be here
775 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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)>1.e-10){
787 std::cout<<" ############### Warning: cumulativeBR["<<i_level<<"]["<<i_level-1<<"]-1 is "<<cumulativeBR[i_level-1]-1<<" ###############"<<std::endl;
788 }
789 }
790 }
791 else{
792 //std::cout<<" ############### WARNING: total GammaRho for level "<<i_level<<" is "<<thisTotalGammaRho<<". We recalculate it allowing all transitions and assuming the E1 PSF ###############"<<std::endl;
793 //NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
794 thisTotalGammaRho=ComputeDecayIntensities(i_level,cumulativeBR,-1,-1,true);
795 }
796
797 return thisTotalGammaRho;
798 }
799
800
801 //retrieves the "lowest" allowed multipolarity:
802 G4int G4NuDEXStatisticalNucleus::GetMultipolarity(Level* theInitialLevel,Level* theFinalLevel){
803
804 if(theInitialLevel->spinx2+theFinalLevel->spinx2==0){
805 return 0;
806 }
807 G4int Lmin=std::abs(theInitialLevel->spinx2-theFinalLevel->spinx2)/2;
808 if(Lmin==0){Lmin=1;}
809 if(Lmin%2==0){
810 if(theInitialLevel->parity==theFinalLevel->parity){
811 return +Lmin;
812 }
813 else{
814 return -Lmin;
815 }
816 }
817 else{
818 if(theInitialLevel->parity==theFinalLevel->parity){
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 spin and parity
831 //if spinx2<0, then retrieves the closest level
832 G4int G4NuDEXStatisticalNucleus::GetClosestLevel(G4double Energy,G4int spinx2,G4bool parity){
833
834 //std::cout<<" XXX finding closest level of spin "<<spinx2/2.<<" and parity "<<parity<<" to "<<Energy<<" MeV"<<std::endl;
835
836 //------------------------------------------------------------------------------
837 // We try to go closer to the solution, otherwise it takes too much time:
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 && theLevels[i].parity==parity) || spinx2<0){
854 break;
855 }
856 }
857 for(G4int i=i_close_down;i>=0;i--){
858 i_down=i;
859 if((theLevels[i].spinx2==spinx2 && theLevels[i].parity==parity) || spinx2<0){
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].Energy-Energy);
869 if((theLevels[i].spinx2==spinx2 && theLevels[i].parity==parity) || spinx2<0){ //then this is a candidate
870 if(EnergyDistance<MinEnergyDistance || MinEnergyDistance<0){
871 MinEnergyDistance=EnergyDistance;
872 result=i;
873 }
874 }
875 }
876 //std::cout<<" XXX found --> "<<result<<std::endl;
877
878
879 return result;
880 }
881
882
883 Level* G4NuDEXStatisticalNucleus::GetLevel(G4int i_level){
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="<<A_Int+1000*Z_Int<<" , requested level i_level="<<i_level<<" does not exist ############"<<std::endl;
893
894 return 0;
895 }
896
897 G4double G4NuDEXStatisticalNucleus::GetLevelEnergy(G4int i_level){
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::ComputeKnownLevelsMissingBR(){
908
909
910 for(G4int i=1;i<NKnownLevels;i++){
911 if(theKnownLevels[i].NGammas!=0){
912 for(G4int j=0;j<theKnownLevels[i].NGammas;j++){
913 theKnownLevels[i].multipolarity[j]=GetMultipolarity(&theLevels[i],&theLevels[theKnownLevels[i].FinalLevelID[j]]);
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].NGammas++;}
921 }
922 if(tmp[0]!=0){theKnownLevels[i].NGammas++;}
923 if(theKnownLevels[i].NGammas<=0){
924 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
925 }
926
927 theKnownLevels[i].FinalLevelID=new G4int[theKnownLevels[i].NGammas];
928 theKnownLevels[i].multipolarity=new G4int[theKnownLevels[i].NGammas];
929 theKnownLevels[i].Eg=new G4double[theKnownLevels[i].NGammas];
930 theKnownLevels[i].cumulPtot=new G4double[theKnownLevels[i].NGammas];
931 theKnownLevels[i].Pg=new G4double[theKnownLevels[i].NGammas];
932 theKnownLevels[i].Pe=new G4double[theKnownLevels[i].NGammas];
933 theKnownLevels[i].Icc=new G4double[theKnownLevels[i].NGammas];
934 G4int cont=0;
935 if(tmp[0]!=0){
936 theKnownLevels[i].FinalLevelID[cont]=0;
937 theKnownLevels[i].Eg[cont]=theKnownLevels[i].Energy-theKnownLevels[0].Energy;
938 theKnownLevels[i].cumulPtot[cont]=tmp[0];
939 theKnownLevels[i].multipolarity[cont]=GetMultipolarity(&theLevels[i],&theLevels[theKnownLevels[i].FinalLevelID[cont]]);
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[i].Energy-theKnownLevels[j].Energy;
946 theKnownLevels[i].cumulPtot[cont]=tmp[j];
947 theKnownLevels[i].multipolarity[cont]=GetMultipolarity(&theLevels[i],&theLevels[theKnownLevels[i].FinalLevelID[cont]]);
948 cont++;
949 }
950 }
951 delete [] tmp;
952 for(G4int j=0;j<theKnownLevels[i].NGammas;j++){
953 theKnownLevels[i].Icc[j]=0;
954 if(ElectronConversionFlag==2){
955 if(theICC){
956 theKnownLevels[i].Icc[j]=theICC->GetICC(theKnownLevels[i].Eg[j],theKnownLevels[i].multipolarity[j]);
957 }
958 }
959 }
960 G4double alpha=theKnownLevels[i].Icc[0];
961 theKnownLevels[i].Pg[0]=theKnownLevels[i].cumulPtot[0]*(1./(alpha+1.));
962 theKnownLevels[i].Pe[0]=theKnownLevels[i].cumulPtot[0]*(alpha/(alpha+1.));
963 for(G4int j=1;j<theKnownLevels[i].NGammas;j++){
964 alpha=theKnownLevels[i].Icc[j];
965 theKnownLevels[i].Pg[j]=(theKnownLevels[i].cumulPtot[j]-theKnownLevels[i].cumulPtot[j-1])*(1./(alpha+1.));
966 theKnownLevels[i].Pe[j]=(theKnownLevels[i].cumulPtot[j]-theKnownLevels[i].cumulPtot[j-1])*(alpha/(alpha+1.));
967 }
968 }
969 }
970
971
972 }
973
974
975
976
977 void G4NuDEXStatisticalNucleus::CreateLevelScheme(){
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 the level scheme
983 NUnknownLevels=0; //will be updated to 1 when creating the capture level
984 }
985 else{
986 //===================================================================
987 //Unknown levels:
988 G4int maxarraysize=EstimateNumberOfLevelsToFill()*1.1/2.+10000;
989 do{
990 maxarraysize*=2;
991 if(theUnkonwnLevels!=0){delete [] theUnkonwnLevels;}
992 //std::cout<<" Max array size of "<<maxarraysize<<std::endl;
993 theUnkonwnLevels=new Level[maxarraysize];
994 NUnknownLevels=GenerateAllUnknownLevels(theUnkonwnLevels,maxarraysize);
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<<" levels in total: "<<NKnownLevels<<" known and "<<NUnknownLevels<<" statistically generated "<<std::endl;
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[NKnownLevels+i]);
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].NLevels;
1017 if(theLevels[i-1].Energy>theLevels[i].Energy){
1018 std::cout<<" ########### Error creating the level scheme ###########"<<std::endl;
1019 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1020 }
1021 }
1022
1023 std::cout<<" NuDEX: Generated statistical nucleus for ZA="<<Z_Int*1000+A_Int<<" up to "<<theLevels[NLevels-1].Energy<<" MeV, with "<<NLevels<<" levels in total: "<<NKnownLevels<<" from the database and "<<NUnknownLevels<<" from statistical models, including bands (without bands --> "<<TotalNIndividualLevels<<" levels). "<<std::endl;
1024
1025 }
1026
1027
1028 void G4NuDEXStatisticalNucleus::CreateThermalCaptureLevel(unsigned int seed){
1029
1030
1031 G4int capturespinx2=((std::fabs(I0)+0.5)*2+0.01); //this spin is always possible ...
1032 G4bool capturepar=true; if(I0<0){capturepar=false;}
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->Integer(4294967295)+1;
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<theThermalCaptureLevel.Energy){
1049 NLevelsBelowThermalCaptureLevel++;
1050 }
1051 }
1052 NLevelsBelowThermalCaptureLevel--; // we exclude the last level, transitions to there have no sense
1053
1054 if(NLevelsBelowThermalCaptureLevel<=0){
1055 NLevelsBelowThermalCaptureLevel=1; //we can go to the ground state (if Sn>0)
1056 }
1057
1058 }
1059
1060 G4int G4NuDEXStatisticalNucleus::GenerateLevelsInSmallRange(G4double Emin,G4double Emax,G4int spinx2,G4bool parity,Level* someLevels,G4int MaxNLevelsToFill){
1061
1062 //If A_Int even/odd --> spinx2 (spin_val*2) should be even/odd
1063 if(((A_Int+spinx2)%2)!=0){
1064 return 0;
1065 }
1066
1067 //Get the level density:
1068 G4double meanNLevels=theLD->Integrate(Emin,Emax,spinx2/2.,parity);
1069
1070 //Sample total number of levels: ????????
1071 G4int thisNLevels=0;
1072 if(meanNLevels>0){
1073 thisNLevels=(G4int)theRandom1->Poisson(meanNLevels);
1074 }
1075
1076 if(thisNLevels>=MaxNLevelsToFill){
1077 std::cout<<" Warning: not enough space to fill levels "<<std::endl;
1078 return -1;
1079 }
1080
1081 //Distribute the levels in the energy interval: ??????
1082 for(G4int i=0;i<thisNLevels;i++){
1083 someLevels[i].Energy=theRandom1->Uniform(Emin,Emax);
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::GenerateLevelsInBigRange(G4double Emin,G4double Emax,G4int spinx2,G4bool parity,Level* someLevels,G4int MaxNLevelsToFill){
1096
1097 G4int TotalNLevels=0;
1098 G4int NIntervals=1000;
1099
1100 // When the LD changes significantly between two levels the Wigner law does not have sense (ot I don't know how to apply it)
1101 // So we will sample the levels according to a Poisson Law when rho(E+<D>)/rho(E) is big and according to Wigner when
1102 // rho(E+<D>)/rho(E) is small, at higher energies. <D> is the mean energy distance between two levels of the same spin and par.
1103 // The difference between "big" and "small" is given by:
1104 G4double WignerRatioThreshold=2;
1105 G4double LevDenThreshold=1; //if there is less than LevDenThreshold/MeV then Wigner has also no sense
1106
1107 for(G4int i=0;i<NIntervals;i++){
1108 G4double emin=Emin+(Emax-Emin)*i/(G4double)NIntervals;
1109 G4double emax=Emin+(Emax-Emin)*(i+1)/(G4double)NIntervals;
1110 G4double meanene=(emin+emax)/2.;
1111 G4double LevDen1=theLD->GetLevelDensity(meanene,spinx2/2.,parity);
1112 if(LevDen1>LevDenThreshold){
1113 G4double LevDen2=theLD->GetLevelDensity(meanene+1./LevDen1,spinx2/2.,parity);
1114 if(LevDen2/LevDen1<WignerRatioThreshold){ //then apply Wigner
1115 //std::cout<<" Wigner way to generate levels abobe "<<emin<<", being E_unk_min = "<<E_unk_min<<std::endl;
1116 G4int newExtraLevels=GenerateWignerLevels(emin,Emax,spinx2,parity,&(someLevels[TotalNLevels]),MaxNLevelsToFill-TotalNLevels);
1117 if(newExtraLevels<0){return -1;}
1118 TotalNLevels+=newExtraLevels;
1119 break;
1120 }
1121 }
1122 //then use Poisson:
1123 G4int newExtraLevels=GenerateLevelsInSmallRange(emin,emax,spinx2,parity,&(someLevels[TotalNLevels]),MaxNLevelsToFill-TotalNLevels);
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*rho*x*x), where x is the energy distance between two levels of the same spin and parity
1133 G4int G4NuDEXStatisticalNucleus::GenerateWignerLevels(G4double Emin,G4double Emax,G4int spinx2,G4bool parity,Level* someLevels,G4int MaxNLevelsToFill){
1134
1135 //If A_Int even/odd --> spinx2 (spin_val*2) should be even/odd
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(previousELevel,spinx2/2.,parity); //levels/MeV
1145 G4double arandgamma=theRandom1->Uniform();
1146 G4double DeltaEMultipliedByLevDen=std::sqrt(-4./3.141592*std::log(1.-arandgamma));
1147 nextELevel=previousELevel+DeltaEMultipliedByLevDen/LevDen;
1148 if(nextELevel<Emax){
1149 someLevels[TotalNLevels].Energy=nextELevel;
1150 someLevels[TotalNLevels].spinx2=spinx2;
1151 someLevels[TotalNLevels].parity=parity;
1152 someLevels[TotalNLevels].seed=0;
1153 someLevels[TotalNLevels].KnownLevelID=-1;
1154 someLevels[TotalNLevels].NLevels=1;
1155 someLevels[TotalNLevels].Width=0;
1156 TotalNLevels++;
1157 if(TotalNLevels>=MaxNLevelsToFill){
1158 std::cout<<" Warning: not enough space to fill levels "<<std::endl;
1159 //NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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, not individually:
1172 G4int G4NuDEXStatisticalNucleus::GenerateBandLevels(G4int bandmin,G4int bandmax,G4int spinx2,G4bool parity,Level* someLevels,G4int MaxNLevelsToFill){
1173
1174 //If A_Int even/odd --> spinx2 (spin_val*2) should be even/odd
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(__LINE__).c_str(),"##### Error in NuDEX #####");
1185 }
1186
1187 for(G4int i=bandmin;i<=bandmax;i++){
1188 G4double emin=Emin+(Emax-Emin)*i/(G4double)NBands;
1189 G4double emax=Emin+(Emax-Emin)*(i+1.)/(G4double)NBands;
1190 G4double AverageNumberOfLevels=theLD->Integrate(emin,emax,spinx2/2.,parity);
1191 G4int NumberOfLevelsInThisBand=0;
1192 if(AverageNumberOfLevels>0){
1193 NumberOfLevelsInThisBand=(G4int)theRandom1->Poisson(AverageNumberOfLevels);
1194 }
1195 if(NumberOfLevelsInThisBand>0){
1196 someLevels[TotalNLevels].Energy=(emax+emin)/2.;
1197 someLevels[TotalNLevels].spinx2=spinx2;
1198 someLevels[TotalNLevels].parity=parity;
1199 someLevels[TotalNLevels].seed=0;
1200 someLevels[TotalNLevels].KnownLevelID=-1;
1201 someLevels[TotalNLevels].NLevels=NumberOfLevelsInThisBand;
1202 someLevels[TotalNLevels].Width=emax-emin;
1203 TotalNLevels++;
1204 if(TotalNLevels>=MaxNLevelsToFill){
1205 std::cout<<" Warning: not enough space to fill levels "<<std::endl;
1206 //NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1207 return -1;
1208 }
1209 }
1210 }
1211
1212 return TotalNLevels;
1213 }
1214
1215
1216
1217 G4int G4NuDEXStatisticalNucleus::GenerateAllUnknownLevels(Level* someLevels,G4int MaxNLevelsToFill){
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_val*2) should be even/odd
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 and E_unk_max
1231 //We will create the levels one by one at low energies and directly in bands at higher energies
1232 //The limit between the two ranges will be given by E_lim_onebyone
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; // band corresponding to E_lim_onebyone
1237
1238 //----------------------------------------------------
1239 //Calculate E_lim_onebyone:
1240 #ifndef GENERATEEXPLICITLYALLLEVELSCHEME
1241 if(NBands>0){
1242 if(MinLevelsPerBand<=0){ // All the level scheme in bands
1243 E_lim_onebyone=0;
1244 i_Band_E_lim_onebyone=0;
1245 }
1246 else{
1247 G4double bandwidth=(Emax_bands-Emin_bands)/NBands;
1248 G4double rho_lim_onebyone=3.*(MinLevelsPerBand+10.)/bandwidth; // above this energy we start the creation of bands without sampling the levels one by one
1249 E_lim_onebyone=theLD->EstimateInverse(rho_lim_onebyone,spinx2/2.,parity);
1250 }
1251 }
1252 if(E_unk_max-Emax_bands>0.001){ //then E_unk_max>Emax_bands and we generate all the levels explicitly
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 between two bands:
1258 if(E_lim_onebyone>E_unk_min && E_lim_onebyone<E_unk_max){
1259 for(G4int i=0;i<NBands;i++){
1260 G4double elow_band=Emin_bands+(Emax_bands-Emin_bands)*i/(G4double)NBands;
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 have to create some of the levels one by one
1273 if(E_lim_onebyone<Emax){
1274 Emax=E_lim_onebyone;
1275 }
1276 NLev=GenerateLevelsInBigRange(Emin,Emax,spinx2,parity,&(someLevels[TotalNLevels]),MaxNLevelsToFill-TotalNLevels);
1277 if(NLev<0){return -1;}
1278 if(NBands>0 && NLev>0){
1279 NLev=CreateBandsFromLevels(NLev,&(someLevels[TotalNLevels]),spinx2,parity);
1280 }
1281 TotalNLevels+=NLev;
1282 }
1283
1284 if(i_Band_E_lim_onebyone<NBands){ //then we have to create some of the levels directly with bands
1285 NLev=GenerateBandLevels(i_Band_E_lim_onebyone,NBands-1,spinx2,parity,&(someLevels[TotalNLevels]),MaxNLevelsToFill-TotalNLevels);
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), ComparisonLevels);
1297
1298 return TotalNLevels;
1299 }
1300
1301 //Junta varios niveles en uno solo, creando distintas bandas, para un spin y paridad determinados.
1302 //Entiende que no hay otros spines ni paridades. Si hay otros, peta.
1303 //Devuelve el numero de niveles actualizado.
1304 G4int G4NuDEXStatisticalNucleus::CreateBandsFromLevels(G4int thisNLevels,Level* someLevels,G4int spinx2,G4bool parity){
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/(G4double)NBands;
1312 G4double emax=Emin+(Emax-Emin)*(i+1.)/(G4double)NBands;
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 || someLevels[j].parity!=parity){
1323 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1324 }
1325 if(someLevels[j].Energy>=emin && someLevels[j].Energy<=emax){
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[j].Energy<=emax){
1332 theBandLevels[i].NLevels+=someLevels[j].NLevels;
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],&theBandLevels[i]);
1346 }
1347 i--;
1348 FinalNBands--;
1349 }
1350 }
1351
1352 //Replace levels with bands and update number of levels:
1353 G4int NbandsCopied=0;
1354 for(G4int i=0;i<thisNLevels;i++){
1355 if(someLevels[i].Energy<0){
1356 if(NbandsCopied<FinalNBands){ //this level is replaced by a band
1357 CopyLevel(&theBandLevels[NbandsCopied],&someLevels[i]);
1358 NbandsCopied++;
1359 }
1360 else{ //there is no band to replace. Copy the last level here
1361 CopyLevel(&someLevels[thisNLevels-1],&someLevels[i]);
1362 i--;
1363 thisNLevels--;
1364 }
1365 }
1366 }
1367
1368 //Simple check:
1369 if(NbandsCopied!=FinalNBands){
1370 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1371 }
1372 delete [] theBandLevels;
1373 return thisNLevels;
1374 }
1375
1376
1377 //Esto lo que hace es intentar reemplazar niveles de la parte estadistica por los que se conocen:
1378 G4int G4NuDEXStatisticalNucleus::InsertHighEnergyKnownLevels(){
1379
1380
1381
1382 G4bool* HasBeenInserted=new G4bool[KnownLevelsVectorSize];
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: first levels with NGammas>0, then the rest ...
1388 for(G4int k=1;k<5;k++){ //loop in the distance between levels condition
1389 G4double MaxEnergyDistance=0.1*k; //The level to replace should be close to it, so first we try with X MeV and then 2X MeV ...
1390 for(G4int i=NKnownLevels;i<KnownLevelsVectorSize;i++){ //loop in the known levels
1391 if(theKnownLevels[i].Energy>Sn){break;}
1392 if(HasBeenInserted[i]==false && (theKnownLevels[i].NGammas>0 || kk>0)){
1393 //--------------------------------------------------------------------------------
1394 G4double MinEnergyDistance=-1,EnergyDistance=0;
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++){ // loop in the unknown levels
1399 if(theLevels[j].spinx2==thespinx2 && theLevels[j].parity==thepar){
1400 EnergyDistance=std::fabs(theKnownLevels[i].Energy-theLevels[j].Energy);
1401 if((EnergyDistance<MinEnergyDistance || MinEnergyDistance<0) && EnergyDistance<MaxEnergyDistance && theLevels[j].KnownLevelID<0){
1402 MinEnergyDistance=EnergyDistance;
1403 unknownLevelID=j;
1404 }
1405 //else if(EnergyDistance>MinEnergyDistance && MinEnergyDistance>0){
1406 //break;
1407 //}
1408 }
1409 }
1410 if(unknownLevelID>0 && theLevels[unknownLevelID].NLevels==1){ //then we replace the stat-level by the known level:
1411 //std::cout<<" Copy Level "<<i<<" with E= "<<theKnownLevels[i].Energy<<" into level "<<unknownLevelID<<" with E="<<theLevels[unknownLevelID].Energy<<std::endl;
1412 CopyLevel(&theKnownLevels[i],&theLevels[unknownLevelID]);
1413 theLevels[unknownLevelID].KnownLevelID=i;
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), ComparisonLevels);
1425
1426
1427
1428 //-----------------------------------------------------------------------------
1429 //we cannot go to from an inserted level with NGammas>0 to the statistical part of the level scheme (the level ID will then be different in the known and unknown level vectors. There is one exception: if the level has been inserted.
1430 //so, if it is the case, we change the final level for the closest one:
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].NGammas;j++){
1435 if(theKnownLevels[knownID].FinalLevelID[j]>=NKnownLevels){//this cannot be
1436 //-----------------------------------------------------
1437 G4int i_finalknownlevel=theKnownLevels[knownID].FinalLevelID[j];
1438 G4double MinEnergyDistance=-1;
1439 G4int i_statlevel=-1;
1440 for(G4int k=0;k<i;k++){
1441 G4double EnergyDistance=std::fabs(theKnownLevels[i_finalknownlevel].Energy-theLevels[k].Energy);
1442 if(EnergyDistance<MinEnergyDistance || MinEnergyDistance<0){
1443 MinEnergyDistance=EnergyDistance;
1444 i_statlevel=k;
1445 }
1446 }
1447 if(MinEnergyDistance<0){
1448 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1449 }
1450 //std::cout<<" Final known level "<<i_finalknownlevel<<" with E = "<<theKnownLevels[i_finalknownlevel].Energy<<" has been replaced by final level "<<i_statlevel<<" with E = "<<theLevels[i_statlevel].Energy<<std::endl;
1451 if(theLevels[i_statlevel].KnownLevelID>=0){ //then is an inserted level, we change only the final level
1452 theKnownLevels[knownID].FinalLevelID[j]=i_statlevel;
1453 }
1454 else{ //is a real statistical level. We change the multipolarity and the alpha:
1455 theKnownLevels[knownID].FinalLevelID[j]=i_statlevel;
1456 theKnownLevels[knownID].multipolarity[j]=GetMultipolarity(&theLevels[i],&theLevels[i_statlevel]);
1457 theKnownLevels[knownID].Eg[j]=theLevels[i].Energy-theLevels[i_statlevel].Energy;
1458 theKnownLevels[knownID].Pg[j]=theKnownLevels[knownID].Pg[j]+theKnownLevels[knownID].Pe[j];
1459 theKnownLevels[knownID].Pe[j]=0;
1460 theKnownLevels[knownID].Icc[j]=0; //set to 0 to avoid problems with the electron conversion
1461 }
1462 //-----------------------------------------------------
1463 }
1464 }
1465 }
1466 }
1467 //-----------------------------------------------------------------------------
1468
1469 return 0;
1470 }
1471
1472
1473
1474 G4int G4NuDEXStatisticalNucleus::EstimateNumberOfLevelsToFill(){
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,meanNLevels;
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)NIntervals;
1495 emax=Emin+(Emax-Emin)*(i+1)/(G4double)NIntervals;
1496 meanEnergy=(emax+emin)/2.;
1497
1498 //Positive parity:
1499 LevDen=theLD->GetLevelDensity(meanEnergy,spinx2/2.,true); //levels/MeV
1500 meanNLevels=LevDen*(emax-emin);
1501 TotalNLevels+=meanNLevels;
1502 TotalNLevesWithSameSpin+=meanNLevels;
1503 if(NBands>0 && meanEnergy>Emin_bands && meanEnergy<Emax_bands){ //if there are bands and these levels go inside:
1504 TotalNLevelsInsideBands+=meanNLevels;
1505 }
1506
1507
1508 //Negative parity:
1509 LevDen=theLD->GetLevelDensity(meanEnergy,spinx2/2.,false); //levels/MeV
1510 meanNLevels=LevDen*(emax-emin);
1511 TotalNLevels+=meanNLevels;
1512 TotalNLevesWithSameSpin+=meanNLevels;
1513 if(NBands>0 && meanEnergy>Emin_bands && meanEnergy<Emax_bands){ //if there are bands and these levels go inside:
1514 TotalNLevelsInsideBands+=meanNLevels;
1515 }
1516
1517 }
1518 if(MaxNLevelsWithSameSpinAndParity<TotalNLevesWithSameSpin){MaxNLevelsWithSameSpinAndParity=TotalNLevesWithSameSpin;}
1519 }
1520
1521 MaxNLevelsWithSameSpinAndParity/=2.;
1522
1523 if(TotalNLevelsInsideBands>0){
1524 //then there are some levels inside bands and the amount of levels needed will be reduced:
1525 G4double TotalNLevelsOutsideBands=TotalNLevels-TotalNLevelsInsideBands;
1526 G4double MaxNLevelsNeededForBands=NBands*maxspinx2*MinLevelsPerBand;
1527 TotalNLevels=MaxNLevelsWithSameSpinAndParity+TotalNLevelsOutsideBands+MaxNLevelsNeededForBands;
1528 }
1529
1530 return (G4int)TotalNLevels;
1531 }
1532
1533
1534 G4double G4NuDEXStatisticalNucleus::TakeTargetNucleiI0(const char* fname,G4int& check){
1535
1536 std::ifstream in(fname);
1537 if(!in.good()){
1538 std::cout<<" ######## Error opening file "<<fname<<" ########"<<std::endl;
1539 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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(__LINE__).c_str(),"##### Error in NuDEX #####");
1557 }
1558 in.ignore(10000,'\n');
1559 G4double spin,par;
1560 in.get(buffer,16);
1561 in.get(buffer,6); spin=std::fabs(atof(buffer)); // some spins are negative ???
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::ReadKnownLevels(const char* fname){
1572
1573
1574 std::ifstream in(fname);
1575 if(!in.good()){
1576 std::cout<<" ######## Error opening file "<<fname<<" ########"<<std::endl;
1577 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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=atoi(buffer);
1587 in.get(buffer,16);
1588 in.get(buffer,13); G4double newSn=atof(buffer);
1589 if(Sn>0 && std::fabs(Sn-newSn)>0.01){
1590 std::cout<<" ######## WARNING: Sn value from the level density file ("<<Sn<<") is different than the one from the known levels file ("<<newSn<<"). We use the first value. ########"<<std::endl;
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[KnownLevelsVectorSize];
1608 for(G4int i=0;i<KnownLevelsVectorSize;i++){theKnownLevels[i].NGammas=0;}
1609 G4double spin,par;
1610 for(G4int i=0;i<KnownLevelsVectorSize;i++){
1611 in.get(buffer,4); theKnownLevels[i].id=atoi(buffer)-1;
1612 in.get(buffer,2);
1613 in.get(buffer,11); theKnownLevels[i].Energy=atof(buffer);
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].Energy<Ecrit){
1618 std::cout<<" ######## WARNING: Spin and parity for level "<<i<<" is s="<<spin<<" p="<<par<<" for Z="<<Z_Int<<", A="<<A_Int<<" ########"<<std::endl;
1619 if(spin<0){
1620 spin=0;
1621 if(i>1){ //Random spin, same as one of the levels below this one:
1622 G4int i_sampleSpin=theRandom1->Integer(i-1);
1623 spin=theKnownLevels[i_sampleSpin].spinx2/2.;
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=atof(buffer);
1633 in.get(buffer,4); theKnownLevels[i].NGammas=atoi(buffer);
1634 if(theKnownLevels[i].NGammas>0){
1635 if(spin<0){
1636 spin=0;
1637 if(i>1){ //Random spin, same as one of the levels below this one:
1638 G4int i_sampleSpin=theRandom1->Integer(i-1);
1639 spin=theKnownLevels[i_sampleSpin].spinx2/2.;
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.01);
1648 if(par>0){theKnownLevels[i].parity=true;}else{theKnownLevels[i].parity=false;}
1649
1650 //---------------------------------
1651 //decay modes:
1652 in.get(buffer,27); //dummy
1653 in.get(buffer,4);
1654 G4int decays=theKnownLevels[i].Ndecays=atoi(buffer);
1655 if(decays>0){
1656 theKnownLevels[i].decayFraction=new G4double[decays];
1657 theKnownLevels[i].decayMode=new std::string[decays];
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(buffer);
1663 in.get(buffer,2);
1664 in.get(buffer,8);
1665 theKnownLevels[i].decayMode[j]=std::string(buffer);
1666 }
1667 //----------------------------------
1668
1669 in.ignore(10000,'\n');
1670
1671 if(theKnownLevels[i].NGammas>0){
1672 theKnownLevels[i].FinalLevelID=new G4int[theKnownLevels[i].NGammas];
1673 theKnownLevels[i].multipolarity=new G4int[theKnownLevels[i].NGammas];
1674 theKnownLevels[i].Eg=new G4double[theKnownLevels[i].NGammas];
1675 theKnownLevels[i].cumulPtot=new G4double[theKnownLevels[i].NGammas];
1676 theKnownLevels[i].Pg=new G4double[theKnownLevels[i].NGammas];
1677 theKnownLevels[i].Pe=new G4double[theKnownLevels[i].NGammas];
1678 theKnownLevels[i].Icc=new G4double[theKnownLevels[i].NGammas];
1679 }
1680
1681 for(G4int j=0;j<theKnownLevels[i].NGammas;j++){
1682 in.get(buffer,40);
1683
1684 in.get(buffer,5); theKnownLevels[i].FinalLevelID[j]=atoi(buffer)-1;
1685 theKnownLevels[i].multipolarity[j]=0;
1686 in.get(buffer,2);
1687 in.get(buffer,11); theKnownLevels[i].Eg[j]=atof(buffer);
1688 in.get(buffer,2);
1689 in.get(buffer,11); theKnownLevels[i].Pg[j]=atof(buffer);
1690 in.get(buffer,2);
1691 in.get(buffer,11); theKnownLevels[i].Pe[j]=atof(buffer);
1692 in.get(buffer,2);
1693 in.get(buffer,11); theKnownLevels[i].Icc[j]=atof(buffer);
1694 theKnownLevels[i].cumulPtot[j]=theKnownLevels[i].Pg[j]*(1+theKnownLevels[i].Icc[j]); //we rely in Pg and Icc, where Icc=Pe/Pg
1695 in.ignore(10000,'\n');
1696 if(theKnownLevels[i].FinalLevelID[j]>=i+1){
1697 std::cout<<" ######## Error reading file "<<fname<<" ########"<<std::endl;
1698 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1699 }
1700 }
1701
1702 if(theKnownLevels[i].id!=i || !in.good()){
1703 std::cout<<" ######## Error reading file "<<fname<<" ########"<<std::endl;
1704 std::cout<<" Level "<<i<<" has id = "<<theKnownLevels[i].id<<std::endl;
1705 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1706 }
1707
1708 //--------------------------------------------------------------------------------------
1709 //normalize, and put cumulPtot as cumulative. Re-calculate Pe
1710 G4double totalEmissionProb=0;
1711 for(G4int j=0;j<theKnownLevels[i].NGammas;j++){
1712 totalEmissionProb+=theKnownLevels[i].cumulPtot[j];
1713 }
1714 //------------------------------------
1715 if(totalEmissionProb==0){//sometimes all the levels have 0 prob.
1716 for(G4int j=0;j<theKnownLevels[i].NGammas;j++){
1717 theKnownLevels[i].cumulPtot[j]=1./theKnownLevels[i].NGammas;
1718 theKnownLevels[i].Pg[j]=theKnownLevels[i].cumulPtot[j]/(1.+theKnownLevels[i].Icc[j]);
1719 }
1720 totalEmissionProb=1;
1721 }
1722 //------------------------------------
1723 for(G4int j=0;j<theKnownLevels[i].NGammas;j++){
1724 theKnownLevels[i].cumulPtot[j]/=totalEmissionProb;
1725 theKnownLevels[i].Pg[j]/=totalEmissionProb;
1726 theKnownLevels[i].Pe[j]=theKnownLevels[i].Icc[j]*theKnownLevels[i].Pg[j];
1727 }
1728 for(G4int j=1;j<theKnownLevels[i].NGammas;j++){
1729 theKnownLevels[i].cumulPtot[j]+=theKnownLevels[i].cumulPtot[j-1];
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(const char* fname){
1755
1756 std::ifstream in(fname);
1757 if(!in.good()){
1758 std::cout<<" ######## Error opening file "<<fname<<" ########"<<std::endl;
1759 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
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>>word>>word>>word>>Ecrit; break;
1769 }
1770 in.ignore(10000,'\n');
1771 }
1772 in.close();
1773
1774 return Ecrit;
1775 }
1776
1777 G4int G4NuDEXStatisticalNucleus::ReadSpecialInputFile(const char* fname){
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,'\n');}
1788 if(word==std::string("END")){break;}
1789 //now we will take values only if they have not been set yet:
1790 else if(word==std::string("LEVELDENSITYTYPE")){if(LevelDensityType<0){in>>LevelDensityType;}}
1791 else if(word==std::string("MAXSPIN")){if(maxspinx2<0){in>>MaxSpin; maxspinx2=(G4int)(2.*MaxSpin+0.01);}}
1792 else if(word==std::string("MINLEVELSPERBAND")){if(MinLevelsPerBand<0){in>>MinLevelsPerBand;}}
1793 else if(word==std::string("BANDWIDTH_MEV")){if(BandWidth==0){in>>BandWidth;}}
1794 else if(word==std::string("MAXEXCENERGY_MEV")){if(MaxExcEnergy==0){in>>MaxExcEnergy;}}
1795 else if(word==std::string("ECRIT_MEV")){if(Ecrit<0){in>>Ecrit;}}
1796 else if(word==std::string("KNOWNLEVELSFLAG")){if(KnownLevelsFlag<0){in>>KnownLevelsFlag;}}
1797
1798 else if(word==std::string("PSF_FLAG")){if(PSFflag<0){in>>PSFflag;}}
1799 else if(word==std::string("BROPTION")){if(BROpt<0){in>>BROpt;}}
1800 else if(word==std::string("SAMPLEGAMMAWIDTHS")){if(SampleGammaWidths<0){in>>SampleGammaWidths;}}
1801
1802 else if(word==std::string("ELECTRONCONVERSIONFLAG")){if(ElectronConversionFlag<0){in>>ElectronConversionFlag;}}
1803 else if(word==std::string("PRIMARYTHCAPGAMNORM")){if(PrimaryGammasIntensityNormFactor<0){in>>PrimaryGammasIntensityNormFactor;}}
1804 else if(word==std::string("PRIMARYGAMMASECUT")){if(PrimaryGammasEcut<0){in>>PrimaryGammasEcut;}}
1805 }
1806 in.close();
1807 return 1;
1808 }
1809
1810 G4int G4NuDEXStatisticalNucleus::ReadGeneralStatNuclParameters(const char* fname){
1811
1812 std::ifstream in(fname);
1813 if(!in.good()){
1814 std::cout<<" ######## Error opening file "<<fname<<" ########"<<std::endl;
1815 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1816 }
1817 char line[1000];
1818 in.getline(line,1000);
1819 in.getline(line,1000);
1820 G4int tmpZ,tmpA,tmpLDtype,tmpPSFflag,tmpMaxSpin,tmpMinLev,tmpBrOption,tmpSampleGW;
1821 G4double tmpBandWidth,tmpMaxExcEnergy;
1822 unsigned int tmpseed1,tmpseed2,tmpseed3;
1823 G4int finalLDtype=0,finalPSFflag=0,finalMaxSpin=0,finalMinLev=0,finalBrOption=0,finalSampleGW=0;
1824 G4double finalBandWidth=0,finalMaxExcEnergy=0;
1825 unsigned int finalseed1=0,finalseed2=0,finalseed3=0;
1826 G4bool NuclDataHasBeenRead=false;
1827 G4bool DefaulDataHasBeenRead=false;
1828 while(in>>tmpZ>>tmpA>>tmpLDtype>>tmpPSFflag>>tmpMaxSpin>>tmpMinLev>>tmpBandWidth>>tmpMaxExcEnergy>>tmpBrOption>>tmpSampleGW>>tmpseed1>>tmpseed2>>tmpseed3){
1829 G4bool TakeData=false;
1830 if(tmpZ==Z_Int && tmpA==A_Int){ //then this is our nucleus
1831 NuclDataHasBeenRead=true;
1832 TakeData=true;
1833 }
1834 else if(tmpZ==0 && tmpA==0 && NuclDataHasBeenRead==false){ //default, only if our nucleus has not been read
1835 DefaulDataHasBeenRead=true;
1836 TakeData=true;
1837 }
1838 if(TakeData){
1839 finalLDtype=tmpLDtype; finalPSFflag=tmpPSFflag; finalMaxSpin=tmpMaxSpin; finalMinLev=tmpMinLev; finalBrOption=tmpBrOption; finalSampleGW=tmpSampleGW;
1840 finalBandWidth=tmpBandWidth; finalMaxExcEnergy=tmpMaxExcEnergy;
1841 finalseed1=tmpseed1; finalseed2=tmpseed2; finalseed3=tmpseed3;
1842 }
1843 }
1844 in.close();
1845
1846 if(NuclDataHasBeenRead==false && DefaulDataHasBeenRead==false){
1847 std::cout<<" ######## Error reading "<<fname<<" ########"<<std::endl;
1848 NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1849 }
1850
1851 // Replace data only if it has not been set by the SetSomeInitalParameters method:
1852 if(LevelDensityType<0){LevelDensityType=finalLDtype;}
1853 if(PSFflag<0){PSFflag=finalPSFflag;}
1854 if(maxspinx2<0){maxspinx2=(G4int)(2.*finalMaxSpin+0.01);}
1855 if(MinLevelsPerBand<0){MinLevelsPerBand=finalMinLev;}
1856 if(BandWidth==0){BandWidth=finalBandWidth;}
1857 if(MaxExcEnergy==0){MaxExcEnergy=finalMaxExcEnergy;}
1858 if(BROpt<0){BROpt=finalBrOption;}
1859 if(SampleGammaWidths<0){SampleGammaWidths=finalSampleGW;}
1860 if(Rand1seedProvided==false){seed1=finalseed1; theRandom1->SetSeed(finalseed1);}
1861 if(Rand2seedProvided==false){seed2=finalseed2; theRandom2->SetSeed(finalseed2);}
1862 if(Rand3seedProvided==false){seed3=finalseed3; theRandom3->SetSeed(finalseed3);}
1863
1864 //-------------------------------------------------------
1865 //Now some checks:
1866 if(maxspinx2<1){
1867 std::cout<<" ######## Error: maximum spin for generating the statistical nucleus with A="<<A_Int<<" and Z="<<Z_Int<<" has been set to "<<maxspinx2/2.<<" ########"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1868 }
1869 if(MinLevelsPerBand<=0 && BandWidth<0){
1870 std::cout<<" ######## Error: MinLevelsPerBand and BandWidth for generating the statistical nucleus with A="<<A_Int<<" and Z="<<Z_Int<<" has been set to MinLevelsPerBand="<<MinLevelsPerBand<<" and BandWidth="<<BandWidth<<" ########"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1871 }
1872 if(BROpt!=0 && BROpt!=1 && BROpt!=2){
1873 std::cout<<" ######## Error: BROpt for generating the statistical nucleus with A="<<A_Int<<" and Z="<<Z_Int<<" has been set to BROpt="<<BROpt<<", and has to be BROpt=0,1 or 2 ########"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1874 }
1875 if(SampleGammaWidths!=0 && SampleGammaWidths!=1){
1876 std::cout<<" ######## Error: SampleGammaWidths parameter for generating the statistical nucleus with A="<<A_Int<<" and Z="<<Z_Int<<" has been set to SampleGammaWidths="<<SampleGammaWidths<<", and has to be SampleGammaWidths=0 or 1 ########"<<std::endl; NuDEXException(__FILE__,std::to_string(__LINE__).c_str(),"##### Error in NuDEX #####");
1877 }
1878 //-------------------------------------------------------
1879
1880
1881 return 0;
1882 }
1883
1884
1885 void G4NuDEXStatisticalNucleus::GenerateThermalCaptureLevelBR(const char* dirname){
1886
1887 char fname[1000];
1888 snprintf(fname,1000,"%s/PrimaryCaptureGammas.dat",dirname);
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 file, if they are there:
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 = "<<ELevel<<" and Sn = "<<Sn<<" for ZA = "<<aZA<<" ##########"<<std::endl;
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; //renormalization
1915 //TotalThI+=ThI[i];
1916 }
1917 break;
1918 }
1919 }
1920 in.ignore(10000,'\n');
1921 }
1922 in.close();
1923
1924 if(theThermalCaptureLevelCumulBR){delete [] theThermalCaptureLevelCumulBR;}
1925 theThermalCaptureLevelCumulBR=new G4double[NLevelsBelowThermalCaptureLevel]; //the final result
1926 for(G4int i=0;i<NLevelsBelowThermalCaptureLevel;i++){
1927 theThermalCaptureLevelCumulBR[i]=0;
1928 }
1929
1930 //------------------------------------------------------------------------------------------------------
1931 //We calculate which transitrions go to an existing level and the total intensity of the transitions:
1932 G4double totalThGInt=0,ENDSFLevelEnergy=0,MinlevelDist=0,MinlevelDist_known=0,LevDist=0;
1933 G4int i_closest=0,i_closest_known=0;
1934 G4double MaxAllowedLevelDistance=0.010; //10 keV
1935 G4bool ComputePrimaryGammasEcut=false;
1936 if(PrimaryGammasEcut==0){ComputePrimaryGammasEcut=true;}
1937
1938 //We take only those intensities going to our levels:
1939 for(G4int i=0;i<ng;i++){
1940 ENDSFLevelEnergy=ELevel-ThEg[i];
1941 if(ComputePrimaryGammasEcut && PrimaryGammasEcut<ENDSFLevelEnergy){
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<NLevelsBelowThermalCaptureLevel;j++){
1949 LevDist=std::fabs(ENDSFLevelEnergy-theLevels[j].Energy);
1950 if(theLevels[j].KnownLevelID>=0){ //then this is a known level. We priorize known levels.
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<MaxAllowedLevelDistance){ // We priorize known levels.
1962 theThermalCaptureLevelCumulBR[i_closest_known]=ThI[i];
1963 totalThGInt+=theThermalCaptureLevelCumulBR[i_closest_known];
1964 }
1965 else if(MinlevelDist<MaxAllowedLevelDistance){
1966 theThermalCaptureLevelCumulBR[i_closest]=ThI[i];
1967 totalThGInt+=theThermalCaptureLevelCumulBR[i_closest];
1968 }
1969 }
1970 //if(TotalThI>0){std::cout<<" NuDEX: Primary thermal gammas for ZA="<<Z_Int*1000+A_Int<<" file: "<<TotalThI<<" accepted: "<<totalThGInt<<" ratio: "<<totalThGInt/TotalThI*100.<<" %"<<std::endl;}
1971 //else{std::cout<<" Primary thermal gammas for "<<Z_Int*1000+A_Int<<" file: "<<TotalThI<<std::endl;}
1972 std::cout<<" NuDEX: Primary thermal gammas for ZA="<<Z_Int*1000+A_Int<<" found in the database: "<<totalThGInt*100.<<" %"<<std::endl;
1973 //------------------------------------------------------------------------------------------------------
1974
1975 //------------------------------------------------------------------------------------------------------
1976 //If we don't have all the info, we take the rest of the BR from one of the levels, but with a proper normalization:
1977 if(totalThGInt<0.95){
1978 G4double TotalNeededIntensity=1.-totalThGInt;
1979 G4double oldInt,TotalOldIntensity=0;
1980 //-------------------------
1981 //We put the capture level replacing one of the levels and calculate the BR:
1982 Level tmpLevel;
1983 CopyLevel(&theLevels[NLevelsBelowThermalCaptureLevel],&tmpLevel);
1984 CopyLevel(&theThermalCaptureLevel,&theLevels[NLevelsBelowThermalCaptureLevel]);
1985 G4double* CumulBR_th_v2=new G4double[NLevelsBelowThermalCaptureLevel];
1986 ComputeDecayIntensities(NLevelsBelowThermalCaptureLevel,CumulBR_th_v2);
1987 CopyLevel(&tmpLevel,&theLevels[NLevelsBelowThermalCaptureLevel]);
1988 //-------------------------
1989
1990 for(G4int i=0;i<NLevelsBelowThermalCaptureLevel;i++){
1991 if(i==0){oldInt=CumulBR_th_v2[i];}else{oldInt=CumulBR_th_v2[i]-CumulBR_th_v2[i-1];}
1992 if(theThermalCaptureLevelCumulBR[i]==0 && theLevels[i].Energy>=PrimaryGammasEcut){TotalOldIntensity+=oldInt;}
1993 }
1994
1995 if(TotalOldIntensity>0){
1996 for(G4int i=0;i<NLevelsBelowThermalCaptureLevel;i++){
1997 if(theThermalCaptureLevelCumulBR[i]==0 && theLevels[i].Energy>=PrimaryGammasEcut){
1998 if(i==0){oldInt=CumulBR_th_v2[i];}else{oldInt=CumulBR_th_v2[i]-CumulBR_th_v2[i-1];}
1999 theThermalCaptureLevelCumulBR[i]=oldInt*TotalNeededIntensity/TotalOldIntensity;
2000 }
2001 }
2002 }
2003 delete [] CumulBR_th_v2;
2004 }
2005 //------------------------------------------------------------------------------------------------------
2006
2007 //theThermalCaptureLevelCumulBR still not cumulative
2008
2009 for(G4int i=1;i<NLevelsBelowThermalCaptureLevel;i++){
2010 theThermalCaptureLevelCumulBR[i]+=theThermalCaptureLevelCumulBR[i-1];
2011 }
2012 for(G4int i=0;i<NLevelsBelowThermalCaptureLevel;i++){
2013 theThermalCaptureLevelCumulBR[i]/=theThermalCaptureLevelCumulBR[NLevelsBelowThermalCaptureLevel-1];
2014 }
2015
2016 }
2017
2018 //Cambiamos las intensidades de los "primary gammas". Del correspondiente al ninvel con energÃa "LevelEnergy"
2019 void G4NuDEXStatisticalNucleus::ChangeThermalCaptureLevelBR(G4double LevelEnergy,G4double absoluteIntensity){
2020
2021 if(!theThermalCaptureLevelCumulBR){return;}
2022 G4int level_id=GetClosestLevel(LevelEnergy,-1,true);
2023 if(level_id<0 || level_id>=NLevelsBelowThermalCaptureLevel){
2024 std::cout<<" ############## WARNING in "<<__FILE__<<", line "<<__LINE__<<" ##############"<<std::endl;
2025 std::cout<<" ---> "<<level_id<<" "<<LevelEnergy<<std::endl;
2026 }
2027
2028 for(G4int i=NLevelsBelowThermalCaptureLevel-1;i>0;i--){
2029 theThermalCaptureLevelCumulBR[i]-=theThermalCaptureLevelCumulBR[i-1];
2030 }
2031 G4double OldIntensity=theThermalCaptureLevelCumulBR[level_id];
2032 theThermalCaptureLevelCumulBR[level_id]=absoluteIntensity*(1.-OldIntensity)/(1.-absoluteIntensity);
2033
2034 for(G4int i=1;i<NLevelsBelowThermalCaptureLevel;i++){
2035 theThermalCaptureLevelCumulBR[i]+=theThermalCaptureLevelCumulBR[i-1];
2036 }
2037 for(G4int i=0;i<NLevelsBelowThermalCaptureLevel;i++){
2038 theThermalCaptureLevelCumulBR[i]/=theThermalCaptureLevelCumulBR[NLevelsBelowThermalCaptureLevel-1];
2039 }
2040 if(level_id==0){
2041 std::cout<<" Thermal primary gammas to level "<<level_id<<", with E="<<theLevels[level_id].Energy<<" MeV changed from "<<OldIntensity<<" to "<<theThermalCaptureLevelCumulBR[level_id]<<std::endl;
2042 }
2043 else{
2044 std::cout<<" Thermal primary gammas to level "<<level_id<<", with E="<<theLevels[level_id].Energy<<" MeV changed from "<<OldIntensity<<" to "<<theThermalCaptureLevelCumulBR[level_id]-theThermalCaptureLevelCumulBR[level_id-1]<<std::endl;
2045 }
2046 }
2047
2048 void G4NuDEXStatisticalNucleus::PrintParameters(std::ostream &out){
2049
2050 out<<" ###################################################################################### "<<std::endl;
2051 out<<" GENERAL_PARS"<<std::endl;
2052 out<<" Z = "<<Z_Int<<" A = "<<A_Int<<std::endl;
2053 out<<" Sn = "<<Sn<<" I0(ZA-1) = "<<I0<<std::endl;
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_unknown_max = "<<E_unk_max<<std::endl;
2059 out<<" maxspin = "<<maxspinx2/2.<<std::endl;
2060 out<<" MaxExcEnergy = "<<MaxExcEnergy<<std::endl;
2061 out<<" NBands = "<<NBands<<" MinLevelsPerBand = "<<MinLevelsPerBand<<" BandWidth = "<<BandWidth<<std::endl;
2062 out<<" Emin_bands = "<<Emin_bands<<" Emax_bands = "<<Emax_bands<<std::endl;
2063 out<<" NLevels = "<<NLevels<<" NKnownLevels = "<<NKnownLevels<<" NUnknownLevels = "<<NUnknownLevels<<std::endl;
2064 out<<" BROpt = "<<BROpt<<" SampleGammaWidths = "<<SampleGammaWidths<<std::endl;
2065 out<<" PrimaryGammasIntensityNormFactor = "<<PrimaryGammasIntensityNormFactor<<" PrimaryGammasEcut = "<<PrimaryGammasEcut<<std::endl;
2066 out<<" KnownLevelsFlag = "<<KnownLevelsFlag<<std::endl;
2067 out<<" ElectronConversionFlag = "<<ElectronConversionFlag<<std::endl;
2068 out<<" ###################################################################################### "<<std::endl;
2069
2070 }
2071
2072 void G4NuDEXStatisticalNucleus::PrintKnownLevels(std::ostream &out){
2073
2074 out<<" ########################################################################################################## "<<std::endl;
2075 out<<" KNOWN_LEVEL_SCHEME "<<std::endl;
2076 out<<" NKnownLevels = "<<NKnownLevels<<std::endl;
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 %10.4g %3d %3d",theKnownLevels[i].id+1,theKnownLevels[i].Energy,theKnownLevels[i].spinx2/2.,2*(G4int)theKnownLevels[i].parity-1,theKnownLevels[i].T12,theKnownLevels[i].NGammas,theKnownLevels[i].Ndecays);
2082 out<<buffer;
2083 for(G4int j=0;j<theKnownLevels[i].Ndecays;j++){
2084 snprintf(buffer,1000," %10.4g %7s",theKnownLevels[i].decayFraction[j],theKnownLevels[i].decayMode[j].c_str());
2085 out<<buffer;
2086 }
2087 out<<std::endl;
2088 for(G4int j=0;j<theKnownLevels[i].NGammas;j++){
2089 snprintf(buffer,1000," %4d %10.4g %10.4g %10.4g %10.4g %10.4g %2d",theKnownLevels[i].FinalLevelID[j]+1,theKnownLevels[i].Eg[j],theKnownLevels[i].Pg[j],theKnownLevels[i].Pe[j],theKnownLevels[i].Icc[j],theKnownLevels[i].cumulPtot[j],theKnownLevels[i].multipolarity[j]);
2090 out<<buffer<<std::endl;
2091 }
2092 }
2093 out<<" ########################################################################################################## "<<std::endl;
2094
2095 }
2096
2097 void G4NuDEXStatisticalNucleus::PrintKnownLevelsInDEGENformat(std::ostream &out){
2098
2099 out<<" ########################################################################################################## "<<std::endl;
2100 out<<" KNOWN_LEVES_DEGEN "<<std::endl;
2101 out<<" NKnownLevels = "<<NKnownLevels<<std::endl;
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;j++){
2108 if(theKnownLevels[i].Pg[j]>MaxIntens){MaxIntens=theKnownLevels[i].Pg[j];}
2109 }
2110 for(G4int j=0;j<theKnownLevels[i].NGammas;j++){
2111 //snprintf(buffer,1000,"%10.3f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f",theKnownLevels[i].Energy*1000.,theKnownLevels[i].spinx2/2.,2.*(G4int)theKnownLevels[i].parity-1,theKnownLevels[i].Eg[j]*1000.,0.,theKnownLevels[i].Pg[j]/MaxIntens*100.,0.,theKnownLevels[i].Icc[j]);
2112 GammaEnergy=theKnownLevels[i].Energy-theKnownLevels[theKnownLevels[i].FinalLevelID[j]].Energy;
2113 snprintf(buffer,1000,"%10.3f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f",theKnownLevels[i].Energy*1000.,theKnownLevels[i].spinx2/2.,2.*(G4int)theKnownLevels[i].parity-1,GammaEnergy*1000.,0.,theKnownLevels[i].Pg[j]/MaxIntens*100.,0.,theKnownLevels[i].Icc[j]);
2114 out<<buffer<<std::endl;
2115 }
2116 }
2117 out<<" ########################################################################################################## "<<std::endl;
2118
2119 }
2120
2121 void G4NuDEXStatisticalNucleus::PrintLevelDensity(std::ostream &out){
2122
2123 if(theLD==0){return;}
2124
2125 G4double Emin=0;
2126 G4double Emax=E_unk_max;
2127 G4int np=100;
2128
2129 out<<" ###################################################################################### "<<std::endl;
2130 out<<" LEVELDENSITY"<<std::endl;
2131 G4double ene,exp=0;
2132 G4double *ld=new G4double[maxspinx2+2];
2133 G4bool *WriteThisSpin=new G4bool[maxspinx2+1];
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<<" "<<Ecrit<<" "<<maxspinx2<<std::endl;
2143 out<<"ENE EXP TOT SUM(J)";
2144 for(G4int spx2=0;spx2<=maxspinx2;spx2++){
2145 if(WriteThisSpin[spx2]){out<<" J="<<spx2/2.;}
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[j].Energy<ene){exp+=theLevels[j].NLevels;}}
2153 out<<ene<<" "<<exp<<" ";
2154 ld[maxspinx2+1]=0;
2155 for(G4int spx2=0;spx2<=maxspinx2;spx2++){
2156 ld[spx2]=2*theLD->GetLevelDensity(ene,spx2/2.,true);
2157 ld[maxspinx2+1]+=ld[spx2];
2158 }
2159 //out<<ld[maxspinx2+1];
2160 out<<theLD->GetLevelDensity(ene,0,true,true)<<" "<<ld[maxspinx2+1];
2161 for(G4int spx2=0;spx2<=maxspinx2;spx2++){
2162 if(WriteThisSpin[spx2]){out<<" "<<ld[spx2];}
2163 }
2164 out<<std::endl;
2165 }
2166 out<<" ###################################################################################### "<<std::endl;
2167
2168 delete [] ld;
2169 delete [] WriteThisSpin;
2170 }
2171
2172 void G4NuDEXStatisticalNucleus::PrintLevelSchemeInDEGENformat(const char* fname,G4int MaxLevelID){
2173
2174 std::ofstream out(fname);
2175 char buffer[1000];
2176 for(G4int i=0;i<NLevels;i++){
2177 if(theLevels[i].Energy>Ecrit && (MaxLevelID>0 && i<=MaxLevelID)){
2178 snprintf(buffer,1000,"%13.5f %17.8f %17.8f ",theLevels[i].Energy*1000.,theLevels[i].spinx2/2.,2.*(G4int)theLevels[i].parity-1);
2179 out<<buffer<<std::endl;
2180 }
2181 }
2182 out.close();
2183
2184 }
2185
2186 void G4NuDEXStatisticalNucleus::PrintLevelScheme(std::ostream &out){
2187 out<<" ###################################################################################### "<<std::endl;
2188 out<<" LEVELSCHEME"<<std::endl;
2189 for(G4int i=0;i<NLevels;i++){
2190 out<<i<<" "<<theLevels[i].Energy<<" "<<theLevels[i].spinx2/2.<<" "<<theLevels[i].parity<<" "<<theLevels[i].KnownLevelID<<" "<<theLevels[i].NLevels<<" "<<theLevels[i].Width<<" "<<theLevels[i].seed<<std::endl;
2191 }
2192 out<<" ###################################################################################### "<<std::endl;
2193 }
2194
2195 void G4NuDEXStatisticalNucleus::PrintThermalPrimaryTransitions(std::ostream &out){
2196
2197 out<<" #################################################### "<<std::endl;
2198 out<<" THERMAL PRIMARY TRANSITIONS"<<std::endl;
2199 out<<" "<<NLevelsBelowThermalCaptureLevel<<std::endl;
2200 if(theThermalCaptureLevelCumulBR!=0){
2201 out<<" "<<0<<" "<<theLevels[0].Energy<<" "<<Sn-theLevels[0].Energy<<" "<<theThermalCaptureLevelCumulBR[0]<<std::endl;
2202 for(G4int i=1;i<NLevelsBelowThermalCaptureLevel;i++){
2203 out<<" "<<i<<" "<<theLevels[i].Energy<<" "<<Sn-theLevels[i].Energy<<" "<<theThermalCaptureLevelCumulBR[i]-theThermalCaptureLevelCumulBR[i-1]<<std::endl;
2204 }
2205 }
2206 out<<" #################################################### "<<std::endl;
2207
2208 G4double ThresholdIntensity=0.01;
2209 out<<" #################################################### "<<std::endl;
2210 out<<" STRONGEST THERMAL PRIMARY TRANSITIONS"<<std::endl;
2211 out<<" "<<NLevelsBelowThermalCaptureLevel<<std::endl;
2212 if(theThermalCaptureLevelCumulBR!=0){
2213 if(theThermalCaptureLevelCumulBR[0]>ThresholdIntensity){out<<" "<<0<<" "<<theLevels[0].Energy<<" "<<Sn-theLevels[0].Energy<<" "<<theThermalCaptureLevelCumulBR[0]<<std::endl;}
2214 for(G4int i=1;i<NLevelsBelowThermalCaptureLevel;i++){
2215 if(theThermalCaptureLevelCumulBR[i]-theThermalCaptureLevelCumulBR[i-1]>ThresholdIntensity){out<<" "<<i<<" "<<theLevels[i].Energy<<" "<<Sn-theLevels[i].Energy<<" "<<theThermalCaptureLevelCumulBR[i]-theThermalCaptureLevelCumulBR[i-1]<<std::endl;}
2216 }
2217 }
2218 out<<" #################################################### "<<std::endl;
2219 }
2220
2221 void G4NuDEXStatisticalNucleus::PrintTotalCumulBR(G4int i_level,std::ostream &out){
2222
2223 if(TotalCumulBR[i_level]!=0){
2224 out<<" #################################################### "<<std::endl;
2225 out<<" CUMULBR FROM LEVEL "<<i_level<<" with ENERGY "<<theLevels[i_level].Energy<<std::endl;
2226 for(G4int i=0;i<i_level;i++){
2227 out<<theLevels[i].Energy<<" "<<theLevels[i].spinx2/2.<<" "<<theLevels[i].parity<<" "<<TotalCumulBR[i_level][i]<<std::endl;
2228 }
2229 out<<" #################################################### "<<std::endl;
2230 }
2231
2232 }
2233
2234 void G4NuDEXStatisticalNucleus::PrintBR(G4int i_level,G4double MaxExcEneToPrint_MeV,std::ostream &out){
2235
2236 if(TotalCumulBR[i_level]!=0){
2237 out<<" #################################################### "<<std::endl;
2238 out<<" BR FROM LEVEL "<<i_level<<" with ENERGY "<<theLevels[i_level].Energy<<std::endl;
2239 for(G4int i=0;i<i_level;i++){
2240 if(theLevels[i].Energy<MaxExcEneToPrint_MeV || MaxExcEneToPrint_MeV<0){
2241 if(i==0){
2242 out<<theLevels[i].Energy<<" "<<theLevels[i].spinx2/2.<<" "<<theLevels[i].parity<<" "<<TotalCumulBR[i_level][i]<<std::endl;
2243 }
2244 else{
2245 out<<theLevels[i].Energy<<" "<<theLevels[i].spinx2/2.<<" "<<theLevels[i].parity<<" "<<TotalCumulBR[i_level][i]-TotalCumulBR[i_level][i-1]<<std::endl;
2246 }
2247 }
2248 }
2249 out<<" #################################################### "<<std::endl;
2250 }
2251
2252
2253 }
2254
2255
2256 void G4NuDEXStatisticalNucleus::PrintPSF(std::ostream &out){
2257
2258 thePSF->PrintPSFParameters(out);
2259
2260 G4int NVals=400;
2261 G4int nEnePSF=(G4int)Sn+1; //number of excitation energies where the PSF are evaluated
2262 G4double EnePSF[200];
2263 G4double Emin=0;
2264 G4double Emax=10;
2265 G4double xval,e1,m1,e2;
2266
2267 out<<" #################################################### "<<std::endl;
2268 out<<" PSF"<<std::endl;
2269 out<<" "<<NVals<<" "<<Emin<<" "<<Emax<<" "<<nEnePSF<<std::endl;
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 "<<std::endl;
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.4E %10.4E",xval,e1,m1,e2);
2288 out<<word<<std::endl;
2289 }
2290 }
2291 out<<" #################################################### "<<std::endl;
2292
2293 }
2294
2295
2296 void G4NuDEXStatisticalNucleus::PrintICC(std::ostream &out){
2297
2298 theICC->PrintICC(out);
2299
2300 }
2301
2302 void G4NuDEXStatisticalNucleus::PrintAll(std::ostream &out){
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(std::ostream &out){
2317
2318 out<<"LEVELDENSITYTYPE "<<LevelDensityType<<std::endl;
2319 out<<"MAXSPIN "<<maxspinx2/2.<<std::endl;
2320 out<<"MINLEVELSPERBAND "<<MinLevelsPerBand<<std::endl;
2321 out<<"BANDWIDTH_MEV "<<BandWidth<<std::endl;
2322 out<<"MAXEXCENERGY_MEV "<<MaxExcEnergy<<std::endl;
2323 out<<"ECRIT_MEV "<<Ecrit<<std::endl;
2324 out<<"KNOWNLEVELSFLAG "<<KnownLevelsFlag<<std::endl;
2325 out<<std::endl;
2326 out<<"PSF_FLAG "<<PSFflag<<std::endl;
2327 out<<"BROPTION "<<BROpt<<std::endl;
2328 out<<"SAMPLEGAMMAWIDTHS "<<SampleGammaWidths<<std::endl;
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 "<<ElectronConversionFlag<<std::endl;
2335 out<<"PRIMARYTHCAPGAMNORM "<<PrimaryGammasIntensityNormFactor<<std::endl;
2336 out<<"PRIMARYGAMMASECUT "<<PrimaryGammasEcut<<std::endl;
2337 out<<std::endl;
2338 theLD->PrintParametersInInputFileFormat(out);
2339 thePSF->PrintPSFParametersInInputFileFormat(out);
2340 out<<std::endl;
2341 out<<"END"<<std::endl;
2342
2343 }
2344
2345 G4int ComparisonLevels(const void* va, const void* vb)
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:-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