00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "G4NeutronHPDeExGammas.hh"
00033 #include "G4SystemOfUnits.hh"
00034
00035 void G4NeutronHPDeExGammas::Init(std::ifstream & aDataFile)
00036 {
00037 G4NeutronHPGamma ** theGammas = new G4NeutronHPGamma * [50];
00038 G4int nGammas = 0;
00039 G4int nBuff = 50;
00040 for(;;)
00041 {
00042 G4NeutronHPGamma * theNew = new G4NeutronHPGamma;
00043 if(!theNew->Init(aDataFile))
00044 {
00045 delete theNew;
00046 break;
00047 }
00048 else
00049 {
00050 if(nGammas==nBuff)
00051 {
00052 nBuff+=50;
00053 G4NeutronHPGamma ** buffer = new G4NeutronHPGamma * [nBuff];
00054 for(G4int i=0;i<nGammas;i++) buffer[i] = theGammas[i];
00055 delete [] theGammas;
00056 theGammas = buffer;
00057 }
00058 theGammas[nGammas] = theNew;
00059 nGammas++;
00060 }
00061 }
00062
00063
00064
00065
00066 G4double currentE = 0;
00067 G4double nextE = 0;
00068 G4int i;
00069 G4double epsilon = 0.01*keV;
00070 for(i=0; i<nGammas; i++)
00071 {
00072 nextE = theGammas[i]->GetLevelEnergy();
00073 if(std::abs(currentE-nextE)>epsilon) nLevels++;
00074 currentE = nextE;
00075 }
00076
00077
00078
00079 theLevels = new G4NeutronHPLevel[nLevels];
00080 levelStart = new G4int [nLevels];
00081 levelSize = new G4int [nLevels];
00082
00083
00084
00085 currentE = 0;
00086 nextE = 0;
00087 G4int levelCounter=-1;
00088 for(i=0; i<nGammas; i++)
00089 {
00090 nextE = theGammas[i]->GetLevelEnergy();
00091 if(std::abs(currentE-nextE)>epsilon)
00092 {
00093 levelCounter++;
00094 levelStart[levelCounter] = i;
00095 levelSize[levelCounter] = 0;
00096 }
00097 levelSize[levelCounter]++;
00098 currentE = nextE;
00099 }
00100
00101 for(i=0; i<nLevels; i++)
00102 {
00103 theLevels[i].SetNumberOfGammas(levelSize[i]);
00104 for(G4int ii=levelStart[i]; ii<levelStart[i]+levelSize[i]; ii++)
00105 {
00106 theLevels[i].SetGamma(ii-levelStart[i], theGammas[ii]);
00107 }
00108 }
00109
00110
00111 G4double levelE, gammaE, currentLevelE;
00112 G4double min;
00113 for(i=0; i<nGammas; i++)
00114 {
00115 G4int it=-1;
00116 gammaE = theGammas[i]->GetGammaEnergy();
00117 currentLevelE = theGammas[i]->GetLevelEnergy();
00118 min = currentLevelE-gammaE-epsilon;
00119 for(G4int ii=0; ii<nLevels; ii++)
00120 {
00121 levelE = theLevels[ii].GetLevelEnergy();
00122 if(std::abs(currentLevelE-(levelE+gammaE))<min)
00123 {
00124 min = std::abs(currentLevelE-(levelE+gammaE));
00125 it = ii;
00126 }
00127 }
00128
00129 if ( it != -1 && currentLevelE == theLevels[it].GetLevelEnergy() )
00130 {
00131
00132
00133
00134 it +=-1;
00135 }
00136
00137 if(it!=-1) theGammas[i]->SetNext(&theLevels[it]);
00138 }
00139
00140
00141 delete [] theGammas;
00142
00143
00144 }