[819] | 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 | // |
---|
[1347] | 26 | // $Id: G4NuclearLevelManager.cc,v 1.14 2010/11/17 16:50:53 vnivanch Exp $ |
---|
| 27 | // GEANT4 tag $Name: geant4-09-04-ref-00 $ |
---|
| 28 | // |
---|
[819] | 29 | // ------------------------------------------------------------------- |
---|
| 30 | // GEANT 4 class file |
---|
| 31 | // |
---|
| 32 | // CERN, Geneva, Switzerland |
---|
| 33 | // |
---|
| 34 | // File name: G4NuclearLevelManager |
---|
| 35 | // |
---|
| 36 | // Author: Maria Grazia Pia (pia@genova.infn.it) |
---|
| 37 | // |
---|
| 38 | // Creation date: 24 October 1998 |
---|
| 39 | // |
---|
| 40 | // Modifications: |
---|
| 41 | // |
---|
| 42 | // 15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it) |
---|
| 43 | // Added half-life, angular momentum, parity, emissioni type |
---|
| 44 | // reading from experimental data. |
---|
| 45 | // 02 May 2003, Vladimir Ivanchenko remove rublic copy constructor |
---|
[1340] | 46 | // 06 Oct 2010, M. Kelsey -- Use object storage, not pointers, drop |
---|
| 47 | // public access to list, simplify list construction |
---|
[819] | 48 | // ------------------------------------------------------------------- |
---|
| 49 | |
---|
| 50 | #include "G4NuclearLevelManager.hh" |
---|
| 51 | |
---|
| 52 | #include "globals.hh" |
---|
| 53 | #include "G4NuclearLevel.hh" |
---|
| 54 | #include "G4ios.hh" |
---|
| 55 | #include "G4HadronicException.hh" |
---|
| 56 | #include <stdlib.h> |
---|
| 57 | #include <fstream> |
---|
| 58 | #include <sstream> |
---|
| 59 | #include <algorithm> |
---|
| 60 | #include "G4HadTmpUtil.hh" |
---|
| 61 | |
---|
| 62 | G4NuclearLevelManager::G4NuclearLevelManager(): |
---|
| 63 | _nucleusA(0), _nucleusZ(0), _fileName(""), _validity(false), |
---|
| 64 | _levels(0), _levelEnergy(0), _gammaEnergy(0), _probability(0) |
---|
| 65 | { } |
---|
| 66 | |
---|
[1347] | 67 | G4NuclearLevelManager::G4NuclearLevelManager(G4int Z, G4int A, const G4String& filename) : |
---|
[1340] | 68 | _nucleusA(A), _nucleusZ(Z), _fileName(filename), _validity(false), |
---|
| 69 | _levels(0), _levelEnergy(0), _gammaEnergy(0), _probability(0) |
---|
[819] | 70 | { |
---|
| 71 | if (A <= 0 || Z <= 0 || Z > A ) |
---|
| 72 | throw G4HadronicException(__FILE__, __LINE__, "==== G4NuclearLevelManager ==== (Z,A) <0, or Z>A"); |
---|
| 73 | |
---|
| 74 | MakeLevels(); |
---|
| 75 | } |
---|
| 76 | |
---|
| 77 | G4NuclearLevelManager::~G4NuclearLevelManager() |
---|
[1340] | 78 | { |
---|
| 79 | ClearLevels(); |
---|
[819] | 80 | } |
---|
| 81 | |
---|
[1347] | 82 | void G4NuclearLevelManager::SetNucleus(G4int Z, G4int A, const G4String& filename) |
---|
[819] | 83 | { |
---|
[1340] | 84 | if (A <= 0 || Z <= 0 || Z > A ) |
---|
| 85 | throw G4HadronicException(__FILE__, __LINE__, "==== G4NuclearLevelManager ==== (Z,A) <0, or Z>A"); |
---|
| 86 | |
---|
| 87 | if (_nucleusZ != Z || _nucleusA != A) |
---|
[819] | 88 | { |
---|
[1340] | 89 | _nucleusA = A; |
---|
| 90 | _nucleusZ = Z; |
---|
| 91 | _fileName = filename; |
---|
| 92 | MakeLevels(); |
---|
[819] | 93 | } |
---|
| 94 | } |
---|
| 95 | |
---|
[1347] | 96 | const G4NuclearLevel* G4NuclearLevelManager::GetLevel(G4int i) const { |
---|
[1340] | 97 | return (i>=0 && i<NumberOfLevels()) ? (*_levels)[i] : 0; |
---|
[819] | 98 | } |
---|
| 99 | |
---|
| 100 | |
---|
[1340] | 101 | const G4NuclearLevel* |
---|
[1347] | 102 | G4NuclearLevelManager::NearestLevel(G4double energy, |
---|
| 103 | G4double eDiffMax) const |
---|
| 104 | { |
---|
[1340] | 105 | if (NumberOfLevels() <= 0) return 0; |
---|
[819] | 106 | |
---|
[1340] | 107 | G4int iNear = -1; |
---|
[819] | 108 | |
---|
[1340] | 109 | G4double diff = 9999. * GeV; |
---|
| 110 | for (unsigned int i=0; i<_levels->size(); i++) |
---|
[819] | 111 | { |
---|
[1340] | 112 | G4double e = GetLevel(i)->Energy(); |
---|
| 113 | G4double eDiff = std::abs(e - energy); |
---|
| 114 | if (eDiff < diff && eDiff <= eDiffMax) |
---|
| 115 | { |
---|
| 116 | diff = eDiff; |
---|
| 117 | iNear = i; |
---|
[819] | 118 | } |
---|
| 119 | } |
---|
[1340] | 120 | |
---|
| 121 | return GetLevel(iNear); // Includes range checking on iNear |
---|
[819] | 122 | } |
---|
| 123 | |
---|
| 124 | |
---|
| 125 | G4double G4NuclearLevelManager::MinLevelEnergy() const |
---|
| 126 | { |
---|
[1340] | 127 | return (NumberOfLevels() > 0) ? _levels->front()->Energy() : 9999.*GeV; |
---|
[819] | 128 | } |
---|
| 129 | |
---|
| 130 | |
---|
| 131 | G4double G4NuclearLevelManager::MaxLevelEnergy() const |
---|
| 132 | { |
---|
[1340] | 133 | return (NumberOfLevels() > 0) ? _levels->back()->Energy() : 0.*GeV; |
---|
[819] | 134 | } |
---|
| 135 | |
---|
| 136 | |
---|
| 137 | const G4NuclearLevel* G4NuclearLevelManager::HighestLevel() const |
---|
| 138 | { |
---|
[1340] | 139 | return (NumberOfLevels() > 0) ? _levels->front() : 0; |
---|
[819] | 140 | } |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | const G4NuclearLevel* G4NuclearLevelManager::LowestLevel() const |
---|
| 144 | { |
---|
[1340] | 145 | return (NumberOfLevels() > 0) ? _levels->back() : 0; |
---|
[819] | 146 | } |
---|
| 147 | |
---|
| 148 | |
---|
[1347] | 149 | G4bool G4NuclearLevelManager::Read(std::ifstream& dataFile) |
---|
| 150 | { |
---|
[1340] | 151 | G4bool goodRead = ReadDataLine(dataFile); |
---|
[819] | 152 | |
---|
[1340] | 153 | if (goodRead) ProcessDataLine(); |
---|
| 154 | return goodRead; |
---|
| 155 | } |
---|
| 156 | |
---|
| 157 | // NOTE: Standard stream I/O generates a 45 byte std::string per item! |
---|
| 158 | |
---|
| 159 | G4bool G4NuclearLevelManager::ReadDataLine(std::ifstream& dataFile) { |
---|
| 160 | /***** DO NOT USE REGULAR STREAM I/O |
---|
[819] | 161 | G4bool result = true; |
---|
| 162 | if (dataFile >> _levelEnergy) |
---|
| 163 | { |
---|
| 164 | dataFile >> _gammaEnergy >> _probability >> _polarity >> _halfLife |
---|
| 165 | >> _angularMomentum >> _totalCC >> _kCC >> _l1CC >> _l2CC |
---|
| 166 | >> _l3CC >> _m1CC >> _m2CC >> _m3CC >> _m4CC >> _m5CC |
---|
| 167 | >> _nPlusCC; |
---|
| 168 | } |
---|
[1340] | 169 | else result = false; |
---|
| 170 | *****/ |
---|
| 171 | |
---|
| 172 | // Each item will return iostream status |
---|
| 173 | return (ReadDataItem(dataFile, _levelEnergy) && |
---|
| 174 | ReadDataItem(dataFile, _gammaEnergy) && |
---|
| 175 | ReadDataItem(dataFile, _probability) && |
---|
| 176 | ReadDataItem(dataFile, _polarity) && |
---|
| 177 | ReadDataItem(dataFile, _halfLife) && |
---|
| 178 | ReadDataItem(dataFile, _angularMomentum) && |
---|
| 179 | ReadDataItem(dataFile, _totalCC) && |
---|
| 180 | ReadDataItem(dataFile, _kCC) && |
---|
| 181 | ReadDataItem(dataFile, _l1CC) && |
---|
| 182 | ReadDataItem(dataFile, _l2CC) && |
---|
| 183 | ReadDataItem(dataFile, _l3CC) && |
---|
| 184 | ReadDataItem(dataFile, _m1CC) && |
---|
| 185 | ReadDataItem(dataFile, _m2CC) && |
---|
| 186 | ReadDataItem(dataFile, _m3CC) && |
---|
| 187 | ReadDataItem(dataFile, _m4CC) && |
---|
| 188 | ReadDataItem(dataFile, _m5CC) && |
---|
| 189 | ReadDataItem(dataFile, _nPlusCC) ); |
---|
[819] | 190 | } |
---|
| 191 | |
---|
[1340] | 192 | G4bool |
---|
[1347] | 193 | G4NuclearLevelManager::ReadDataItem(std::istream& dataFile, G4double& x) |
---|
| 194 | { |
---|
[1340] | 195 | G4bool okay = (dataFile >> buffer); // Get next token |
---|
| 196 | if (okay) x = strtod(buffer, NULL); |
---|
[819] | 197 | |
---|
[1340] | 198 | return okay; |
---|
| 199 | } |
---|
| 200 | |
---|
[1347] | 201 | void G4NuclearLevelManager::ProcessDataLine() |
---|
| 202 | { |
---|
[1340] | 203 | const G4double minProbability = 1e-8; |
---|
| 204 | |
---|
| 205 | // Assign units for dimensional quantities |
---|
| 206 | _levelEnergy *= keV; |
---|
| 207 | _gammaEnergy *= keV; |
---|
| 208 | _halfLife *= second; |
---|
| 209 | |
---|
| 210 | // The following adjustment is needed to take care of anomalies in |
---|
| 211 | // data files, where some transitions show up with relative probability |
---|
| 212 | // zero |
---|
| 213 | if (_probability < minProbability) _probability = minProbability; |
---|
| 214 | // the folowwing is to convert icc probability to accumulative ones |
---|
| 215 | _l1CC += _kCC; |
---|
| 216 | _l2CC += _l1CC; |
---|
| 217 | _l3CC += _l2CC; |
---|
| 218 | _m1CC += _l3CC; |
---|
| 219 | _m2CC += _m1CC; |
---|
| 220 | _m3CC += _m2CC; |
---|
| 221 | _m4CC += _m3CC; |
---|
| 222 | _m5CC += _m4CC; |
---|
| 223 | _nPlusCC += _m5CC; |
---|
| 224 | |
---|
| 225 | if (_nPlusCC!=0) { // Normalize to probabilities |
---|
| 226 | _kCC /= _nPlusCC; |
---|
| 227 | _l1CC /= _nPlusCC; |
---|
| 228 | _l2CC /= _nPlusCC; |
---|
| 229 | _l3CC /= _nPlusCC; |
---|
| 230 | _m1CC /= _nPlusCC; |
---|
| 231 | _m2CC /= _nPlusCC; |
---|
| 232 | _m3CC /= _nPlusCC; |
---|
| 233 | _m4CC /= _nPlusCC; |
---|
| 234 | _m5CC /= _nPlusCC; |
---|
| 235 | _nPlusCC /= _nPlusCC; |
---|
| 236 | } else { // Total was zero, reset to unity |
---|
| 237 | _kCC = 1; |
---|
| 238 | _l1CC = 1; |
---|
| 239 | _l2CC = 1; |
---|
| 240 | _l3CC = 1; |
---|
| 241 | _m1CC = 1; |
---|
| 242 | _m2CC = 1; |
---|
| 243 | _m3CC = 1; |
---|
| 244 | _m4CC = 1; |
---|
| 245 | _m5CC = 1; |
---|
| 246 | _nPlusCC = 1; |
---|
| 247 | } |
---|
| 248 | |
---|
| 249 | // G4cout << "Read " << _levelEnergy << " " << _gammaEnergy << " " << _probability << G4endl; |
---|
| 250 | } |
---|
| 251 | |
---|
| 252 | |
---|
| 253 | void G4NuclearLevelManager::ClearLevels() |
---|
| 254 | { |
---|
| 255 | _validity = false; |
---|
| 256 | |
---|
| 257 | if (NumberOfLevels() > 0) { |
---|
| 258 | std::for_each(_levels->begin(), _levels->end(), DeleteLevel()); |
---|
| 259 | _levels->clear(); |
---|
| 260 | } |
---|
| 261 | |
---|
| 262 | delete _levels; |
---|
| 263 | _levels = 0; |
---|
| 264 | } |
---|
| 265 | |
---|
[819] | 266 | void G4NuclearLevelManager::MakeLevels() |
---|
| 267 | { |
---|
| 268 | _validity = false; |
---|
[1340] | 269 | if (NumberOfLevels() > 0) ClearLevels(); // Discard existing data |
---|
| 270 | |
---|
[819] | 271 | std::ifstream inFile(_fileName, std::ios::in); |
---|
| 272 | if (! inFile) |
---|
| 273 | { |
---|
| 274 | #ifdef GAMMAFILEWARNING |
---|
| 275 | if (_nucleusZ > 10) G4cout << " G4NuclearLevelManager: nuclide (" |
---|
| 276 | << _nucleusZ << "," << _nucleusA |
---|
| 277 | << ") does not have a gamma levels file" << G4endl; |
---|
| 278 | #endif |
---|
| 279 | return; |
---|
| 280 | } |
---|
| 281 | |
---|
| 282 | _levels = new G4PtrLevelVector; |
---|
[1340] | 283 | |
---|
| 284 | // Read individual gamma data and fill levels for this nucleus |
---|
[819] | 285 | |
---|
[1340] | 286 | G4NuclearLevel* thisLevel = 0; |
---|
| 287 | G4int nData = 0; |
---|
| 288 | |
---|
| 289 | while (Read(inFile)) { |
---|
| 290 | thisLevel = UseLevelOrMakeNew(thisLevel); // May create new pointer |
---|
| 291 | AddDataToLevel(thisLevel); |
---|
| 292 | nData++; // For debugging purposes |
---|
| 293 | } |
---|
| 294 | |
---|
| 295 | FinishLevel(thisLevel); // Final must be completed by hand |
---|
[819] | 296 | |
---|
| 297 | // ---- MGP ---- Don't forget to close the file |
---|
| 298 | inFile.close(); |
---|
| 299 | |
---|
| 300 | // G4cout << " ==== MakeLevels ===== " << nData << " data read " << G4endl; |
---|
| 301 | |
---|
[1340] | 302 | G4PtrSort<G4NuclearLevel>(_levels); |
---|
| 303 | |
---|
| 304 | _validity = true; |
---|
| 305 | |
---|
| 306 | return; |
---|
| 307 | } |
---|
[819] | 308 | |
---|
[1340] | 309 | G4NuclearLevel* |
---|
[1347] | 310 | G4NuclearLevelManager::UseLevelOrMakeNew(G4NuclearLevel* level) |
---|
| 311 | { |
---|
[1340] | 312 | if (level && _levelEnergy == level->Energy()) return level; // No change |
---|
[819] | 313 | |
---|
[1340] | 314 | if (level) FinishLevel(level); // Save what we have up to now |
---|
| 315 | |
---|
| 316 | // G4cout << "Making a new level... " << _levelEnergy << G4endl; |
---|
| 317 | return new G4NuclearLevel(_levelEnergy, _halfLife, _angularMomentum); |
---|
[819] | 318 | } |
---|
| 319 | |
---|
[1347] | 320 | void G4NuclearLevelManager::AddDataToLevel(G4NuclearLevel* level) |
---|
| 321 | { |
---|
[1340] | 322 | if (!level) return; // Sanity check |
---|
[819] | 323 | |
---|
[1340] | 324 | level->_energies.push_back(_gammaEnergy); |
---|
| 325 | level->_weights.push_back(_probability); |
---|
| 326 | level->_polarities.push_back(_polarity); |
---|
| 327 | level->_kCC.push_back(_kCC); |
---|
| 328 | level->_l1CC.push_back(_l1CC); |
---|
| 329 | level->_l2CC.push_back(_l2CC); |
---|
| 330 | level->_l3CC.push_back(_l3CC); |
---|
| 331 | level->_m1CC.push_back(_m1CC); |
---|
| 332 | level->_m2CC.push_back(_m2CC); |
---|
| 333 | level->_m3CC.push_back(_m3CC); |
---|
| 334 | level->_m4CC.push_back(_m4CC); |
---|
| 335 | level->_m5CC.push_back(_m5CC); |
---|
| 336 | level->_nPlusCC.push_back(_nPlusCC); |
---|
| 337 | level->_totalCC.push_back(_totalCC); |
---|
| 338 | } |
---|
| 339 | |
---|
[1347] | 340 | void G4NuclearLevelManager::FinishLevel(G4NuclearLevel* level) |
---|
| 341 | { |
---|
[1340] | 342 | if (!level || !_levels) return; // Sanity check |
---|
| 343 | |
---|
| 344 | level->Finalize(); |
---|
| 345 | _levels->push_back(level); |
---|
| 346 | } |
---|
| 347 | |
---|
| 348 | |
---|
[819] | 349 | void G4NuclearLevelManager::PrintAll() |
---|
| 350 | { |
---|
[1340] | 351 | G4int nLevels = NumberOfLevels(); |
---|
[819] | 352 | |
---|
[1340] | 353 | G4cout << " ==== G4NuclearLevelManager ==== (" << _nucleusZ << ", " << _nucleusA |
---|
| 354 | << ") has " << nLevels << " levels" << G4endl |
---|
| 355 | << "Highest level is at energy " << MaxLevelEnergy() << " MeV " |
---|
| 356 | << G4endl << "Lowest level is at energy " << MinLevelEnergy() |
---|
| 357 | << " MeV " << G4endl; |
---|
[819] | 358 | |
---|
[1340] | 359 | for (G4int i=0; i<nLevels; i++) |
---|
| 360 | GetLevel(i)->PrintAll(); |
---|
[819] | 361 | } |
---|
| 362 | |
---|
| 363 | G4NuclearLevelManager::G4NuclearLevelManager(const G4NuclearLevelManager &right) |
---|
| 364 | { |
---|
| 365 | _levelEnergy = right._levelEnergy; |
---|
| 366 | _gammaEnergy = right._gammaEnergy; |
---|
| 367 | _probability = right._probability; |
---|
| 368 | _polarity = right._polarity; |
---|
| 369 | _halfLife = right._halfLife; |
---|
| 370 | _angularMomentum = right._angularMomentum; |
---|
| 371 | _kCC = right._kCC; |
---|
| 372 | _l1CC = right._l1CC; |
---|
| 373 | _l2CC = right._l2CC; |
---|
| 374 | _l3CC = right._l3CC; |
---|
| 375 | _m1CC = right._m1CC; |
---|
| 376 | _m2CC = right._m2CC; |
---|
| 377 | _m3CC = right._m3CC; |
---|
| 378 | _m4CC = right._m4CC; |
---|
| 379 | _m5CC = right._m5CC; |
---|
| 380 | _nPlusCC = right._nPlusCC; |
---|
| 381 | _totalCC = right._totalCC; |
---|
| 382 | _nucleusA = right._nucleusA; |
---|
| 383 | _nucleusZ = right._nucleusZ; |
---|
| 384 | _fileName = right._fileName; |
---|
| 385 | _validity = right._validity; |
---|
[1340] | 386 | |
---|
[819] | 387 | if (right._levels != 0) |
---|
| 388 | { |
---|
| 389 | _levels = new G4PtrLevelVector; |
---|
| 390 | G4int n = right._levels->size(); |
---|
| 391 | G4int i; |
---|
| 392 | for (i=0; i<n; i++) |
---|
| 393 | { |
---|
[1340] | 394 | _levels->push_back(new G4NuclearLevel(*(right.GetLevel(i)))); |
---|
[819] | 395 | } |
---|
| 396 | |
---|
| 397 | G4PtrSort<G4NuclearLevel>(_levels); |
---|
| 398 | } |
---|
| 399 | else |
---|
| 400 | { |
---|
| 401 | _levels = 0; |
---|
| 402 | } |
---|
| 403 | } |
---|
| 404 | |
---|