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

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

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