source: trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4NuclearLevelManager.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 11.5 KB
Line 
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// $Id: G4NuclearLevelManager.cc,v 1.14 2010/11/17 16:50:53 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
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
46//        06 Oct 2010, M. Kelsey -- Use object storage, not pointers, drop
47//              public access to list, simplify list construction
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
62G4NuclearLevelManager::G4NuclearLevelManager():
63    _nucleusA(0), _nucleusZ(0), _fileName(""), _validity(false), 
64    _levels(0), _levelEnergy(0), _gammaEnergy(0), _probability(0)
65{ }
66
67G4NuclearLevelManager::G4NuclearLevelManager(G4int Z, G4int A, const G4String& filename) :
68    _nucleusA(A), _nucleusZ(Z), _fileName(filename), _validity(false), 
69    _levels(0), _levelEnergy(0), _gammaEnergy(0), _probability(0)
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
77G4NuclearLevelManager::~G4NuclearLevelManager()
78{
79  ClearLevels();
80}
81
82void G4NuclearLevelManager::SetNucleus(G4int Z, G4int A, const G4String& filename)
83{
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)
88    {
89      _nucleusA = A;
90      _nucleusZ = Z;
91      _fileName = filename;
92      MakeLevels();
93    }
94}
95
96const G4NuclearLevel* G4NuclearLevelManager::GetLevel(G4int i) const {
97  return (i>=0 && i<NumberOfLevels()) ? (*_levels)[i] : 0;
98}
99
100
101const G4NuclearLevel* 
102G4NuclearLevelManager::NearestLevel(G4double energy, 
103                                    G4double eDiffMax) const 
104{
105  if (NumberOfLevels() <= 0) return 0;
106
107  G4int iNear = -1;
108 
109  G4double diff = 9999. * GeV;
110  for (unsigned int i=0; i<_levels->size(); i++)
111    {
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;
118        }
119    }
120 
121  return GetLevel(iNear);       // Includes range checking on iNear
122}
123
124
125G4double G4NuclearLevelManager::MinLevelEnergy() const
126{
127  return (NumberOfLevels() > 0) ? _levels->front()->Energy() : 9999.*GeV;
128}
129
130
131G4double G4NuclearLevelManager::MaxLevelEnergy() const
132{
133  return (NumberOfLevels() > 0) ? _levels->back()->Energy() : 0.*GeV;
134}
135
136
137const G4NuclearLevel* G4NuclearLevelManager::HighestLevel() const
138{
139  return (NumberOfLevels() > 0) ? _levels->front() : 0;
140}
141
142
143const G4NuclearLevel* G4NuclearLevelManager::LowestLevel() const
144{
145  return (NumberOfLevels() > 0) ? _levels->back() : 0;
146}
147
148
149G4bool G4NuclearLevelManager::Read(std::ifstream& dataFile) 
150{
151  G4bool goodRead = ReadDataLine(dataFile);
152 
153  if (goodRead) ProcessDataLine();
154  return goodRead;
155}
156
157// NOTE:  Standard stream I/O generates a 45 byte std::string per item!
158
159G4bool G4NuclearLevelManager::ReadDataLine(std::ifstream& dataFile) {
160  /***** DO NOT USE REGULAR STREAM I/O
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    }
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) );
190}
191
192G4bool
193G4NuclearLevelManager::ReadDataItem(std::istream& dataFile, G4double& x) 
194{
195  G4bool okay = (dataFile >> buffer);           // Get next token
196  if (okay) x = strtod(buffer, NULL);
197
198  return okay;
199}
200
201void G4NuclearLevelManager::ProcessDataLine() 
202{
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
253void 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
266void G4NuclearLevelManager::MakeLevels()
267{
268  _validity = false;
269  if (NumberOfLevels() > 0) ClearLevels();      // Discard existing data
270
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;
283
284  // Read individual gamma data and fill levels for this nucleus
285 
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
296 
297  // ---- MGP ---- Don't forget to close the file
298  inFile.close();
299       
300  //  G4cout << " ==== MakeLevels ===== " << nData << " data read " << G4endl;
301
302  G4PtrSort<G4NuclearLevel>(_levels);
303 
304  _validity = true;
305 
306  return;
307}
308
309G4NuclearLevel* 
310G4NuclearLevelManager::UseLevelOrMakeNew(G4NuclearLevel* level) 
311{
312  if (level && _levelEnergy == level->Energy()) return level;   // No change
313
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);
318}
319
320void G4NuclearLevelManager::AddDataToLevel(G4NuclearLevel* level) 
321{
322  if (!level) return;           // Sanity check
323
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
340void G4NuclearLevelManager::FinishLevel(G4NuclearLevel* level) 
341{
342  if (!level || !_levels) return;               // Sanity check
343
344  level->Finalize();
345  _levels->push_back(level);
346}
347
348
349void G4NuclearLevelManager::PrintAll()
350{
351  G4int nLevels = NumberOfLevels();
352   
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;
358   
359  for (G4int i=0; i<nLevels; i++)
360    GetLevel(i)->PrintAll();
361}
362
363G4NuclearLevelManager::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;
386
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          {
394            _levels->push_back(new G4NuclearLevel(*(right.GetLevel(i))));
395          }
396       
397        G4PtrSort<G4NuclearLevel>(_levels);
398      }
399    else 
400      {
401        _levels = 0;
402      }
403}
404
Note: See TracBrowser for help on using the repository browser.