Geant4_10
G4NuclideTable.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 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
27 //
28 // MODULE: G4NuclideTable.cc
29 //
30 // Date: 10/10/13
31 // Author: T.Koi
32 //
33 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34 //
35 // HISTORY
36 // Based on G4IsomerTable
38 //
39 #include "G4NuclideTable.hh"
40 
41 #include "G4ios.hh"
42 #include "globals.hh"
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include <iomanip>
46 #include <fstream>
47 #include <sstream>
48 
49 const G4double G4NuclideTable::levelTolerance = 1.0*eV;
50 // const G4double G4NuclideTable::levelTolerance = 1.0e-3*eV;
51 // torelance for excitation energy
52 
53 
55 G4NuclideTable* G4NuclideTable::GetInstance() {
56  static G4NuclideTable instance;
57  return &instance;
58 }
59 
61 G4NuclideTable::G4NuclideTable()
62  :G4VIsotopeTable("Isomer"),
63  threshold_of_half_life(1000.0*ns),
64  fUserDefinedList(NULL),
65  fIsotopeList(0)
66 {
67  SetVerboseLevel(G4ParticleTable::GetParticleTable()->GetVerboseLevel());
68  FillHardCodeList();
69 }
70 
72 G4NuclideTable::~G4NuclideTable()
73 {
74  if (fIsotopeList!=0) {
75  for (size_t i = 0 ; i<fIsotopeList->size(); i++) {
76  delete (*fIsotopeList)[i];
77  }
78  fIsotopeList->clear();
79  delete fIsotopeList;
80  fIsotopeList = 0;
81  }
82 
83  for ( std::map< G4int , std::multimap< G4double , G4IsotopeProperty* > >::iterator
84  it = map_pre_load_list.begin(); it != map_pre_load_list.end(); it++ ) {
85  it->second.clear();
86  }
87  map_pre_load_list.clear();
88 
89  for ( std::map< G4int , std::multimap< G4double , G4IsotopeProperty* > >::iterator
90  it = map_hard_code_list.begin(); it != map_hard_code_list.end(); it++ ) {
91  for ( std::multimap< G4double , G4IsotopeProperty* >::iterator
92  itt = it->second.begin(); itt != it->second.end(); itt++ ) {
93  delete itt->second;
94  }
95  it->second.clear();
96  }
97  map_hard_code_list.clear();
98 
99  for ( std::map< G4int , std::multimap< G4double , G4IsotopeProperty* > >::iterator
100  it = map_full_list.begin(); it != map_full_list.end(); it++ ) {
101  for ( std::multimap< G4double , G4IsotopeProperty* >::iterator
102  itt = it->second.begin(); itt != it->second.end(); itt++ ) {
103  delete itt->second;
104  }
105  it->second.clear();
106  }
107  map_full_list.clear();
108 }
109 
111 //
112 G4IsotopeProperty* G4NuclideTable::GetIsotope(G4int Z, G4int A, G4double E)
113 {
114 
115  G4IsotopeProperty* fProperty = 0;
116  G4int ionCode = 1000*Z + A;
117 
118  //Serching pre-load
119  //Note: isomer level is properly set only for pre_load_list.
120  if ( map_pre_load_list.find( ionCode ) != map_pre_load_list.end() ) {
121 
122  std::multimap< G4double , G4IsotopeProperty* >::iterator lower_bound_itr =
123  map_pre_load_list.find( ionCode ) -> second.lower_bound ( E - levelTolerance/2 );
124 
125  //std::multimap< G4double , G4IsotopeProperty* >::iterator upper_bound_itr =
126  //map_pre_load_list.find( ionCode ) -> second.upper_bound ( E );
127 
128  G4double levelE = DBL_MAX;
129  if ( lower_bound_itr != map_pre_load_list.find( ionCode ) -> second.end() ) {
130  levelE = lower_bound_itr->first;
131  if ( levelE - levelTolerance/2 <= E && E < levelE + levelTolerance/2 ) {
132  return lower_bound_itr->second; // found
133  }
134  }
135  }
136 
137  //Searching hard-code
138  if ( map_hard_code_list.find( ionCode ) != map_hard_code_list.end() ) {
139  std::multimap< G4double , G4IsotopeProperty* >::iterator lower_bound_itr =
140  map_hard_code_list.find( ionCode ) -> second.lower_bound ( E - levelTolerance/2 );
141 
142  //std::multimap< G4double , G4IsotopeProperty* >::iterator upper_bound_itr =
143  //map_pre_load_list.find( ionCode ) -> second.upper_bound ( E );
144 
145  G4double levelE = DBL_MAX;
146  if ( lower_bound_itr != map_hard_code_list.find( ionCode ) -> second.end() ) {
147  levelE = lower_bound_itr->first;
148  if ( levelE - levelTolerance/2 <= E && E < levelE + levelTolerance/2 ) {
149  return lower_bound_itr->second; // found
150  }
151  }
152  }
153 
154  //Searching big-list
155  char* path = getenv("G4ENSDFSTATEDATA");
156 
157  if ( !path ) {
158  return fProperty; // not found;
159  }
160 
161  if ( map_full_list.find( ionCode ) == map_full_list.end() ) {
162 
163  std::multimap<G4double, G4IsotopeProperty*> aMultiMap;
164  map_full_list.insert( std::pair< G4int , std::multimap< G4double , G4IsotopeProperty* > > ( ionCode , aMultiMap ) );
165 
166  std::fstream ifs;
167  G4String filename(path);
168  filename += "/ENSDFSTATE.dat";
169  ifs.open( filename.c_str() );
170 
171  G4bool reading_target = false;
172 
173  G4int ionZ;
174  G4int ionA;
175  G4double ionE;
176  G4double ionLife;
177  G4int ionJ;
178  G4double ionMu;
179 
180  ifs >> ionZ >> ionA >> ionE >> ionLife >> ionJ >> ionMu;
181 
182  while ( ifs.good() ) {
183 
184  if ( ionZ == Z && ionA == A ) {
185 
186  reading_target = true;
187 
188  ionE *= keV;
189  ionLife *= ns;
190  ionMu *= (joule/tesla);
191 
192  G4IsotopeProperty* property = new G4IsotopeProperty();
193 
194  G4int iLevel=9;
195  property->SetAtomicNumber(ionZ);
196  property->SetAtomicMass(ionA);
197  property->SetIsomerLevel(iLevel);
198  property->SetEnergy(ionE);
199  property->SetiSpin(ionJ);
200  property->SetLifeTime(ionLife);
201  property->SetMagneticMoment(ionMu);
202 
203  map_full_list.find ( ionCode ) -> second.insert( std::pair< G4double, G4IsotopeProperty* >( ionE , property ) );
204 
205  } else if ( reading_target == true ) {
206  ifs.close();
207  break;
208  }
209 
210  ifs >> ionZ >> ionA >> ionE >> ionLife >> ionJ >> ionMu;
211  }
212 
213  ifs.close();
214  }
215 
216 
217  if ( map_full_list.find( ionCode ) != map_full_list.end() ) {
218 
219  std::multimap< G4double , G4IsotopeProperty* >::iterator lower_bound_itr =
220  map_full_list.find( ionCode ) -> second.lower_bound ( E - levelTolerance/2 );
221 
222  //std::multimap< G4double , G4IsotopeProperty* >::iterator upper_bound_itr =
223  //map_full_list.find( ionCode ) -> second.upper_bound ( E - levelTolerance/2 );
224 
225  G4double levelE = DBL_MAX;
226  if ( lower_bound_itr != map_full_list.find( ionCode ) -> second.end() ) {
227  levelE = lower_bound_itr->first;
228  if ( levelE - levelTolerance/2 < E && E < levelE + levelTolerance/2 ) {
229  return lower_bound_itr->second; // found
230  }
231  }
232  }
233 
234  return fProperty; // not found;
235 }
236 
238 G4IsotopeProperty*
239  G4NuclideTable::GetIsotopeByIsoLvl(G4int Z, G4int A, G4int lvl)
240 {
241  if(lvl==0) return GetIsotope(Z,A,0.0);
242  return (G4IsotopeProperty*)0;
243 }
244 
246 void G4NuclideTable::FillHardCodeList()
247 {
248  for (size_t i=0; i<nEntries_ground_state; i++) {
249 
250  G4int ionZ = (G4int)groundStateTable[i][idxZ];
251  G4int ionA = (G4int)groundStateTable[i][idxA];
252  G4int lvl = 0; // ground state
253  G4double ionE = groundStateTable[i][idxEnergy]*keV;
254  G4double ionLife = groundStateTable[i][idxLife]*ns;
255  G4int ionJ = (G4int)(groundStateTable[i][idxSpin]);
256  G4double ionMu = groundStateTable[i][idxMu]*(joule/tesla);
257 
258  G4int ionCode = 1000*ionZ + ionA;
259 
260  G4IsotopeProperty* fProperty = new G4IsotopeProperty();
261 
262  // Set Isotope Property
263  fProperty->SetAtomicNumber(ionZ);
264  fProperty->SetAtomicMass(ionA);
265  fProperty->SetIsomerLevel(lvl);
266  fProperty->SetEnergy(ionE);
267  fProperty->SetiSpin(ionJ);
268  fProperty->SetLifeTime(ionLife);
269  fProperty->SetDecayTable(0);
270  fProperty->SetMagneticMoment(ionMu);
271 
272  if ( map_hard_code_list.find ( ionCode ) == map_hard_code_list.end() ) {
273  std::multimap<G4double, G4IsotopeProperty*> aMultiMap;
274  map_hard_code_list.insert( std::pair< G4int , std::multimap< G4double , G4IsotopeProperty* > > ( ionCode , aMultiMap ) );
275  }
276  map_hard_code_list.find ( ionCode ) -> second.insert( std::pair< G4double, G4IsotopeProperty* >( ionE , fProperty ) );
277 
278  }
279 
280  for (size_t i=0; i<nEntries_excite_state; i++) {
281 
282  G4int ionZ = (G4int)exciteStateTable[i][idxZ];
283  G4int ionA = (G4int)exciteStateTable[i][idxA];
284  G4double ionE = exciteStateTable[i][idxEnergy]*keV;
285  G4double ionLife = exciteStateTable[i][idxLife]*ns;
286  G4int ionJ = (G4int)(exciteStateTable[i][idxSpin]);
287  G4double ionMu = exciteStateTable[i][idxMu]*(joule/tesla);
288 
289  G4int ionCode = 1000*ionZ + ionA;
290 
291  G4IsotopeProperty* fProperty = new G4IsotopeProperty();
292 
293  // Set Isotope Property
294  fProperty->SetAtomicNumber(ionZ);
295  fProperty->SetAtomicMass(ionA);
296  fProperty->SetIsomerLevel(9);
297  fProperty->SetEnergy(ionE);
298  fProperty->SetiSpin(ionJ);
299  fProperty->SetLifeTime(ionLife);
300  fProperty->SetDecayTable(0);
301  fProperty->SetMagneticMoment(ionMu);
302 
303  if ( map_hard_code_list.find ( ionCode ) == map_hard_code_list.end() ) {
304  std::multimap<G4double, G4IsotopeProperty*> aMultiMap;
305  map_hard_code_list.insert( std::pair< G4int , std::multimap< G4double , G4IsotopeProperty* > > ( ionCode , aMultiMap ) );
306  }
307  map_hard_code_list.find ( ionCode ) -> second.insert( std::pair< G4double, G4IsotopeProperty* >( ionE , fProperty ) );
308 
309  }
310 }
311 
313 void G4NuclideTable::GenerateNuclide()
314 {
315 
316  if( fIsotopeList !=0 ) return;
317  fIsotopeList = new G4IsotopeList();
318 
319  for (size_t i=0; i<nEntries_ground_state; i++) {
320 
321  G4int ionZ = (G4int)groundStateTable[i][idxZ];
322  G4int ionA = (G4int)groundStateTable[i][idxA];
323  G4int lvl = 0; // ground state
324  G4double ionE = groundStateTable[i][idxEnergy]*keV;
325  G4double ionLife = groundStateTable[i][idxLife]*ns;
326  G4int ionJ = (G4int)(groundStateTable[i][idxSpin]);
327  G4double ionMu = groundStateTable[i][idxMu]*(joule/tesla);
328 
329  if ( ionLife < 0.0 || ionLife*std::log(2.0) > threshold_of_half_life ) {
330 
331  G4IsotopeProperty* fProperty = new G4IsotopeProperty();
332 
333  // Set Isotope Property
334  fProperty->SetAtomicNumber(ionZ);
335  fProperty->SetAtomicMass(ionA);
336  fProperty->SetIsomerLevel(lvl);
337  fProperty->SetEnergy(ionE);
338  fProperty->SetiSpin(ionJ);
339  fProperty->SetLifeTime(ionLife);
340  fProperty->SetDecayTable(0);
341  fProperty->SetMagneticMoment(ionMu);
342 
343  //G4cout << ionZ << " " << ionA << " " << lvl << " " << ionE/keV << " [keV]" << G4endl;
344  fIsotopeList->push_back(fProperty);
345 
346  G4int ionCode = 1000*ionZ + ionA;
347  if ( map_pre_load_list.find ( ionCode ) == map_pre_load_list.end() ) {
348  std::multimap<G4double, G4IsotopeProperty*> aMultiMap;
349  map_pre_load_list.insert( std::pair< G4int , std::multimap< G4double , G4IsotopeProperty* > > ( ionCode , aMultiMap ) );
350  }
351  map_pre_load_list.find ( ionCode ) -> second.insert( std::pair< G4double, G4IsotopeProperty* >( ionE , fProperty ) );
352 
353  }
354  }
355 
356  if ( threshold_of_half_life >= 1.0*ns ) {
357 
358  G4int ionCode=0;
359  G4int iLevel=0;
360  G4double previousE=0.0;
361 
362  for (size_t i=0; i<nEntries_excite_state; i++) {
363 
364  G4int ionZ = (G4int)exciteStateTable[i][idxZ];
365  G4int ionA = (G4int)exciteStateTable[i][idxA];
366  if ( ionCode != 1000*ionZ + ionA ) {
367  previousE=0.0;
368  iLevel = 0;
369  ionCode = 1000*ionZ + ionA;
370  }
371 
372  G4double ionE = exciteStateTable[i][idxEnergy]*keV;
373  G4double ionLife = exciteStateTable[i][idxLife]*ns;
374  G4int ionJ = (G4int)(exciteStateTable[i][idxSpin]);
375  G4double ionMu = exciteStateTable[i][idxMu]*(joule/tesla);
376 
377  if (( ionLife < 0.0 || ionLife > threshold_of_half_life )
378  && (ionE > levelTolerance+previousE)) {
379  previousE = ionE;
380  iLevel++;
381  if ( iLevel > 9 ) iLevel=9;
382  //G4cout << ionZ << " " << ionA << " " << iLevel << " " << ionE/keV << " [keV]" << G4endl;
383 
384  G4IsotopeProperty* fProperty = new G4IsotopeProperty();
385 
386  // Set Isotope Property
387  fProperty->SetAtomicNumber(ionZ);
388  fProperty->SetAtomicMass(ionA);
389  fProperty->SetIsomerLevel(iLevel);
390  fProperty->SetEnergy(ionE);
391  fProperty->SetiSpin(ionJ);
392  fProperty->SetLifeTime(ionLife);
393  fProperty->SetDecayTable(0);
394  fProperty->SetMagneticMoment(ionMu);
395 
396  fIsotopeList->push_back(fProperty);
397 
398  if ( map_pre_load_list.find ( ionCode ) == map_pre_load_list.end() ) {
399  std::multimap<G4double, G4IsotopeProperty*> aMultiMap;
400  map_pre_load_list.insert( std::pair< G4int , std::multimap< G4double , G4IsotopeProperty* > > ( ionCode , aMultiMap ) );
401  }
402  map_pre_load_list.find ( ionCode ) -> second.insert( std::pair< G4double, G4IsotopeProperty* >( ionE , fProperty ) );
403 
404  }
405  }
406  } else {
407 
408  char* path = getenv("G4ENSDFSTATEDATA");
409 
410  if ( !path ) {
411  G4Exception("G4NuclideTable", "PART70000",
412  FatalException, "G4ENSDFSTATEDATA environment variable must be set");
413  }
414 
415  std::fstream ifs;
416  G4String filename(path);
417  filename += "/ENSDFSTATE.dat";
418 
419  ifs.open( filename.c_str() );
420 
421  if ( !ifs.good() ) {
422  G4Exception("G4NuclideTable", "PART70001",
423  FatalException, "ENSDFSTATE.dat is not found.");
424  }
425 
426 
427  G4int ionCode=0;
428  G4int iLevel=0;
429 
430  G4int ionZ;
431  G4int ionA;
432  G4double ionE;
433  G4double ionLife;
434  G4int ionJ;
435  G4double ionMu;
436 
437  ifs >> ionZ >> ionA >> ionE >> ionLife >> ionJ >> ionMu;
438 
439  while ( ifs.good() ) {
440 
441  if ( ionCode != 1000*ionZ + ionA ) {
442  iLevel = 0;
443  ionCode = 1000*ionZ + ionA;
444  }
445 
446  ionE *= keV;
447  ionLife *= ns;
448  ionMu *= (joule/tesla);
449 
450  //if ( ionLife == -1 || ionLife > threshold_of_half_life ) {
451  if ( ionLife*std::log(2.0) > threshold_of_half_life && ionE != 0 ) {
452 
453  iLevel++;
454  if ( iLevel > 9 ) iLevel=9;
455  //G4cout << ionZ << " " << ionA << " " << iLevel << " " << ionE/keV << " [keV]" << G4endl;
456 
457  G4IsotopeProperty* fProperty = new G4IsotopeProperty();
458 
459  // Set Isotope Property
460  fProperty->SetAtomicNumber(ionZ);
461  fProperty->SetAtomicMass(ionA);
462  fProperty->SetIsomerLevel(iLevel);
463  fProperty->SetEnergy(ionE);
464  fProperty->SetiSpin(ionJ);
465  fProperty->SetLifeTime(ionLife);
466  fProperty->SetDecayTable(0);
467  fProperty->SetMagneticMoment(ionMu);
468 
469  fIsotopeList->push_back(fProperty);
470 
471  if ( map_pre_load_list.find ( ionCode ) == map_pre_load_list.end() ) {
472  std::multimap<G4double, G4IsotopeProperty*> aMultiMap;
473  map_pre_load_list.insert( std::pair< G4int , std::multimap< G4double , G4IsotopeProperty* > > ( ionCode , aMultiMap ) );
474  }
475  map_pre_load_list.find ( ionCode ) -> second.insert( std::pair< G4double, G4IsotopeProperty* >( ionE , fProperty ) );
476 
477  }
478 
479  ifs >> ionZ >> ionA >> ionE >> ionLife >> ionJ >> ionMu;
480  }
481 
482  }
483 
484  if ( fUserDefinedList != NULL ) {
485  for ( G4IsotopeList::iterator it = fUserDefinedList->begin() ; it != fUserDefinedList->end() ; it++ ) {
486  fIsotopeList->push_back( *it );
487  }
488  }
489 
490 }
491 
492 void G4NuclideTable::AddState( G4int ionZ, G4int ionA, G4double ionE, G4double ionLife, G4int ionJ=0, G4double ionMu=0.0)
493 {
494  if ( fUserDefinedList == NULL ) fUserDefinedList = new G4IsotopeList();
495 
496  G4IsotopeProperty* fProperty = new G4IsotopeProperty();
497 
498  // Set Isotope Property
499  fProperty->SetAtomicNumber(ionZ);
500  fProperty->SetAtomicMass(ionA);
501  fProperty->SetIsomerLevel(9);
502  fProperty->SetEnergy(ionE);
503  fProperty->SetiSpin(ionJ);
504  fProperty->SetLifeTime(ionLife);
505  fProperty->SetDecayTable(0);
506  fProperty->SetMagneticMoment(ionMu);
507 
508  fUserDefinedList->push_back(fProperty);
509 
510 }
511 
G4IsotopeProperty::SetAtomicMass
void SetAtomicMass(G4int A)
Definition: G4IsotopeProperty.hh:133
G4NuclideTable::G4IsotopeList
std::vector< G4IsotopeProperty * > G4IsotopeList
Definition: G4NuclideTable.hh:70
G4int
int G4int
Definition: G4Types.hh:78
G4IsotopeProperty::SetiSpin
void SetiSpin(G4int J)
Definition: G4IsotopeProperty.hh:145
G4NuclideTable::GetIsotopeByIsoLvl
virtual G4IsotopeProperty * GetIsotopeByIsoLvl(G4int Z, G4int A, G4int lvl=0)
Definition: G4NuclideTable.cc:239
G4IsotopeProperty::SetMagneticMoment
void SetMagneticMoment(G4double M)
Definition: G4IsotopeProperty.hh:157
G4IsotopeProperty::SetLifeTime
void SetLifeTime(G4double T)
Definition: G4IsotopeProperty.hh:192
Z
Float_t Z
Definition: plot.C:39
G4bool
bool G4bool
Definition: G4Types.hh:79
G4NuclideTable::GetIsotope
virtual G4IsotopeProperty * GetIsotope(G4int Z, G4int A, G4double E)
Definition: G4NuclideTable.cc:112
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhysicalConstants.hh
G4IsotopeProperty::SetEnergy
void SetEnergy(G4double E)
Definition: G4IsotopeProperty.hh:169
G4IsotopeProperty::SetDecayTable
void SetDecayTable(G4DecayTable *table)
Definition: G4IsotopeProperty.hh:204
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:95
G4VIsotopeTable::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VIsotopeTable.hh:105
G4ios.hh
FatalException
Definition: G4ExceptionSeverity.hh:60
G4NuclideTable::AddState
void AddState(G4int, G4int, G4double, G4double, G4int, G4double)
Definition: G4NuclideTable.cc:492
G4NuclideTable::GetInstance
static G4NuclideTable * GetInstance()
Definition: G4NuclideTable.cc:55
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4VIsotopeTable::SetVerboseLevel
void SetVerboseLevel(G4int level)
Definition: G4VIsotopeTable.hh:111
G4NuclideTable::~G4NuclideTable
virtual ~G4NuclideTable()
Definition: G4NuclideTable.cc:72
G4IsotopeProperty::SetAtomicNumber
void SetAtomicNumber(G4int Z)
Definition: G4IsotopeProperty.hh:121
ns
#define ns
Definition: xmlparse.cc:597
G4NuclideTable.hh
G4IsotopeProperty::SetIsomerLevel
void SetIsomerLevel(G4int level)
Definition: G4IsotopeProperty.hh:180
Generated on Sat Dec 14 2013 14:34:21 for Geant4_10 by   doxygen 1.8.5