G4NuclearLevelManager.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //      GEANT 4 class file 
00030 //
00031 //      CERN, Geneva, Switzerland
00032 //
00033 //      File name:     G4NuclearLevelManager
00034 //
00035 //      Author:        Maria Grazia Pia (pia@genova.infn.it)
00036 // 
00037 //      Creation date: 24 October 1998
00038 //
00039 //      Modifications: 
00040 //      
00041 //        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
00042 //              Added half-life, angular momentum, parity, emissioni type
00043 //              reading from experimental data. 
00044 //        02 May 2003,   Vladimir Ivanchenko remove rublic copy constructor
00045 //        06 Oct 2010, M. Kelsey -- Use object storage, not pointers, drop
00046 //              public access to list, simplify list construction
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 G4NuclearLevelManager::G4NuclearLevelManager():
00063     _nucleusA(0), _nucleusZ(0), _fileName(""), _validity(false), 
00064     _levels(0), _levelEnergy(0), _gammaEnergy(0), _probability(0)
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   //G4cout << "G4NuclearLevelManager::NearestLevel E(MeV)= " 
00110   //     << energy/MeV << " dEmax(MeV)= " << eDiffMax/MeV << G4endl;
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       //G4cout << i << ".   eDiff(MeV)= " << eDiff/MeV << G4endl;
00118       if (eDiff < diff && eDiff <= eDiffMax)
00119         { 
00120           diff = eDiff; 
00121           iNear = i;
00122         }
00123     }
00124   
00125   return GetLevel(iNear);       // Includes range checking on 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 // NOTE:  Standard stream I/O generates a 45 byte std::string per item!
00162 
00163 G4bool G4NuclearLevelManager::ReadDataLine(std::ifstream& dataFile) {
00164   /***** DO NOT USE REGULAR STREAM I/O
00165   G4bool result = true;
00166   if (dataFile >> _levelEnergy)
00167     {
00168       dataFile >> _gammaEnergy >> _probability >> _polarity >> _halfLife
00169                >> _angularMomentum  >> _totalCC >> _kCC >> _l1CC >> _l2CC 
00170                >> _l3CC >> _m1CC >> _m2CC >> _m3CC >> _m4CC >> _m5CC
00171                >> _nPlusCC;
00172     }
00173   else result = false;
00174   *****/
00175 
00176   // Each item will return iostream status
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);           // Get next token
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   // Assign units for dimensional quantities
00210   _levelEnergy *= keV;
00211   _gammaEnergy *= keV;
00212   _halfLife *= second;
00213   
00214   // The following adjustment is needed to take care of anomalies in 
00215   // data files, where some transitions show up with relative probability
00216   // zero
00217   if (_probability < minProbability) _probability = minProbability;
00218   // the folowwing is to convert icc probability to accumulative ones
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) {    // Normalize to probabilities
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 {              // Total was zero, reset to unity
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   // G4cout << "Read " << _levelEnergy << " " << _gammaEnergy << " " << _probability << G4endl;
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();      // Discard existing data
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   // Read individual gamma data and fill levels for this nucleus
00289  
00290   G4NuclearLevel* thisLevel = 0;
00291   G4int nData = 0;
00292 
00293   while (Read(inFile)) {
00294     thisLevel = UseLevelOrMakeNew(thisLevel);   // May create new pointer
00295     AddDataToLevel(thisLevel);
00296     nData++;                                    // For debugging purposes
00297   }
00298 
00299   FinishLevel(thisLevel);               // Final  must be completed by hand
00300   
00301   // ---- MGP ---- Don't forget to close the file 
00302   inFile.close();
00303         
00304   //  G4cout << " ==== MakeLevels ===== " << nData << " data read " << G4endl;
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;   // No change
00317 
00318   if (level) FinishLevel(level);        // Save what we have up to now
00319 
00320   //  G4cout << "Making a new level... " << _levelEnergy << G4endl;
00321   return new G4NuclearLevel(_levelEnergy, _halfLife, _angularMomentum);
00322 }
00323 
00324 void G4NuclearLevelManager::AddDataToLevel(G4NuclearLevel* level) 
00325 {
00326   if (!level) return;           // Sanity check
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;               // Sanity check
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 

Generated on Mon May 27 17:49:06 2013 for Geant4 by  doxygen 1.4.7