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
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048 #include <stdlib.h>
00049 #include <fstream>
00050 #include <sstream>
00051 #include <algorithm>
00052
00053 #include "G4NuclearLevelManager.hh"
00054
00055 #include "globals.hh"
00056 #include "G4SystemOfUnits.hh"
00057 #include "G4NuclearLevel.hh"
00058 #include "G4ios.hh"
00059 #include "G4HadronicException.hh"
00060 #include "G4HadTmpUtil.hh"
00061
00062
00063
00064
00065
00066
00067 G4NuclearLevelManager::G4NuclearLevelManager(G4int Z, G4int A, const G4String& filename) :
00068 _nucleusA(A), _nucleusZ(Z), _fileName(filename), _validity(false),
00069 _levels(0), _levelEnergy(0), _gammaEnergy(0), _probability(0)
00070 {
00071 if (A <= 0 || Z <= 0 || Z > A )
00072 throw G4HadronicException(__FILE__, __LINE__, "==== G4NuclearLevelManager ==== (Z,A) <0, or Z>A");
00073
00074 MakeLevels();
00075 }
00076
00077 G4NuclearLevelManager::~G4NuclearLevelManager()
00078 {
00079 ClearLevels();
00080 }
00081
00082 void G4NuclearLevelManager::SetNucleus(G4int Z, G4int A, const G4String& filename)
00083 {
00084 if (A <= 0 || Z <= 0 || Z > A )
00085 throw G4HadronicException(__FILE__, __LINE__, "==== G4NuclearLevelManager ==== (Z,A) <0, or Z>A");
00086
00087 if (_nucleusZ != Z || _nucleusA != A)
00088 {
00089 _nucleusA = A;
00090 _nucleusZ = Z;
00091 _fileName = filename;
00092 MakeLevels();
00093 }
00094 }
00095
00096 const G4NuclearLevel* G4NuclearLevelManager::GetLevel(G4int i) const {
00097 return (i>=0 && i<NumberOfLevels()) ? (*_levels)[i] : 0;
00098 }
00099
00100
00101 const G4NuclearLevel*
00102 G4NuclearLevelManager::NearestLevel(G4double energy,
00103 G4double eDiffMax) const
00104 {
00105 if (NumberOfLevels() <= 0) return 0;
00106
00107 G4int iNear = -1;
00108
00109
00110
00111
00112 G4double diff = 9999. * GeV;
00113 for (unsigned int i=0; i<_levels->size(); i++)
00114 {
00115 G4double e = GetLevel(i)->Energy();
00116 G4double eDiff = std::abs(e - energy);
00117
00118 if (eDiff < diff && eDiff <= eDiffMax)
00119 {
00120 diff = eDiff;
00121 iNear = i;
00122 }
00123 }
00124
00125 return GetLevel(iNear);
00126 }
00127
00128
00129 G4double G4NuclearLevelManager::MinLevelEnergy() const
00130 {
00131 return (NumberOfLevels() > 0) ? _levels->front()->Energy() : 9999.*GeV;
00132 }
00133
00134
00135 G4double G4NuclearLevelManager::MaxLevelEnergy() const
00136 {
00137 return (NumberOfLevels() > 0) ? _levels->back()->Energy() : 0.*GeV;
00138 }
00139
00140
00141 const G4NuclearLevel* G4NuclearLevelManager::HighestLevel() const
00142 {
00143 return (NumberOfLevels() > 0) ? _levels->front() : 0;
00144 }
00145
00146
00147 const G4NuclearLevel* G4NuclearLevelManager::LowestLevel() const
00148 {
00149 return (NumberOfLevels() > 0) ? _levels->back() : 0;
00150 }
00151
00152
00153 G4bool G4NuclearLevelManager::Read(std::ifstream& dataFile)
00154 {
00155 G4bool goodRead = ReadDataLine(dataFile);
00156
00157 if (goodRead) ProcessDataLine();
00158 return goodRead;
00159 }
00160
00161
00162
00163 G4bool G4NuclearLevelManager::ReadDataLine(std::ifstream& dataFile) {
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177 return (ReadDataItem(dataFile, _levelEnergy) &&
00178 ReadDataItem(dataFile, _gammaEnergy) &&
00179 ReadDataItem(dataFile, _probability) &&
00180 ReadDataItem(dataFile, _polarity) &&
00181 ReadDataItem(dataFile, _halfLife) &&
00182 ReadDataItem(dataFile, _angularMomentum) &&
00183 ReadDataItem(dataFile, _totalCC) &&
00184 ReadDataItem(dataFile, _kCC) &&
00185 ReadDataItem(dataFile, _l1CC) &&
00186 ReadDataItem(dataFile, _l2CC) &&
00187 ReadDataItem(dataFile, _l3CC) &&
00188 ReadDataItem(dataFile, _m1CC) &&
00189 ReadDataItem(dataFile, _m2CC) &&
00190 ReadDataItem(dataFile, _m3CC) &&
00191 ReadDataItem(dataFile, _m4CC) &&
00192 ReadDataItem(dataFile, _m5CC) &&
00193 ReadDataItem(dataFile, _nPlusCC) );
00194 }
00195
00196 G4bool
00197 G4NuclearLevelManager::ReadDataItem(std::istream& dataFile, G4double& x)
00198 {
00199 G4bool okay = (dataFile >> buffer);
00200 if (okay) x = strtod(buffer, NULL);
00201
00202 return okay;
00203 }
00204
00205 void G4NuclearLevelManager::ProcessDataLine()
00206 {
00207 const G4double minProbability = 1e-8;
00208
00209
00210 _levelEnergy *= keV;
00211 _gammaEnergy *= keV;
00212 _halfLife *= second;
00213
00214
00215
00216
00217 if (_probability < minProbability) _probability = minProbability;
00218
00219 _l1CC += _kCC;
00220 _l2CC += _l1CC;
00221 _l3CC += _l2CC;
00222 _m1CC += _l3CC;
00223 _m2CC += _m1CC;
00224 _m3CC += _m2CC;
00225 _m4CC += _m3CC;
00226 _m5CC += _m4CC;
00227 _nPlusCC += _m5CC;
00228
00229 if (_nPlusCC!=0) {
00230 _kCC /= _nPlusCC;
00231 _l1CC /= _nPlusCC;
00232 _l2CC /= _nPlusCC;
00233 _l3CC /= _nPlusCC;
00234 _m1CC /= _nPlusCC;
00235 _m2CC /= _nPlusCC;
00236 _m3CC /= _nPlusCC;
00237 _m4CC /= _nPlusCC;
00238 _m5CC /= _nPlusCC;
00239 _nPlusCC /= _nPlusCC;
00240 } else {
00241 _kCC = 1;
00242 _l1CC = 1;
00243 _l2CC = 1;
00244 _l3CC = 1;
00245 _m1CC = 1;
00246 _m2CC = 1;
00247 _m3CC = 1;
00248 _m4CC = 1;
00249 _m5CC = 1;
00250 _nPlusCC = 1;
00251 }
00252
00253
00254 }
00255
00256
00257 void G4NuclearLevelManager::ClearLevels()
00258 {
00259 _validity = false;
00260
00261 if (NumberOfLevels() > 0) {
00262 std::for_each(_levels->begin(), _levels->end(), DeleteLevel());
00263 _levels->clear();
00264 }
00265
00266 delete _levels;
00267 _levels = 0;
00268 }
00269
00270 void G4NuclearLevelManager::MakeLevels()
00271 {
00272 _validity = false;
00273 if (NumberOfLevels() > 0) ClearLevels();
00274
00275 std::ifstream inFile(_fileName, std::ios::in);
00276 if (! inFile)
00277 {
00278 #ifdef GAMMAFILEWARNING
00279 if (_nucleusZ > 10) G4cout << " G4NuclearLevelManager: nuclide ("
00280 << _nucleusZ << "," << _nucleusA
00281 << ") does not have a gamma levels file" << G4endl;
00282 #endif
00283 return;
00284 }
00285
00286 _levels = new G4PtrLevelVector;
00287
00288
00289
00290 G4NuclearLevel* thisLevel = 0;
00291 G4int nData = 0;
00292
00293 while (Read(inFile)) {
00294 thisLevel = UseLevelOrMakeNew(thisLevel);
00295 AddDataToLevel(thisLevel);
00296 nData++;
00297 }
00298
00299 FinishLevel(thisLevel);
00300
00301
00302 inFile.close();
00303
00304
00305
00306 G4PtrSort<G4NuclearLevel>(_levels);
00307
00308 _validity = true;
00309
00310 return;
00311 }
00312
00313 G4NuclearLevel*
00314 G4NuclearLevelManager::UseLevelOrMakeNew(G4NuclearLevel* level)
00315 {
00316 if (level && _levelEnergy == level->Energy()) return level;
00317
00318 if (level) FinishLevel(level);
00319
00320
00321 return new G4NuclearLevel(_levelEnergy, _halfLife, _angularMomentum);
00322 }
00323
00324 void G4NuclearLevelManager::AddDataToLevel(G4NuclearLevel* level)
00325 {
00326 if (!level) return;
00327
00328 level->_energies.push_back(_gammaEnergy);
00329 level->_weights.push_back(_probability);
00330 level->_polarities.push_back(_polarity);
00331 level->_kCC.push_back(_kCC);
00332 level->_l1CC.push_back(_l1CC);
00333 level->_l2CC.push_back(_l2CC);
00334 level->_l3CC.push_back(_l3CC);
00335 level->_m1CC.push_back(_m1CC);
00336 level->_m2CC.push_back(_m2CC);
00337 level->_m3CC.push_back(_m3CC);
00338 level->_m4CC.push_back(_m4CC);
00339 level->_m5CC.push_back(_m5CC);
00340 level->_nPlusCC.push_back(_nPlusCC);
00341 level->_totalCC.push_back(_totalCC);
00342 }
00343
00344 void G4NuclearLevelManager::FinishLevel(G4NuclearLevel* level)
00345 {
00346 if (!level || !_levels) return;
00347
00348 level->Finalize();
00349 _levels->push_back(level);
00350 }
00351
00352
00353 void G4NuclearLevelManager::PrintAll()
00354 {
00355 G4int nLevels = NumberOfLevels();
00356
00357 G4cout << " ==== G4NuclearLevelManager ==== (" << _nucleusZ << ", " << _nucleusA
00358 << ") has " << nLevels << " levels" << G4endl
00359 << "Highest level is at energy " << MaxLevelEnergy() << " MeV "
00360 << G4endl << "Lowest level is at energy " << MinLevelEnergy()
00361 << " MeV " << G4endl;
00362
00363 for (G4int i=0; i<nLevels; i++)
00364 GetLevel(i)->PrintAll();
00365 }
00366
00367 G4NuclearLevelManager::G4NuclearLevelManager(const G4NuclearLevelManager &right)
00368 {
00369 _levelEnergy = right._levelEnergy;
00370 _gammaEnergy = right._gammaEnergy;
00371 _probability = right._probability;
00372 _polarity = right._polarity;
00373 _halfLife = right._halfLife;
00374 _angularMomentum = right._angularMomentum;
00375 _kCC = right._kCC;
00376 _l1CC = right._l1CC;
00377 _l2CC = right._l2CC;
00378 _l3CC = right._l3CC;
00379 _m1CC = right._m1CC;
00380 _m2CC = right._m2CC;
00381 _m3CC = right._m3CC;
00382 _m4CC = right._m4CC;
00383 _m5CC = right._m5CC;
00384 _nPlusCC = right._nPlusCC;
00385 _totalCC = right._totalCC;
00386 _nucleusA = right._nucleusA;
00387 _nucleusZ = right._nucleusZ;
00388 _fileName = right._fileName;
00389 _validity = right._validity;
00390
00391 if (right._levels != 0)
00392 {
00393 _levels = new G4PtrLevelVector;
00394 G4int n = right._levels->size();
00395 G4int i;
00396 for (i=0; i<n; ++i)
00397 {
00398 _levels->push_back(new G4NuclearLevel(*(right.GetLevel(i))));
00399 }
00400
00401 G4PtrSort<G4NuclearLevel>(_levels);
00402 }
00403 else
00404 {
00405 _levels = 0;
00406 }
00407 for (G4int i=0; i<30; ++i) {
00408 buffer[i] = right.buffer[i];
00409 }
00410 }
00411