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

Last change on this file since 1245 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 13.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//
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//     
45// -------------------------------------------------------------------
46
47#include "G4NuclearLevelManager.hh"
48
49#include "globals.hh"
50#include "G4NuclearLevel.hh"
51#include "G4ios.hh"
52#include "G4HadronicException.hh"
53#include <stdlib.h>
54#include <fstream>
55#include <sstream>
56#include <algorithm>
57#include "G4HadTmpUtil.hh"
58
59G4NuclearLevelManager::G4NuclearLevelManager():
60    _nucleusA(0), _nucleusZ(0), _fileName(""), _validity(false), 
61    _levels(0), _levelEnergy(0), _gammaEnergy(0), _probability(0)
62{ }
63
64G4NuclearLevelManager::G4NuclearLevelManager(const G4int Z, const G4int A, const G4String& filename) :
65    _nucleusA(A), _nucleusZ(Z), _fileName(filename)
66{ 
67    if (A <= 0 || Z <= 0 || Z > A )
68        throw G4HadronicException(__FILE__, __LINE__, "==== G4NuclearLevelManager ==== (Z,A) <0, or Z>A");
69
70    _levels = 0;
71
72    MakeLevels();
73}
74
75G4NuclearLevelManager::~G4NuclearLevelManager()
76{ 
77  if ( _levels ) {
78    if (_levels->size()>0) 
79      {
80        std::for_each(_levels->begin(), _levels->end(), DeleteLevel());
81        _levels->clear();
82      }
83    delete _levels;
84    _levels = 0;
85  }
86}
87
88void G4NuclearLevelManager::SetNucleus(const G4int Z, const G4int A, const G4String& filename)
89{
90    if (_nucleusZ != Z || _nucleusA != A)
91    {
92        _nucleusA = A;
93        _nucleusZ = Z;
94        _fileName = filename;
95        MakeLevels();
96    }
97   
98}
99
100G4bool G4NuclearLevelManager::IsValid() const
101{
102    return _validity;
103}
104
105
106G4int G4NuclearLevelManager::NumberOfLevels() const
107{
108    G4int n = 0;
109    if (_levels != 0) n = _levels->size();
110    return n;
111}
112
113
114const G4PtrLevelVector* G4NuclearLevelManager::GetLevels() const
115{
116    return _levels;
117}
118
119
120const G4NuclearLevel* G4NuclearLevelManager::
121NearestLevel(const G4double energy, const G4double eDiffMax) const
122{
123    G4int iNear = -1;
124 
125    G4double diff = 9999. * GeV;
126    if (_levels != 0)
127    {
128        unsigned int i = 0;
129        for (i=0; i<_levels->size(); i++)
130        {
131            G4double e = _levels->operator[](i)->Energy();
132            G4double eDiff = std::abs(e - energy);
133            if (eDiff < diff && eDiff <= eDiffMax)
134            { 
135                diff = eDiff; 
136                iNear = i;
137            }
138        }
139    }
140    if (_levels != 0 && iNear >= 0 && iNear < static_cast<G4int>(_levels->size()) )
141    { 
142        return _levels->operator[](iNear); 
143    }
144    else
145    { 
146        return 0; 
147    }
148}
149
150
151G4double G4NuclearLevelManager::MinLevelEnergy() const
152{
153    G4double eMin = 9999.*GeV;
154    if (_levels != 0)
155    {
156        if (_levels->size() > 0) eMin = _levels->front()->Energy(); 
157    }
158    return eMin;
159}
160
161
162G4double G4NuclearLevelManager::MaxLevelEnergy() const
163{
164    G4double eMax = 0.;
165    if (_levels != 0)
166    {
167        if (_levels->size() > 0) eMax = _levels->back()->Energy(); 
168    }
169    return eMax;
170}
171
172
173const G4NuclearLevel* G4NuclearLevelManager::HighestLevel() const
174{
175    if (_levels!= 0 && _levels->size() > 0) return _levels->front(); 
176    else return 0; 
177}
178
179
180const G4NuclearLevel* G4NuclearLevelManager::LowestLevel() const
181{
182    if (_levels != 0 && _levels->size() > 0) return _levels->back();
183    else return 0;
184}
185
186
187G4bool G4NuclearLevelManager::Read(std::ifstream& dataFile)
188{
189  const G4double minProbability = 1e-8;
190 
191  G4bool result = true;
192
193  if (dataFile >> _levelEnergy)
194    {
195      dataFile >> _gammaEnergy >> _probability >> _polarity >> _halfLife
196               >> _angularMomentum  >> _totalCC >> _kCC >> _l1CC >> _l2CC
197               >> _l3CC >> _m1CC >> _m2CC >> _m3CC >> _m4CC >> _m5CC
198               >> _nPlusCC;
199      _levelEnergy *= keV;
200      _gammaEnergy *= keV;
201      _halfLife *= second;
202       
203      // The following adjustment is needed to take care of anomalies in
204      // data files, where some transitions show up with relative probability
205      // zero
206      if (_probability < minProbability) _probability = minProbability;
207      // the folowwing is to convert icc probability to accumulative ones
208      _l1CC += _kCC;
209      _l2CC += _l1CC;
210      _l3CC += _l2CC;
211      _m1CC += _l3CC;
212      _m2CC += _m1CC;
213      _m3CC += _m2CC;
214      _m4CC += _m3CC;
215      _m5CC += _m4CC;
216      _nPlusCC += _m5CC;
217      if(_nPlusCC!=0)
218      {
219        _kCC /= _nPlusCC;
220        _l1CC /= _nPlusCC;
221        _l2CC /= _nPlusCC;
222        _l3CC /= _nPlusCC;
223        _m1CC /= _nPlusCC;
224        _m2CC /= _nPlusCC;
225        _m3CC /= _nPlusCC;
226        _m4CC /= _nPlusCC;
227        _m5CC /= _nPlusCC;
228        _nPlusCC /= _nPlusCC; 
229      }
230      else
231      {
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    else
247    {
248        result = false;
249    }
250   
251    return result;
252}
253
254
255void G4NuclearLevelManager::MakeLevels()
256{
257  _validity = false;
258  std::ifstream inFile(_fileName, std::ios::in);
259  if (! inFile) 
260    {
261#ifdef GAMMAFILEWARNING
262      if (_nucleusZ > 10) G4cout << " G4NuclearLevelManager: nuclide (" 
263                                 << _nucleusZ << "," << _nucleusA
264                                 << ") does not have a gamma levels file" << G4endl;
265#endif
266      return;
267    }
268 
269  if (_levels != 0)
270    {
271      if (_levels->size()>0) 
272        {
273          std::vector<G4NuclearLevel*>::iterator pos;
274          for(pos=_levels->begin(); pos!=_levels->end(); pos++)
275            if (*pos) delete *pos;
276          _levels->clear();
277        }
278      delete _levels;
279    }
280  else 
281    {
282      _validity = true;
283    }
284
285  _levels = new G4PtrLevelVector;
286 
287  std::vector<G4double> eLevel;
288  std::vector<G4double> eGamma;
289  std::vector<G4double> wGamma;
290  std::vector<G4double> pGamma; // polarity
291  std::vector<G4double> hLevel; // half life
292  std::vector<G4double> aLevel; // angular momentum
293  std::vector<G4double> kConve; //  internal convertion coefficiencies
294  std::vector<G4double> l1Conve;
295  std::vector<G4double> l2Conve;
296  std::vector<G4double> l3Conve;
297  std::vector<G4double> m1Conve;
298  std::vector<G4double> m2Conve;
299  std::vector<G4double> m3Conve;
300  std::vector<G4double> m4Conve;
301  std::vector<G4double> m5Conve;
302  std::vector<G4double> npConve;
303  std::vector<G4double> toConve;
304 
305       
306  while (Read(inFile))
307    {
308      eLevel.push_back(_levelEnergy);
309      eGamma.push_back(_gammaEnergy);
310      wGamma.push_back(_probability);
311      pGamma.push_back(_polarity);
312      hLevel.push_back(_halfLife);
313      aLevel.push_back(_angularMomentum);
314      kConve.push_back(_kCC);
315      l1Conve.push_back(_l1CC);
316      l2Conve.push_back(_l2CC);
317      l3Conve.push_back(_l3CC);
318      m1Conve.push_back(_m1CC);
319      m2Conve.push_back(_m2CC);
320      m3Conve.push_back(_m3CC);
321      m4Conve.push_back(_m4CC);
322      m5Conve.push_back(_m5CC);
323      npConve.push_back(_nPlusCC);
324      toConve.push_back(_totalCC);
325    }
326 
327  // ---- MGP ---- Don't forget to close the file
328  inFile.close();
329       
330  G4int nData = eLevel.size();
331       
332  //  G4cout << " ==== MakeLevels ===== " << nData << " data read " << G4endl;
333       
334  G4double thisLevelEnergy = eLevel[0];
335  G4double thisLevelHalfLife = 0.;
336  G4double thisLevelAngMom = 0.;
337  std::vector<G4double> thisLevelEnergies;
338  std::vector<G4double> thisLevelWeights;
339  std::vector<G4double> thisLevelPolarities;
340  std::vector<G4double> thisLevelkCC;
341  std::vector<G4double> thisLevell1CC;
342  std::vector<G4double> thisLevell2CC;
343  std::vector<G4double> thisLevell3CC;
344  std::vector<G4double> thisLevelm1CC;
345  std::vector<G4double> thisLevelm2CC;
346  std::vector<G4double> thisLevelm3CC;
347  std::vector<G4double> thisLevelm4CC;
348  std::vector<G4double> thisLevelm5CC;
349  std::vector<G4double> thisLevelnpCC;
350  std::vector<G4double> thisLeveltoCC;
351 
352  G4double e = -1.;
353  G4int i;
354  for (i=0; i<nData; i++)
355    {
356      e = eLevel[i];
357      if (e != thisLevelEnergy)
358        {
359          //      G4cout << "Making a new level... " << e << " "
360          //             << thisLevelEnergies.entries() << " "
361          //             << thisLevelWeights.entries() << G4endl;
362         
363          G4NuclearLevel* newLevel = new G4NuclearLevel(thisLevelEnergy,
364                                                        thisLevelHalfLife,
365                                                        thisLevelAngMom,
366                                                        thisLevelEnergies,
367                                                        thisLevelWeights,
368                                                        thisLevelPolarities,
369                                                        thisLevelkCC,
370                                                        thisLevell1CC,
371                                                        thisLevell2CC,
372                                                        thisLevell3CC,
373                                                        thisLevelm1CC,
374                                                        thisLevelm2CC,
375                                                        thisLevelm3CC,
376                                                        thisLevelm4CC,
377                                                        thisLevelm5CC,
378                                                        thisLevelnpCC,
379                                                        thisLeveltoCC );
380                                                       
381          _levels->push_back(newLevel);
382          // Reset data vectors
383          thisLevelEnergies.clear();
384          thisLevelWeights.clear();
385          thisLevelPolarities.clear();
386          thisLevelkCC.clear();
387          thisLevell1CC.clear();
388          thisLevell2CC.clear();
389          thisLevell3CC.clear();
390          thisLevelm1CC.clear();
391          thisLevelm2CC.clear();
392          thisLevelm3CC.clear();
393          thisLevelm4CC.clear();
394          thisLevelm5CC.clear();
395          thisLevelnpCC.clear();
396          thisLeveltoCC.clear();
397          thisLevelEnergy = e;
398        }
399      // Append current data
400      thisLevelEnergies.push_back(eGamma[i]);
401      thisLevelWeights.push_back(wGamma[i]);
402      thisLevelPolarities.push_back(pGamma[i]);
403      thisLevelkCC.push_back(kConve[i]);
404      thisLevell1CC.push_back(l1Conve[i]);
405      thisLevell2CC.push_back(l2Conve[i]);
406      thisLevell3CC.push_back(l3Conve[i]);
407      thisLevelm1CC.push_back(m1Conve[i]);
408      thisLevelm2CC.push_back(m2Conve[i]);
409      thisLevelm3CC.push_back(m3Conve[i]);
410      thisLevelm4CC.push_back(m4Conve[i]);
411      thisLevelm5CC.push_back(m5Conve[i]);
412      thisLevelnpCC.push_back(npConve[i]);
413      thisLeveltoCC.push_back(toConve[i]);
414      thisLevelHalfLife = hLevel[i];
415      thisLevelAngMom = aLevel[i];
416    }
417    // Make last level
418    if (e > 0.)
419      {
420        G4NuclearLevel* newLevel = new G4NuclearLevel(e,thisLevelHalfLife,
421                                                      thisLevelAngMom,
422                                                      thisLevelEnergies,
423                                                      thisLevelWeights,
424                                                      thisLevelPolarities,
425                                                      thisLevelkCC,
426                                                      thisLevell1CC,
427                                                      thisLevell2CC,
428                                                      thisLevell3CC,
429                                                      thisLevelm1CC,
430                                                      thisLevelm2CC,
431                                                      thisLevelm3CC,
432                                                      thisLevelm4CC,
433                                                      thisLevelm5CC,
434                                                      thisLevelnpCC,
435                                                      thisLeveltoCC );
436
437        _levels->push_back(newLevel);
438    }
439
440        G4PtrSort<G4NuclearLevel>(_levels);
441
442    return;
443}
444
445
446void G4NuclearLevelManager::PrintAll()
447{
448    G4int nLevels = 0;
449    if (_levels != 0) nLevels = _levels->size();
450   
451    G4cout << " ==== G4NuclearLevelManager ==== (" << _nucleusZ << ", " << _nucleusA
452           << ") has " << nLevels << " levels" << G4endl
453           << "Highest level is at energy " << MaxLevelEnergy() << " MeV "
454           << G4endl << "Lowest level is at energy " << MinLevelEnergy()
455           << " MeV " << G4endl;
456   
457    G4int i = 0;
458    for (i=0; i<nLevels; i++)
459    { _levels->operator[](i)->PrintAll(); }
460}
461
462
463G4NuclearLevelManager::G4NuclearLevelManager(const G4NuclearLevelManager &right)
464{
465    _levelEnergy = right._levelEnergy;
466    _gammaEnergy = right._gammaEnergy;
467    _probability = right._probability;
468    _polarity = right._polarity;
469    _halfLife = right._halfLife;
470    _angularMomentum = right._angularMomentum;
471    _kCC = right._kCC;
472    _l1CC = right._l1CC;
473    _l2CC = right._l2CC;
474    _l3CC = right._l3CC;
475    _m1CC = right._m1CC;
476    _m2CC = right._m2CC;
477    _m3CC = right._m3CC;
478    _m4CC = right._m4CC;
479    _m5CC = right._m5CC;
480    _nPlusCC = right._nPlusCC;
481    _totalCC = right._totalCC;
482    _nucleusA = right._nucleusA;
483    _nucleusZ = right._nucleusZ;
484    _fileName = right._fileName;
485    _validity = right._validity;
486    if (right._levels != 0)   
487      {
488        _levels = new G4PtrLevelVector;
489        G4int n = right._levels->size();
490        G4int i;
491        for (i=0; i<n; i++)
492          {
493            _levels->push_back(new G4NuclearLevel(*(right._levels->operator[](i))));
494          }
495       
496        G4PtrSort<G4NuclearLevel>(_levels);
497      }
498    else 
499      {
500        _levels = 0;
501      }
502}
503
Note: See TracBrowser for help on using the repository browser.