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