Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4ParticleHPIsoProbabilityTable_NJOY.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//
29// Geant4 source file
30//
31// File name: G4ParticleHPIsoProbabilityTable_NJOY.cc
32//
33// Authors: Marek Zmeskal (CTU, Czech Technical University in Prague, Czech Republic)
34// Loic Thulliez (CEA France)
35//
36// Creation date: 4 June 2024
37//
38// Description: Class for the probability table of the given isotope
39// and for the given temperature generated with NJOY.
40// It reads the files with probability tables and
41// finds the correct cross-section.
42//
43// Modifications:
44//
45// -------------------------------------------------------------------
46//
47//
48
50#include "G4SystemOfUnits.hh"
51#include "G4DynamicParticle.hh"
52#include "G4ParticleHPVector.hh"
53#include "Randomize.hh"
54#include "G4NucleiProperties.hh"
55#include "G4ReactionProduct.hh"
59#include "G4Nucleus.hh"
60#include "G4Element.hh"
61
62#include <string>
63#include <sstream>
64
65///--------------------------------------------------------------------------------------
70
71///--------------------------------------------------------------------------------------
73
74///--------------------------------------------------------------------------------------
75void G4ParticleHPIsoProbabilityTable_NJOY::Init( G4int theZ, G4int theA, G4int them, G4double theT, const G4String& dirName ) {
76 Z = theZ;
77 A = theA;
78 m = them;
79 T = theT;
80 G4cout << "The NJOY probability tables are being initialized for Z=" << Z << " A=" << A << " and T=" << T << " K." << G4endl;
81 std::string strZ = std::to_string(Z);
82 std::string strA = std::to_string(A);
83 filename = strZ + "_" + strA;
84 if ( m != 0 ) {
85 std::string strm = std::to_string(m);
86 filename += "_m" + strm;
87 }
88 G4String fullPathFileName = dirName + filename + "." + std::to_string( (G4int)(T) ) + ".pt";
89 std::istringstream theData( std::ios::in );
90 G4ParticleHPManager::GetInstance()->GetDataStream( fullPathFileName, theData );
91 if ( theData.good() ) {
92 G4double emin;
93 G4double emax;
94 G4double energy;
95 theData >> emin >> emax;
96 Emin = emin * eV;
97 Emax = emax * eV;
98 theData >> nEnergies >> tableOrder >> lssf_flag;
100 theProbabilities = new std::vector< std::vector< G4double >* >;
101 theElasticData = new std::vector< std::vector< G4double >* >;
102 theCaptureData = new std::vector< std::vector< G4double >* >;
103 theFissionData = new std::vector< std::vector< G4double >* >;
104 G4double probability, total, elastic, capture, fission;
105 for ( G4int i = 0; i < nEnergies; i++ ) {
106 theData >> energy;
107 theEnergies->SetEnergy( i, energy * eV );
108 std::vector< G4double >* vecprob = new std::vector< G4double >;
109 std::vector< G4double >* vecela = new std::vector< G4double >;
110 std::vector< G4double >* veccap = new std::vector< G4double >;
111 std::vector< G4double >* vecfiss = new std::vector< G4double >;
112 for ( G4int j = 0; j < tableOrder; j++ ) {
113 theData >> probability >> total >> elastic >> capture >> fission;
114 vecprob->push_back( probability );
115 vecela->push_back( elastic );
116 veccap->push_back( capture );
117 vecfiss->push_back( fission );
118 }
119 theProbabilities->push_back( vecprob );
120 theElasticData->push_back( vecela );
121 theCaptureData->push_back( veccap );
122 theFissionData->push_back( vecfiss );
123 }
124 G4cout << "Probability tables found and succesfully read from " << Emin / keV << " keV to " << Emax / keV << " keV." << G4endl;
125 } else {
126 G4cout << "No probability tables found for this isotope and temperature, smooth cross section will be used instead." << G4endl;
127 }
128}
129
130///--------------------------------------------------------------------------------------
132 const G4Element* ele, G4double& kineticEnergy, G4double& random_number, std::thread::id& id ) {
133 if ( kineticEnergy == energy_cache[id] ) {
134 if ( MTnumber == 2 ) { // elastic cross section
135 return xsela_cache[id];
136 } else if ( MTnumber == 102 ) { // radiative capture cross section
137 return xscap_cache[id];
138 } else if ( MTnumber == 18 ) { // fission cross section
139 return xsfiss_cache[id];
140 }
141 }
142 energy_cache[id] = kineticEnergy;
144
145 if ( kineticEnergy < Emin || kineticEnergy > Emax ) {
146 // if the kinetic energy is outside of the URR limits for the given isotope, it finds the smooth cross section
147 G4int indexEl = (G4int)ele->GetIndex();
148 G4int isotopeJ = -1; // index of isotope in the given element
149 for ( G4int j = 0; j < (G4int)ele->GetNumberOfIsotopes(); j++ ) {
150 if ( A == (G4int)ele->GetIsotope(j)->GetN() ) {
151 isotopeJ = j;
152 break;
153 }
154 }
155 if (isotopeJ == -1) { return 0.0; }
156 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
157 G4double weightedelasticXS;
158 G4double weightedcaptureXS;
159 if ( manHP->GetNeglectDoppler() ) {
160 weightedelasticXS = (*manHP->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
161 weightedcaptureXS = (*manHP->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
162 } else {
163 weightedelasticXS = GetDopplerBroadenedElasticXS( dp, indexEl, isotopeJ );
164 weightedcaptureXS = GetDopplerBroadenedCaptureXS( dp, indexEl, isotopeJ );
165 }
166 xsela_cache[id] = weightedelasticXS / frac;
167 xscap_cache[id] = weightedcaptureXS / frac;
168 if ( Z < 88 ) {
169 xsfiss_cache[id] = 0.0;
170 } else {
171 if ( manHP->GetNeglectDoppler() ) {
172 G4double weightedfissionXS = (*manHP->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
173 xsfiss_cache[id] = weightedfissionXS / frac;
174 } else {
175 G4double weightedfissionXS = GetDopplerBroadenedFissionXS( dp, indexEl, isotopeJ );
176 xsfiss_cache[id] = weightedfissionXS / frac;
177 }
178 }
179 } else if ( lssf_flag == 0 ) {
180 G4int indexE = theEnergies->GetEnergyIndex( kineticEnergy );
181 std::vector< G4double >* theProbability1 = theProbabilities->at(indexE - 1);
182 std::vector< G4double >* theProbability2 = theProbabilities->at(indexE);
183 G4double ene1 = theEnergies->GetEnergy(indexE - 1);
184 G4double ene2 = theEnergies->GetEnergy(indexE);
185 G4double rand = random_number;
186 G4int indexP1;
187 G4int indexP2;
188 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
189 if ( rand <= theProbability1->at(indexP1) ) break;
190 }
191 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
192 if ( rand <= theProbability2->at(indexP2) ) break;
193 }
194 G4double ela1 = theElasticData->at(indexE - 1)->at(indexP1);
195 G4double ela2 = theElasticData->at(indexE)->at(indexP2);
196 G4double cap1 = theCaptureData->at(indexE - 1)->at(indexP1);
197 G4double cap2 = theCaptureData->at(indexE)->at(indexP2);
198 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * barn;
199 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * barn;
200 if ( Z < 88 ) {
201 xsfiss_cache[id] = 0.0;
202 } else {
203 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
204 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
205 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * barn;
206 }
207 } else if ( lssf_flag == 1 ) {
208 G4int indexE = theEnergies->GetEnergyIndex( kineticEnergy );
209 std::vector< G4double >* theProbability1 = theProbabilities->at(indexE - 1);
210 std::vector< G4double >* theProbability2 = theProbabilities->at(indexE);
211 G4double ene1 = theEnergies->GetEnergy(indexE - 1);
212 G4double ene2 = theEnergies->GetEnergy(indexE);
213 G4double rand = random_number;
214 G4int indexP1;
215 G4int indexP2;
216 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
217 if ( rand <= theProbability1->at(indexP1) ) break;
218 }
219 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
220 if ( rand <= theProbability2->at(indexP2) ) break;
221 }
222 G4int indexEl = (G4int)ele->GetIndex();
223 G4int isotopeJ = -1; // index of isotope in the given element
224 for ( G4int j = 0; j < (G4int)ele->GetNumberOfIsotopes(); j++ ) {
225 if ( A == (G4int)ele->GetIsotope(j)->GetN() ) {
226 isotopeJ = j;
227 break;
228 }
229 }
230 if (isotopeJ == -1) { return 0.0; }
231 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
232 G4double weightedelasticXS;
233 G4double weightedcaptureXS;
234 if ( manHP->GetNeglectDoppler() ) {
235 weightedelasticXS = (*manHP->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
236 weightedcaptureXS = (*manHP->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
237 } else {
238 weightedelasticXS = GetDopplerBroadenedElasticXS( dp, indexEl, isotopeJ );
239 weightedcaptureXS = GetDopplerBroadenedCaptureXS( dp, indexEl, isotopeJ );
240 }
241 G4double ela1 = theElasticData->at(indexE - 1)->at(indexP1);
242 G4double ela2 = theElasticData->at(indexE)->at(indexP2);
243 G4double cap1 = theCaptureData->at(indexE - 1)->at(indexP1);
244 G4double cap2 = theCaptureData->at(indexE)->at(indexP2);
245 G4double elasticXS = weightedelasticXS / frac;
246 G4double captureXS = weightedcaptureXS / frac;
247 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * elasticXS;
248 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * captureXS;
249 if ( Z < 88 ) {
250 xsfiss_cache[id] = 0.0;
251 } else {
252 if ( manHP->GetNeglectDoppler() ) {
253 G4double weightedfissionXS = (*manHP->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
254 G4double fissionXS = weightedfissionXS / frac;
255 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
256 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
257 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
258 } else {
259 G4double weightedfissionXS = GetDopplerBroadenedFissionXS( dp, indexEl, isotopeJ );
260 G4double fissionXS = weightedfissionXS / frac;
261 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
262 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
263 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
264 }
265 }
266 }
267 if ( MTnumber == 2 ) { // elastic cross section
268 return xsela_cache[id];
269 } else if ( MTnumber == 102 ) { // radiative capture cross section
270 return xscap_cache[id];
271 } else if ( MTnumber == 18 ) { // fission cross section
272 return xsfiss_cache[id];
273 } else {
274 //G4cout << "Reaction was not found, returns 0." << G4endl;
275 return 0;
276 }
277}
278
279///--------------------------------------------------------------------------------------
281 const G4Element* ele, G4double& kineticEnergy, std::map< std::thread::id, G4double >& random_number_cache, std::thread::id& id ) {
282 energy_cache[id] = kineticEnergy;
283 if ( kineticEnergy < Emin || kineticEnergy > Emax ) {
284 // if the kinetic energy is outside of the URR limits for the given isotope, it finds the smooth cross section
285 G4int indexEl = (G4int)ele->GetIndex();
286 G4int isotopeJ = 0; // index of isotope in the given element
287 for ( G4int j = 0; j < (G4int)ele->GetNumberOfIsotopes(); j++ ) {
288 if ( A == (G4int)ele->GetIsotope(j)->GetN() ) {
289 isotopeJ = j;
290 break;
291 }
292 }
293 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
294 G4double weightedelasticXS;
295 G4double weightedcaptureXS;
296 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
297 weightedelasticXS = (*G4ParticleHPManager::GetInstance()->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
298 weightedcaptureXS = (*G4ParticleHPManager::GetInstance()->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
299 } else {
300 weightedelasticXS = this->GetDopplerBroadenedElasticXS( dp, indexEl, isotopeJ );
301 weightedcaptureXS = this->GetDopplerBroadenedCaptureXS( dp, indexEl, isotopeJ );
302 }
303 xsela_cache[id] = weightedelasticXS / frac;
304 xscap_cache[id] = weightedcaptureXS / frac;
305 if ( Z < 88 ) {
306 xsfiss_cache[id] = 0.0;
307 } else {
308 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
309 G4double weightedfissionXS = (*G4ParticleHPManager::GetInstance()->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
310 xsfiss_cache[id] = weightedfissionXS / frac;
311 } else {
312 G4double weightedfissionXS = this->GetDopplerBroadenedFissionXS( dp, indexEl, isotopeJ );
313 xsfiss_cache[id] = weightedfissionXS / frac;
314 }
315 }
316 } else if ( lssf_flag == 0 ) {
317 G4int indexE = theEnergies->GetEnergyIndex(kineticEnergy);
318 std::vector< G4double >* theProbability1 = theProbabilities->at(indexE - 1);
319 std::vector< G4double >* theProbability2 = theProbabilities->at(indexE);
320 G4double ene1 = theEnergies->GetEnergy(indexE - 1);
321 G4double ene2 = theEnergies->GetEnergy(indexE);
322 G4double rand = G4UniformRand();
323 random_number_cache[id] = rand;
324 G4int indexP1;
325 G4int indexP2;
326 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
327 if ( rand <= theProbability1->at(indexP1) ) break;
328 }
329 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
330 if ( rand <= theProbability2->at(indexP2) ) break;
331 }
332 G4double ela1 = theElasticData->at(indexE - 1)->at(indexP1);
333 G4double ela2 = theElasticData->at(indexE)->at(indexP2);
334 G4double cap1 = theCaptureData->at(indexE - 1)->at(indexP1);
335 G4double cap2 = theCaptureData->at(indexE)->at(indexP2);
336 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * barn;
337 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * barn;
338 if ( Z < 88 ) {
339 xsfiss_cache[id] = 0.0;
340 } else {
341 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
342 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
343 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * barn;
344 }
345 } else if ( lssf_flag == 1 ) {
346 G4int indexE = theEnergies->GetEnergyIndex( kineticEnergy );
347 std::vector< G4double >* theProbability1 = theProbabilities->at(indexE - 1);
348 std::vector< G4double >* theProbability2 = theProbabilities->at(indexE);
349 G4double ene1 = theEnergies->GetEnergy(indexE - 1);
350 G4double ene2 = theEnergies->GetEnergy(indexE);
351 G4double rand = G4UniformRand();
352 random_number_cache[id] = rand;
353 G4int indexP1;
354 G4int indexP2;
355 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
356 if ( rand <= theProbability1->at(indexP1) ) break;
357 }
358 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
359 if ( rand <= theProbability2->at(indexP2) ) break;
360 }
361 G4int indexEl = (G4int)ele->GetIndex();
362 G4int isotopeJ = 0; // index of isotope in the given element
363 for ( G4int j = 0; j < (G4int)ele->GetNumberOfIsotopes(); j++ ) {
364 if ( A == (G4int)ele->GetIsotope(j)->GetN() ) {
365 isotopeJ = j;
366 break;
367 }
368 }
369 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
370 G4double weightedelasticXS;
371 G4double weightedcaptureXS;
372 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
373 weightedelasticXS = (*G4ParticleHPManager::GetInstance()->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
374 weightedcaptureXS = (*G4ParticleHPManager::GetInstance()->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
375 } else {
376 weightedelasticXS = this->GetDopplerBroadenedElasticXS( dp, indexEl, isotopeJ );
377 weightedcaptureXS = this->GetDopplerBroadenedCaptureXS( dp, indexEl, isotopeJ );
378 }
379 G4double ela1 = theElasticData->at(indexE - 1)->at(indexP1);
380 G4double ela2 = theElasticData->at(indexE)->at(indexP2);
381 G4double cap1 = theCaptureData->at(indexE - 1)->at(indexP1);
382 G4double cap2 = theCaptureData->at(indexE)->at(indexP2);
383 G4double elasticXS = weightedelasticXS / frac;
384 G4double captureXS = weightedcaptureXS / frac;
385 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * elasticXS;
386 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * captureXS;
387 if ( Z < 88 ) {
388 xsfiss_cache[id] = 0.0;
389 } else {
390 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
391 G4double weightedfissionXS = (*G4ParticleHPManager::GetInstance()->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
392 G4double fissionXS = weightedfissionXS / frac;
393 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
394 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
395 xsfiss_cache[id] = theInt.Lin(kineticEnergy, ene1, ene2, fiss1, fiss2) * fissionXS;
396 } else {
397 G4double weightedfissionXS = this->GetDopplerBroadenedFissionXS( dp, indexEl, isotopeJ );
398 G4double fissionXS = weightedfissionXS / frac;
399 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
400 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
401 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
402 }
403 }
404 }
405 if ( MTnumber == 2 ) { // elastic cross section
406 return xsela_cache[id];
407 } else if ( MTnumber == 102 ) { // radiative capture cross section
408 return xscap_cache[id];
409 } else if ( MTnumber == 18 ) { // fission cross section
410 return xsfiss_cache[id];
411 } else {
412 G4cout << "Reaction was not found, returns 0." << G4endl;
413 return 0;
414 }
415}
double G4double
Definition G4Types.hh:83
int G4int
Definition G4Types.hh:85
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition Randomize.hh:52
G4double GetIsoCrossSectionPT(const G4DynamicParticle *, G4int, const G4Element *, G4double &, std::map< std::thread::id, G4double > &, std::thread::id &) override
G4double GetCorrelatedIsoCrossSectionPT(const G4DynamicParticle *, G4int, const G4Element *, G4double &, G4double &, std::thread::id &) override
void Init(G4int, G4int, G4int, G4double, const G4String &) override
std::map< std::thread::id, G4double > xsela_cache
std::map< std::thread::id, G4double > xscap_cache
std::map< std::thread::id, G4double > energy_cache
G4double GetDopplerBroadenedFissionXS(const G4DynamicParticle *, G4int, G4int)
std::vector< std::vector< G4double > * > * theElasticData
std::vector< std::vector< G4double > * > * theFissionData
std::map< std::thread::id, G4double > xsfiss_cache
std::vector< std::vector< G4double > * > * theCaptureData
G4double GetDopplerBroadenedCaptureXS(const G4DynamicParticle *, G4int, G4int)
G4double GetDopplerBroadenedElasticXS(const G4DynamicParticle *, G4int, G4int)
std::vector< std::vector< G4double > * > * theProbabilities
std::vector< G4ParticleHPChannel * > * GetFissionFinalStates() const
std::vector< G4ParticleHPChannel * > * GetElasticFinalStates() const
std::vector< G4ParticleHPChannel * > * GetCaptureFinalStates() const
void GetDataStream(const G4String &, std::istringstream &iss)
static G4ParticleHPManager * GetInstance()
std::string to_string(G4FermiAtomicMass mass)