Geant4_10
G4NuclearLevelManager.cc
Go to the documentation of this file.
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 77025 2013-11-20 16:11:51Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 // GEANT 4 class file
30 //
31 // CERN, Geneva, Switzerland
32 //
33 // File name: G4NuclearLevelManager
34 //
35 // Author: Maria Grazia Pia (pia@genova.infn.it)
36 //
37 // Creation date: 24 October 1998
38 //
39 // Modifications:
40 //
41 // 15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
42 // Added half-life, angular momentum, parity, emissioni type
43 // reading from experimental data.
44 // 02 May 2003, Vladimir Ivanchenko remove rublic copy constructor
45 // 06 Oct 2010, M. Kelsey -- Use object storage, not pointers, drop
46 // public access to list, simplify list construction
47 // -------------------------------------------------------------------
48 #include <stdlib.h>
49 #include <fstream>
50 #include <sstream>
51 #include <algorithm>
52 
53 #include "G4NuclearLevelManager.hh"
54 
55 #include "globals.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "G4NuclearLevel.hh"
58 #include "G4ios.hh"
59 #include "G4HadronicException.hh"
60 #include "G4HadTmpUtil.hh"
61 
63  const G4String& filename) :
64  _nucleusA(A), _nucleusZ(Z), _fileName(filename), _validity(false),
65  _levels(0), _levelEnergy(0), _gammaEnergy(0), _probability(0)
66 {
67  if (A <= 0 || Z <= 0 || Z > A ) {
68  throw G4HadronicException(__FILE__, __LINE__,
69  "==== G4NuclearLevelManager ==== (Z,A) <0, or Z>A");
70  }
71  for(G4int i=0; i<30; ++i) { buffer[i] = 0; }
72  MakeLevels();
73 }
74 
76 {
77  ClearLevels();
78 }
79 
81 {
82  if (A <= 0 || Z <= 0 || Z > A ) {
83  throw G4HadronicException(__FILE__, __LINE__,
84  "==== G4NuclearLevelManager ==== (Z,A) <0, or Z>A");
85  }
86  if (_nucleusZ != Z || _nucleusA != A)
87  {
88  _nucleusA = A;
89  _nucleusZ = Z;
90  _fileName = filename;
91  MakeLevels();
92  }
93 }
94 
96  return (i>=0 && i<NumberOfLevels()) ? (*_levels)[i] : 0;
97 }
98 
99 const G4NuclearLevel*
101  G4double eDiffMax) const
102 {
103  if (NumberOfLevels() <= 0) return 0;
104 
105  G4int iNear = -1;
106 
107  //G4cout << "G4NuclearLevelManager::NearestLevel E(MeV)= "
108  // << energy/MeV << " dEmax(MeV)= " << eDiffMax/MeV << G4endl;
109 
110  G4double diff = 9999. * GeV;
111  for (unsigned int i=0; i<_levels->size(); i++)
112  {
113  G4double e = GetLevel(i)->Energy();
114  G4double eDiff = std::abs(e - energy);
115  //G4cout << i << ". eDiff(MeV)= " << eDiff/MeV << G4endl;
116  if (eDiff < diff && eDiff <= eDiffMax)
117  {
118  diff = eDiff;
119  iNear = i;
120  }
121  }
122 
123  return GetLevel(iNear); // Includes range checking on iNear
124 }
125 
127 {
128  return (NumberOfLevels() > 0) ? _levels->front()->Energy() : 9999.*GeV;
129 }
130 
132 {
133  return (NumberOfLevels() > 0) ? _levels->back()->Energy() : 0.*GeV;
134 }
135 
137 {
138  return (NumberOfLevels() > 0) ? _levels->front() : 0;
139 }
140 
142 {
143  return (NumberOfLevels() > 0) ? _levels->back() : 0;
144 }
145 
146 G4bool G4NuclearLevelManager::Read(std::ifstream& dataFile)
147 {
148  G4bool goodRead = ReadDataLine(dataFile);
149 
150  if (goodRead) ProcessDataLine();
151  return goodRead;
152 }
153 
154 // NOTE: Standard stream I/O generates a 45 byte std::string per item!
155 
156 G4bool G4NuclearLevelManager::ReadDataLine(std::ifstream& dataFile) {
157  /***** DO NOT USE REGULAR STREAM I/O
158  G4bool result = true;
159  if (dataFile >> _levelEnergy)
160  {
161  dataFile >> _gammaEnergy >> _probability >> _polarity >> _halfLife
162  >> _angularMomentum >> _totalCC >> _kCC >> _l1CC >> _l2CC
163  >> _l3CC >> _m1CC >> _m2CC >> _m3CC >> _m4CC >> _m5CC
164  >> _nPlusCC;
165  }
166  else result = false;
167  *****/
168 
169  // Each item will return iostream status
170  return (ReadDataItem(dataFile, _levelEnergy) &&
171  ReadDataItem(dataFile, _gammaEnergy) &&
172  ReadDataItem(dataFile, _probability) &&
173  ReadDataItem(dataFile, _polarity) &&
174  ReadDataItem(dataFile, _halfLife) &&
175  ReadDataItem(dataFile, _angularMomentum) &&
176  ReadDataItem(dataFile, _totalCC) &&
177  ReadDataItem(dataFile, _kCC) &&
178  ReadDataItem(dataFile, _l1CC) &&
179  ReadDataItem(dataFile, _l2CC) &&
180  ReadDataItem(dataFile, _l3CC) &&
181  ReadDataItem(dataFile, _m1CC) &&
182  ReadDataItem(dataFile, _m2CC) &&
183  ReadDataItem(dataFile, _m3CC) &&
184  ReadDataItem(dataFile, _m4CC) &&
185  ReadDataItem(dataFile, _m5CC) &&
186  ReadDataItem(dataFile, _nPlusCC) );
187 }
188 
189 G4bool
190 G4NuclearLevelManager::ReadDataItem(std::istream& dataFile, G4double& x)
191 {
192  // G4bool okay = (dataFile >> buffer) != 0; // Get next token
193  // if (okay) x = strtod(buffer, NULL);
194  G4bool okay = true;
195  dataFile >> buffer;
196  if(dataFile.fail()) { okay = false; }
197  else { x = strtod(buffer, NULL); }
198 
199  return okay;
200 }
201 
202 void G4NuclearLevelManager::ProcessDataLine()
203 {
204  const G4double minProbability = 1e-8;
205 
206  // Assign units for dimensional quantities
207  _levelEnergy *= keV;
208  _gammaEnergy *= keV;
209  _halfLife *= second;
210 
211  // The following adjustment is needed to take care of anomalies in
212  // data files, where some transitions show up with relative probability
213  // zero
214  if (_probability < minProbability) _probability = minProbability;
215  // the folowwing is to convert icc probability to accumulative ones
216  _l1CC += _kCC;
217  _l2CC += _l1CC;
218  _l3CC += _l2CC;
219  _m1CC += _l3CC;
220  _m2CC += _m1CC;
221  _m3CC += _m2CC;
222  _m4CC += _m3CC;
223  _m5CC += _m4CC;
224  _nPlusCC += _m5CC;
225 
226  if (_nPlusCC!=0) { // Normalize to probabilities
227  _kCC /= _nPlusCC;
228  _l1CC /= _nPlusCC;
229  _l2CC /= _nPlusCC;
230  _l3CC /= _nPlusCC;
231  _m1CC /= _nPlusCC;
232  _m2CC /= _nPlusCC;
233  _m3CC /= _nPlusCC;
234  _m4CC /= _nPlusCC;
235  _m5CC /= _nPlusCC;
236  _nPlusCC /= _nPlusCC;
237  } else { // Total was zero, reset to unity
238  _kCC = 1;
239  _l1CC = 1;
240  _l2CC = 1;
241  _l3CC = 1;
242  _m1CC = 1;
243  _m2CC = 1;
244  _m3CC = 1;
245  _m4CC = 1;
246  _m5CC = 1;
247  _nPlusCC = 1;
248  }
249 
250  // G4cout << "Read " << _levelEnergy << " " << _gammaEnergy << " " << _probability << G4endl;
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 
266 void 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 
310 G4NuclearLevelManager::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 
320 void 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 
340 void G4NuclearLevelManager::FinishLevel(G4NuclearLevel* level)
341 {
342  if (!level || !_levels) return; // Sanity check
343 
344  level->Finalize();
345  _levels->push_back(level);
346 }
347 
348 
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 
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 
ifstream in
Definition: comparison.C:7
G4NuclearLevelManager(G4int Z, G4int A, const G4String &filename)
G4double Energy() const
const G4NuclearLevel * HighestLevel() const
G4double MinLevelEnergy() const
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
std::vector< G4NuclearLevel * > G4PtrLevelVector
Char_t n[5]
G4GLOB_DLL std::ostream G4cout
double energy
Definition: plottest35.C:25
void PrintAll() const
void SetNucleus(G4int Z, G4int A, const G4String &filename)
Float_t Z
Definition: plot.C:39
const G4NuclearLevel * LowestLevel() const
const G4NuclearLevel * GetLevel(G4int i) const
bool G4bool
Definition: G4Types.hh:79
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
const G4NuclearLevel * NearestLevel(G4double energy, G4double eDiffMax=9999.*CLHEP::GeV) const
G4double MaxLevelEnergy() const