ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NuclideTable.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NuclideTable.cc
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 #include "G4NuclideTableMessenger.hh"
41 
42 #include "G4ios.hh"
43 #include "G4String.hh"
44 #include "globals.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include <iomanip>
48 #include <fstream>
49 #include <sstream>
50 
51 // const G4double G4NuclideTable::levelTolerance =
52 // torelance for excitation energy
53 
54 
56 G4NuclideTable* G4NuclideTable::GetInstance() {
57  static G4NuclideTable instance;
58  return &instance;
59 }
60 
62 G4NuclideTable::G4NuclideTable()
63  :G4VIsotopeTable("Isomer"),
64  threshold_of_half_life(1000.0*ns),
65  minimum_threshold_of_half_life(DBL_MAX),
66  fUserDefinedList(nullptr),
67  fIsotopeList(nullptr),
68  flevelTolerance(1.0*eV)
69 {
70  //SetVerboseLevel(G4ParticleTable::GetParticleTable()->GetVerboseLevel());
71  //FillHardCodeList();
72  fMessenger = new G4NuclideTableMessenger(this);
73  fIsotopeList = new G4IsotopeList();
74  GenerateNuclide();
75 }
76 
78 G4NuclideTable::~G4NuclideTable()
79 {
80 
81  for ( std::map< G4int , std::multimap< G4double , G4IsotopeProperty* > >::iterator
82  it = map_pre_load_list.begin(); it != map_pre_load_list.end(); it++ ) {
83  it->second.clear();
84  }
85  map_pre_load_list.clear();
86 
87 
88  for ( std::map< G4int , std::multimap< G4double , G4IsotopeProperty* > >::iterator
89  it = map_full_list.begin(); it != map_full_list.end(); it++ ) {
90  it->second.clear();
91  }
92  map_full_list.clear();
93 
94  if (fIsotopeList != nullptr) {
95  for (size_t i = 0 ; i<fIsotopeList->size(); i++) {
96  delete (*fIsotopeList)[i];
97  }
98  fIsotopeList->clear();
99  delete fIsotopeList;
100  fIsotopeList = nullptr;
101  }
102 
103 }
104 
106 //
107 G4IsotopeProperty* G4NuclideTable::GetIsotope(G4int Z, G4int A, G4double E,
108  G4Ions::G4FloatLevelBase flb)
109 {
110 
111  G4IsotopeProperty* fProperty = nullptr;
112 
113  // At first searching UserDefined
114  if ( fUserDefinedList ) {
115  for ( G4IsotopeList::iterator it = fUserDefinedList->begin() ; it != fUserDefinedList->end() ; it++ ) {
116 
117  if ( Z == (*it)->GetAtomicNumber() && A == (*it)->GetAtomicMass() ) {
118  G4double levelE = (*it)->GetEnergy();
119  if ( levelE - flevelTolerance/2 <= E && E < levelE + flevelTolerance/2 ) {
120  if( flb == (*it)->GetFloatLevelBase() )
121  { return *it; } //found
122  }
123  }
124  }
125  }
126 
127  //Serching pre-load
128  //Note: isomer level is properly set only for pre_load_list.
129  //
130  G4int ionCode = 1000*Z + A;
131  std::map< G4int , std::multimap< G4double , G4IsotopeProperty* > >::iterator itf = map_pre_load_list.find( ionCode );
132 
133  if ( itf != map_pre_load_list.end() ) {
134  std::multimap< G4double , G4IsotopeProperty* >::iterator lower_bound_itr = itf -> second.lower_bound ( E - flevelTolerance/2 );
135  G4double levelE = DBL_MAX;
136 /*
137  if ( lower_bound_itr != itf -> second.end() ) {
138  levelE = lower_bound_itr->first;
139  if ( levelE - flevelTolerance/2 <= E && E < levelE + flevelTolerance/2 ) {
140  if( flb == (lower_bound_itr->second)->GetFloatLevelBase() )
141  { return lower_bound_itr->second; } // found
142  }
143  }
144 */
145  while ( lower_bound_itr != itf -> second.end() ) {
146  levelE = lower_bound_itr->first;
147  if ( levelE - flevelTolerance/2 <= E && E < levelE + flevelTolerance/2 ) {
148  if ( flb == (lower_bound_itr->second)->GetFloatLevelBase() ) return lower_bound_itr->second; // found
149  } else {
150  break;
151  }
152  lower_bound_itr++;
153  }
154 
155  }
156 
157  return fProperty; // not found
158 }
159 
161 G4IsotopeProperty* G4NuclideTable::GetIsotopeByIsoLvl(G4int Z, G4int A, G4int lvl)
162 {
163  if(lvl==0) return GetIsotope(Z,A,0.0);
164  return nullptr;
165 }
166 
168 void G4NuclideTable::FillHardCodeList()
169 {
170 }
171 
173 void G4NuclideTable::GenerateNuclide()
174 {
175  if ( threshold_of_half_life < minimum_threshold_of_half_life ) {
176  // Need to update full list
177  char* path = std::getenv("G4ENSDFSTATEDATA");
178 
179  if ( !path ) {
180  G4Exception("G4NuclideTable", "PART70000",
181  FatalException, "G4ENSDFSTATEDATA environment variable must be set");
182  return;
183  }
184 
185  std::ifstream ifs;
186  G4String filename(path);
187  filename += "/ENSDFSTATE.dat";
188 
189  ifs.open( filename.c_str() );
190  if ( !ifs.good() ) {
191  G4Exception("G4NuclideTable", "PART70001",
192  FatalException, "ENSDFSTATE.dat is not found.");
193  return;
194  }
195 
196  G4int ionCode=0;
197  G4int iLevel=0;
198  G4int ionZ;
199  G4int ionA;
200  G4double ionE;
201  G4String ionFL;
202  G4double ionLife;
203  G4int ionJ;
204  G4double ionMu;
205 
206  ifs >> ionZ >> ionA >> ionE >> ionFL >> ionLife >> ionJ >> ionMu;
207 
208  while ( ifs.good() ) {// Loop checking, 09.08.2015, K.Kurashige
209  if ( ionCode != 1000*ionZ + ionA ) {
210  iLevel = 0;
211  ionCode = 1000*ionZ + ionA;
212  }
213 
214  ionE *= keV;
215  G4Ions::G4FloatLevelBase flb = StripFloatLevelBase( ionFL );
216  ionLife *= ns;
217  ionMu *= (joule/tesla);
218 
219  if ( ( ionE == 0 && flb == G4Ions::G4FloatLevelBase::no_Float ) //
220  || ( threshold_of_half_life <= ionLife*std::log(2.0) && ionLife*std::log(2.0) < minimum_threshold_of_half_life ) ) {
221  if ( ionE > 0 ) iLevel++;
222  if ( iLevel > 9 ) iLevel=9;
223 
224  G4IsotopeProperty* fProperty = new G4IsotopeProperty();
225  // Set Isotope Property
226  fProperty->SetAtomicNumber(ionZ);
227  fProperty->SetAtomicMass(ionA);
228  fProperty->SetIsomerLevel(iLevel);
229  fProperty->SetEnergy(ionE);
230  fProperty->SetiSpin(ionJ);
231  fProperty->SetLifeTime(ionLife);
232  fProperty->SetDecayTable(nullptr);
233  fProperty->SetMagneticMoment(ionMu);
234  fProperty->SetFloatLevelBase( flb );
235 
236  fIsotopeList->push_back(fProperty);
237 
238  std::map< G4int , std::multimap< G4double , G4IsotopeProperty* > >::iterator itf = map_full_list.find( ionCode );
239  if ( itf == map_full_list.end() ) {
240  std::multimap<G4double, G4IsotopeProperty*> aMultiMap;
241  itf = ( map_full_list.insert( std::pair< G4int , std::multimap< G4double , G4IsotopeProperty* > > ( ionCode , aMultiMap ) ) ).first;
242  }
243  itf -> second.insert( std::pair< G4double, G4IsotopeProperty* >( ionE , fProperty ) );
244  }
245 
246  ifs >> ionZ >> ionA >> ionE >> ionFL >> ionLife >> ionJ >> ionMu;
247  }
248 
249  minimum_threshold_of_half_life = threshold_of_half_life;
250 
251  }
252 
253 
254  // Clear current map
255  for ( std::map< G4int , std::multimap< G4double , G4IsotopeProperty* > >::iterator
256  it = map_pre_load_list.begin(); it != map_pre_load_list.end(); it++ ) {
257  it->second.clear();
258  }
259  map_pre_load_list.clear();
260 
261  // Build map based on current threshold value
262  for ( std::map< G4int , std::multimap< G4double , G4IsotopeProperty* > >::iterator
263  it = map_full_list.begin(); it != map_full_list.end(); it++ ) {
264 
265  G4int ionCode = it->first;
266  std::map< G4int , std::multimap< G4double , G4IsotopeProperty* > >::iterator itf = map_pre_load_list.find( ionCode );
267  if ( itf == map_pre_load_list.end() ) {
268  std::multimap<G4double, G4IsotopeProperty*> aMultiMap;
269  itf = ( map_pre_load_list.insert( std::pair< G4int , std::multimap< G4double , G4IsotopeProperty* > > ( ionCode , aMultiMap ) ) ).first;
270  }
271  G4int iLevel = 0;
272  for ( std::multimap< G4double , G4IsotopeProperty* >::iterator
273  itt = it->second.begin(); itt != it->second.end(); itt++ ) {
274 
275  G4double exEnergy = itt->first;
276  G4double meanLife = itt->second->GetLifeTime();
277  if ( exEnergy == 0.0
278  || meanLife*std::log(2.0) > threshold_of_half_life ) {
279  if ( itt->first != 0.0 ) iLevel++;
280  if ( iLevel > 9 ) iLevel=9;
281  itt->second->SetIsomerLevel( iLevel );
282  itf -> second.insert( std::pair< G4double, G4IsotopeProperty* >( exEnergy , itt->second ) );
283  }
284  }
285  }
286 
287 }
288 
289 void G4NuclideTable::AddState( G4int ionZ, G4int ionA, G4double ionE, G4double ionLife, G4int ionJ, G4double ionMu )
290 {
291  if ( G4Threading::IsMasterThread() ) {
292  G4int flbIndex = 0;
293  ionE = StripFloatLevelBase( ionE, flbIndex );
294  AddState(ionZ,ionA,ionE,flbIndex,ionLife,ionJ,ionMu);
295  }
296 }
297 
298 void G4NuclideTable::AddState( G4int ionZ, G4int ionA, G4double ionE,
299  G4int flbIndex, G4double ionLife, G4int ionJ, G4double ionMu )
300 {
301  if ( G4Threading::IsMasterThread() ) {
302 
303  if ( fUserDefinedList == NULL ) fUserDefinedList = new G4IsotopeList();
304 
305  G4IsotopeProperty* fProperty = new G4IsotopeProperty();
306 
307  // Set Isotope Property
308  fProperty->SetAtomicNumber(ionZ);
309  fProperty->SetAtomicMass(ionA);
310  fProperty->SetIsomerLevel(9);
311  fProperty->SetEnergy(ionE);
312  fProperty->SetiSpin(ionJ);
313  fProperty->SetLifeTime(ionLife);
314  fProperty->SetDecayTable(nullptr);
315  fProperty->SetMagneticMoment(ionMu);
316  fProperty->SetFloatLevelBase(flbIndex);
317 
318  fUserDefinedList->push_back(fProperty);
319  fIsotopeList->push_back(fProperty);
320 
321  }
322 }
323 
324 void G4NuclideTable::AddState( G4int ionZ, G4int ionA, G4double ionE,
325  G4Ions::G4FloatLevelBase flb, G4double ionLife, G4int ionJ, G4double ionMu )
326 {
327  if ( G4Threading::IsMasterThread() ) {
328  if ( fUserDefinedList == nullptr ) fUserDefinedList = new G4IsotopeList();
329 
330  G4IsotopeProperty* fProperty = new G4IsotopeProperty();
331  // Set Isotope Property
332  fProperty->SetAtomicNumber(ionZ);
333  fProperty->SetAtomicMass(ionA);
334  fProperty->SetIsomerLevel(9);
335  fProperty->SetEnergy(ionE);
336  fProperty->SetiSpin(ionJ);
337  fProperty->SetLifeTime(ionLife);
338  fProperty->SetDecayTable(0);
339  fProperty->SetMagneticMoment(ionMu);
340  fProperty->SetFloatLevelBase( flb );
341 
342  fUserDefinedList->push_back(fProperty);
343  fIsotopeList->push_back(fProperty);
344 
345  }
346 }
347 
348 //#include "G4Threading.hh"
349 void G4NuclideTable::SetThresholdOfHalfLife( G4double t )
350 {
351  if ( G4Threading::IsMasterThread() ) {
352  threshold_of_half_life=t;
353  GenerateNuclide();
354  }
355 }
356 
357 G4double G4NuclideTable::StripFloatLevelBase(G4double E, G4int& flbIndex)
358 {
359  G4double rem = std::fmod(E/(1.0E-3*eV),10.0);
360  flbIndex = int(rem);
361  return E-rem;
362 }
363 
364 G4Ions::G4FloatLevelBase G4NuclideTable::StripFloatLevelBase( G4String sFLB )
365 {
366  if ( sFLB.size() < 1 || 2 < sFLB.size() ) {
367  G4String text;
368  text += sFLB;
369  text += " is not valid indicator of G4Ions::G4FloatLevelBase. You may use a wrong version of ENSDFSTATE data. Please use G4ENSDFSTATE2.0 or later.";
370 
371  G4Exception( "G4NuclideTable" , "PART70002" ,
372  FatalException , text );
373  }
374  G4Ions::G4FloatLevelBase flb = noFloat;
375  if ( !(sFLB == '-') ) {
376  flb = G4Ions::FloatLevelBase( sFLB.back() );
377  }
378  return flb;
379 }