Geant4  10.02.p01
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 93357 2015-10-19 13:40:13Z 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 #include "G4NucleiProperties.hh"
62 #include "G4PhysicalConstants.hh"
63 
81 
83  const G4String& filename) :
84  _nucleusA(A), _nucleusZ(Z), _validity(false),
85  _levels(0)
86 {
87  if (A <= 0 || Z <= 0 || Z > A ) {
88  throw G4HadronicException(__FILE__, __LINE__,
89  "==== G4NuclearLevelManager ==== (Z,A) <0, or Z>A");
90  }
91  MakeLevels(filename);
92  /*
93  G4cout << "+++++ New G4NuclearLevelManager Z= " << Z << " A= " << A
94  << " " << _levels;
95  if(_levels) { G4cout << " Nlevels= " << _levels->size(); }
96  G4cout << G4endl;
97  */
98 }
99 
101 {
102  ClearLevels();
103 }
104 
106 {
107  if (A <= 0 || Z <= 0 || Z > A ) {
108  throw G4HadronicException(__FILE__, __LINE__,
109  "==== G4NuclearLevelManager ==== (Z,A) <0, or Z>A");
110  }
111  if (_nucleusZ != Z || _nucleusA != A)
112  {
113  _nucleusA = A;
114  _nucleusZ = Z;
115  MakeLevels(filename);
116  }
117 }
118 
120  return (i>=0 && i<NumberOfLevels()) ? (*_levels)[i] : 0;
121 }
122 
123 const G4NuclearLevel*
125 {
126  if (NumberOfLevels() <= 0) return 0;
127 
128  G4int iNear = -1;
129 
130  //G4cout << "G4NuclearLevelManager::NearestLevel E(MeV)= "
131  // << energy/MeV << " dEmax(MeV)= " << eDiffMax/MeV << G4endl;
132 
133  G4double diff = 1.e+10;
134  for (unsigned int i=0; i<_levels->size(); ++i)
135  {
136  G4double e = GetLevel(i)->Energy();
137  G4double eDiff = std::fabs(e - energy);
138  //G4cout << i << ". eDiff(MeV)= " << eDiff/MeV << G4endl;
139  if (eDiff <= diff)
140  {
141  diff = eDiff;
142  iNear = i;
143  }
144  }
145 
146  return GetLevel(iNear); // Includes range checking on iNear
147 }
148 
150 {
151  return (NumberOfLevels() > 0) ? _levels->front()->Energy() : 9999.*GeV;
152 }
153 
155 {
156  return (NumberOfLevels() > 0) ? _levels->back()->Energy() : 0.*GeV;
157 }
158 
160 {
161  return (NumberOfLevels() > 0) ? _levels->front() : 0;
162 }
163 
165 {
166  return (NumberOfLevels() > 0) ? _levels->back() : 0;
167 }
168 
169 G4bool G4NuclearLevelManager::Read(std::ifstream& dataFile)
170 {
171  G4bool goodRead = ReadDataLine(dataFile);
172 
173  if (goodRead) ProcessDataLine();
174  return goodRead;
175 }
176 
177 // NOTE: Standard stream I/O generates a 45 byte std::string per item!
178 
179 G4bool G4NuclearLevelManager::ReadDataLine(std::ifstream& dataFile) {
180  /***** DO NOT USE REGULAR STREAM I/O
181  G4bool result = true;
182  if (dataFile >> _levelEnergy)
183  {
184  dataFile >> _gammaEnergy >> _probability >> _polarity >> _halfLife
185  >> _angularMomentum >> _totalCC >> _kCC >> _l1CC >> _l2CC
186  >> _l3CC >> _m1CC >> _m2CC >> _m3CC >> _m4CC >> _m5CC
187  >> _nPlusCC;
188  }
189  else result = false;
190  *****/
191 
192  // Each item will return iostream status
193  return (ReadDataItem(dataFile, _levelEnergy) &&
194  ReadDataItem(dataFile, _gammaEnergy) &&
195  ReadDataItem(dataFile, _probability) &&
196  ReadDataItem(dataFile, _polarity) &&
197  ReadDataItem(dataFile, _halfLife) &&
198  ReadDataItem(dataFile, _angularMomentum) &&
199  ReadDataItem(dataFile, _totalCC) &&
200  ReadDataItem(dataFile, _kCC) &&
201  ReadDataItem(dataFile, _l1CC) &&
202  ReadDataItem(dataFile, _l2CC) &&
203  ReadDataItem(dataFile, _l3CC) &&
204  ReadDataItem(dataFile, _m1CC) &&
205  ReadDataItem(dataFile, _m2CC) &&
206  ReadDataItem(dataFile, _m3CC) &&
207  ReadDataItem(dataFile, _m4CC) &&
208  ReadDataItem(dataFile, _m5CC) &&
209  ReadDataItem(dataFile, _nPlusCC) );
210 }
211 
212 G4bool
214 {
215  // G4bool okay = (dataFile >> buffer) != 0; // Get next token
216  // if (okay) x = strtod(buffer, NULL);
217  char buffer[30];
218  for(G4int i=0; i<30; ++i) { buffer[i] = 0; }
219  G4bool okay = true;
220  dataFile >> buffer;
221  if(dataFile.fail()) { okay = false; }
222  else { x = strtod(buffer, NULL); }
223 
224  return okay;
225 }
226 
228 {
229  const G4double minProbability = 1e-8;
230 
231  // Assign units for dimensional quantities
232  _levelEnergy *= keV;
233  _gammaEnergy *= keV;
234  _halfLife *= second;
235 
236  // The following adjustment is needed to take care of anomalies in
237  // data files, where some transitions show up with relative probability
238  // zero
239  if (_probability < minProbability) _probability = minProbability;
240  // the folowwing is to convert icc probability to accumulative ones
241  _l1CC += _kCC;
242  _l2CC += _l1CC;
243  _l3CC += _l2CC;
244  _m1CC += _l3CC;
245  _m2CC += _m1CC;
246  _m3CC += _m2CC;
247  _m4CC += _m3CC;
248  _m5CC += _m4CC;
249  _nPlusCC += _m5CC;
250 
251  if (_nPlusCC!=0) { // Normalize to probabilities
252  _kCC /= _nPlusCC;
253  _l1CC /= _nPlusCC;
254  _l2CC /= _nPlusCC;
255  _l3CC /= _nPlusCC;
256  _m1CC /= _nPlusCC;
257  _m2CC /= _nPlusCC;
258  _m3CC /= _nPlusCC;
259  _m4CC /= _nPlusCC;
260  _m5CC /= _nPlusCC;
261  _nPlusCC /= _nPlusCC;
262  } else { // Total was zero, reset to unity
263  _kCC = 1;
264  _l1CC = 1;
265  _l2CC = 1;
266  _l3CC = 1;
267  _m1CC = 1;
268  _m2CC = 1;
269  _m3CC = 1;
270  _m4CC = 1;
271  _m5CC = 1;
272  _nPlusCC = 1;
273  }
274 
275  // G4cout << "Read " << _levelEnergy << " " << _gammaEnergy << " " << _probability << G4endl;
276 }
277 
279 {
280  _validity = false;
281 
282  if (NumberOfLevels() > 0) {
283  std::for_each(_levels->begin(), _levels->end(), DeleteLevel());
284  _levels->clear();
285  }
286 
287  delete _levels;
288  _levels = 0;
289 }
290 
292 {
293  _validity = false;
294  if (NumberOfLevels() > 0) ClearLevels(); // Discard existing data
295 
296  std::ifstream inFile(filename, std::ios::in);
297  if (! inFile)
298  {
299 #ifdef GAMMAFILEWARNING
300  if (_nucleusZ > 10) G4cout << " G4NuclearLevelManager: nuclide ("
301  << _nucleusZ << "," << _nucleusA
302  << ") does not have a gamma levels file" << G4endl;
303 #endif
304  return;
305  }
306 
308 
309  // Read individual gamma data and fill levels for this nucleus
310 
311  G4NuclearLevel* thisLevel = 0;
312  G4int nData = 0;
313 
314  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
315  while (Read(inFile)) {
316  thisLevel = UseLevelOrMakeNew(thisLevel); // May create new pointer
317  AddDataToLevel(thisLevel);
318  nData++; // For debugging purposes
319  }
320 
321  FinishLevel(thisLevel); // Final must be completed by hand
322 
323  // ---- MGP ---- Don't forget to close the file
324  inFile.close();
325 
326  // G4cout << " ==== MakeLevels ===== " << nData << " data read " << G4endl;
327 
328  G4PtrSort<G4NuclearLevel>(_levels);
329 
330  _validity = true;
331 
332  return;
333 }
334 
337 {
338  if (level && _levelEnergy == level->Energy()) return level; // No change
339 
340  if (level) FinishLevel(level); // Save what we have up to now
341 
342  // G4cout << "Making a new level... " << _levelEnergy << G4endl;
344 }
345 
347 {
348  if (!level) return; // Sanity check
349 
350  level->_energies.push_back(_gammaEnergy);
351  level->_weights.push_back(_probability);
352  level->_polarities.push_back(_polarity);
353  level->_kCC.push_back(_kCC);
354  level->_l1CC.push_back(_l1CC);
355  level->_l2CC.push_back(_l2CC);
356  level->_l3CC.push_back(_l3CC);
357  level->_m1CC.push_back(_m1CC);
358  level->_m2CC.push_back(_m2CC);
359  level->_m3CC.push_back(_m3CC);
360  level->_m4CC.push_back(_m4CC);
361  level->_m5CC.push_back(_m5CC);
362  level->_nPlusCC.push_back(_nPlusCC);
363  level->_totalCC.push_back(_totalCC);
364 }
365 
367 {
368  if (!level || !_levels) return; // Sanity check
369 
370  level->Finalize();
371  _levels->push_back(level);
372 }
373 
374 
376 {
377  G4int nLevels = NumberOfLevels();
378 
379  G4cout << " ==== G4NuclearLevelManager ==== (" << _nucleusZ << ", " << _nucleusA
380  << ") has " << nLevels << " levels" << G4endl
381  << "Highest level is at energy " << MaxLevelEnergy() << " MeV "
382  << G4endl << "Lowest level is at energy " << MinLevelEnergy()
383  << " MeV " << G4endl;
384 
385  for (G4int i=0; i<nLevels; ++i) {
386  GetLevel(i)->PrintAll();
387  }
388 }
389 
391 {
392  G4int nLevels = NumberOfLevels();
394  + neutron_mass_c2 -
396 
397  G4cout << "Z= " << _nucleusZ << " A= " << _nucleusA
398  << " " << nLevels << " levels"
399  << " Efermi(MeV)= " << efermi << G4endl;
400 
401  for (G4int i=0; i<nLevels; ++i) {
402  GetLevel(i)->PrintLevels();
403  }
404 }
405 
407 {
408  _nucleusA = right._nucleusA;
409  _nucleusZ = right._nucleusZ;
410  _validity = right._validity;
411 
412  if (right._levels != 0)
413  {
415  G4int n = right._levels->size();
416  G4int i;
417  for (i=0; i<n; ++i)
418  {
419  _levels->push_back(new G4NuclearLevel(*(right.GetLevel(i))));
420  }
421  G4PtrSort<G4NuclearLevel>(_levels);
422  }
423  else
424  {
425  _levels = 0;
426  }
427 }
428 
std::vector< G4double > _m3CC
static G4double GetNuclearMass(const G4double A, const G4double Z)
std::vector< G4double > _m1CC
G4NuclearLevelManager(G4int Z, G4int A, const G4String &filename)
G4double Energy() const
const G4NuclearLevel * HighestLevel() const
G4double MinLevelEnergy() const
std::vector< G4double > _l3CC
#define buffer
Definition: xmlparse.cc:628
G4NuclearLevel * UseLevelOrMakeNew(G4NuclearLevel *level)
std::vector< G4double > _weights
std::vector< G4double > _polarities
G4bool Read(std::ifstream &aDataFile)
int G4int
Definition: G4Types.hh:78
G4bool ReadDataItem(std::istream &dataFile, G4double &x)
std::vector< G4NuclearLevel * > G4PtrLevelVector
std::vector< G4double > _m5CC
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
std::vector< G4double > _l1CC
void PrintAll() const
void SetNucleus(G4int Z, G4int A, const G4String &filename)
const G4NuclearLevel * LowestLevel() const
const G4NuclearLevel * GetLevel(G4int i) const
bool G4bool
Definition: G4Types.hh:79
std::vector< G4double > _nPlusCC
std::vector< G4double > _energies
void MakeLevels(const G4String &filename)
std::vector< G4double > _m2CC
static const double second
Definition: G4SIunits.hh:156
static const double GeV
Definition: G4SIunits.hh:214
const G4int n
G4bool ReadDataLine(std::ifstream &dataFile)
std::vector< G4double > _kCC
G4double energy(const ThreeVector &p, const G4double m)
const G4double x[NPOINTSGL]
void FinishLevel(G4NuclearLevel *level)
#define G4endl
Definition: G4ios.hh:61
static const double keV
Definition: G4SIunits.hh:213
std::vector< G4double > _m4CC
double G4double
Definition: G4Types.hh:76
void AddDataToLevel(G4NuclearLevel *level)
std::vector< G4double > _l2CC
G4double MaxLevelEnergy() const
const G4NuclearLevel * NearestLevel(G4double energy, G4double eDiffMax=1.e+8) const
std::vector< G4double > _totalCC
void PrintLevels() const