| Geant4 Cross Reference |
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
49 #include "G4ParticleHPIsoProbabilityTable_NJOY.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4DynamicParticle.hh"
52 #include "G4ParticleHPVector.hh"
53 #include "Randomize.hh"
54 #include "G4NucleiProperties.hh"
55 #include "G4ReactionProduct.hh"
56 #include "G4ParticleHPManager.hh"
57 #include "G4ParticleHPChannel.hh"
58 #include "G4ParticleHPChannelList.hh"
59 #include "G4Nucleus.hh"
60 #include "G4Element.hh"
61
62 #include <string>
63 #include <sstream>
64
65 ///--------------------------------------------------------------------------------------
66 G4ParticleHPIsoProbabilityTable_NJOY::G4ParticleHPIsoProbabilityTable_NJOY() {
67 tableOrder = 0;
68 lssf_flag = -1;
69 }
70
71 ///--------------------------------------------------------------------------------------
72 G4ParticleHPIsoProbabilityTable_NJOY::~G4ParticleHPIsoProbabilityTable_NJOY() {}
73
74 ///--------------------------------------------------------------------------------------
75 void 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;
99 theEnergies = new G4ParticleHPVector(nEnergies);
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 ///--------------------------------------------------------------------------------------
131 G4double G4ParticleHPIsoProbabilityTable_NJOY::GetCorrelatedIsoCrossSectionPT( const G4DynamicParticle* dp, G4int MTnumber,
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;
143 if ( kineticEnergy < Emin || kineticEnergy > Emax ) {
144 // if the kinetic energy is outside of the URR limits for the given isotope, it finds the smooth cross section
145 G4int indexEl = (G4int)ele->GetIndex();
146 G4int isotopeJ = -1; // index of isotope in the given element
147 for ( G4int j = 0; j < (G4int)ele->GetNumberOfIsotopes(); j++ ) {
148 if ( A == (G4int)ele->GetIsotope(j)->GetN() ) {
149 isotopeJ = j;
150 break;
151 }
152 }
153 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
154 G4double weightedelasticXS;
155 G4double weightedcaptureXS;
156 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
157 weightedelasticXS = (*G4ParticleHPManager::GetInstance()->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
158 weightedcaptureXS = (*G4ParticleHPManager::GetInstance()->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
159 } else {
160 weightedelasticXS = this->GetDopplerBroadenedElasticXS( dp, indexEl, isotopeJ );
161 weightedcaptureXS = this->GetDopplerBroadenedCaptureXS( dp, indexEl, isotopeJ );
162 }
163 xsela_cache[id] = weightedelasticXS / frac;
164 xscap_cache[id] = weightedcaptureXS / frac;
165 if ( Z < 88 ) {
166 xsfiss_cache[id] = 0.0;
167 } else {
168 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
169 G4double weightedfissionXS = (*G4ParticleHPManager::GetInstance()->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
170 xsfiss_cache[id] = weightedfissionXS / frac;
171 } else {
172 G4double weightedfissionXS = this->GetDopplerBroadenedFissionXS( dp, indexEl, isotopeJ );
173 xsfiss_cache[id] = weightedfissionXS / frac;
174 }
175 }
176 } else if ( lssf_flag == 0 ) {
177 G4int indexE = theEnergies->GetEnergyIndex( kineticEnergy );
178 std::vector< G4double >* theProbability1 = theProbabilities->at(indexE - 1);
179 std::vector< G4double >* theProbability2 = theProbabilities->at(indexE);
180 G4double ene1 = theEnergies->GetEnergy(indexE - 1);
181 G4double ene2 = theEnergies->GetEnergy(indexE);
182 G4double rand = random_number;
183 G4int indexP1;
184 G4int indexP2;
185 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
186 if ( rand <= theProbability1->at(indexP1) ) break;
187 }
188 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
189 if ( rand <= theProbability2->at(indexP2) ) break;
190 }
191 G4double ela1 = theElasticData->at(indexE - 1)->at(indexP1);
192 G4double ela2 = theElasticData->at(indexE)->at(indexP2);
193 G4double cap1 = theCaptureData->at(indexE - 1)->at(indexP1);
194 G4double cap2 = theCaptureData->at(indexE)->at(indexP2);
195 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * barn;
196 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * barn;
197 if ( Z < 88 ) {
198 xsfiss_cache[id] = 0.0;
199 } else {
200 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
201 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
202 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * barn;
203 }
204 } else if ( lssf_flag == 1 ) {
205 G4int indexE = theEnergies->GetEnergyIndex( kineticEnergy );
206 std::vector< G4double >* theProbability1 = theProbabilities->at(indexE - 1);
207 std::vector< G4double >* theProbability2 = theProbabilities->at(indexE);
208 G4double ene1 = theEnergies->GetEnergy(indexE - 1);
209 G4double ene2 = theEnergies->GetEnergy(indexE);
210 G4double rand = random_number;
211 G4int indexP1;
212 G4int indexP2;
213 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
214 if ( rand <= theProbability1->at(indexP1) ) break;
215 }
216 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
217 if ( rand <= theProbability2->at(indexP2) ) break;
218 }
219 G4int indexEl = (G4int)ele->GetIndex();
220 G4int isotopeJ = -1; // index of isotope in the given element
221 for ( G4int j = 0; j < (G4int)ele->GetNumberOfIsotopes(); j++ ) {
222 if ( A == (G4int)ele->GetIsotope(j)->GetN() ) {
223 isotopeJ = j;
224 break;
225 }
226 }
227 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
228 G4double weightedelasticXS;
229 G4double weightedcaptureXS;
230 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
231 weightedelasticXS = (*G4ParticleHPManager::GetInstance()->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
232 weightedcaptureXS = (*G4ParticleHPManager::GetInstance()->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
233 } else {
234 weightedelasticXS = this->GetDopplerBroadenedElasticXS( dp, indexEl, isotopeJ );
235 weightedcaptureXS = this->GetDopplerBroadenedCaptureXS( dp, indexEl, isotopeJ );
236 }
237 G4double ela1 = theElasticData->at(indexE - 1)->at(indexP1);
238 G4double ela2 = theElasticData->at(indexE)->at(indexP2);
239 G4double cap1 = theCaptureData->at(indexE - 1)->at(indexP1);
240 G4double cap2 = theCaptureData->at(indexE)->at(indexP2);
241 G4double elasticXS = weightedelasticXS / frac;
242 G4double captureXS = weightedcaptureXS / frac;
243 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * elasticXS;
244 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * captureXS;
245 if ( Z < 88 ) {
246 xsfiss_cache[id] = 0.0;
247 } else {
248 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
249 G4double weightedfissionXS = (*G4ParticleHPManager::GetInstance()->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
250 G4double fissionXS = weightedfissionXS / frac;
251 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
252 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
253 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
254 } else {
255 G4double weightedfissionXS = this->GetDopplerBroadenedFissionXS( dp, indexEl, isotopeJ );
256 G4double fissionXS = weightedfissionXS / frac;
257 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
258 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
259 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
260 }
261 }
262 }
263 if ( MTnumber == 2 ) { // elastic cross section
264 return xsela_cache[id];
265 } else if ( MTnumber == 102 ) { // radiative capture cross section
266 return xscap_cache[id];
267 } else if ( MTnumber == 18 ) { // fission cross section
268 return xsfiss_cache[id];
269 } else {
270 G4cout << "Reaction was not found, returns 0." << G4endl;
271 return 0;
272 }
273 }
274
275 ///--------------------------------------------------------------------------------------
276 G4double G4ParticleHPIsoProbabilityTable_NJOY::GetIsoCrossSectionPT( const G4DynamicParticle* dp, G4int MTnumber,
277 const G4Element* ele, G4double& kineticEnergy, std::map< std::thread::id, G4double >& random_number_cache, std::thread::id& id ) {
278 energy_cache[id] = kineticEnergy;
279 if ( kineticEnergy < Emin || kineticEnergy > Emax ) {
280 // if the kinetic energy is outside of the URR limits for the given isotope, it finds the smooth cross section
281 G4int indexEl = (G4int)ele->GetIndex();
282 G4int isotopeJ = -1; // index of isotope in the given element
283 for ( G4int j = 0; j < (G4int)ele->GetNumberOfIsotopes(); j++ ) {
284 if ( A == (G4int)ele->GetIsotope(j)->GetN() ) {
285 isotopeJ = j;
286 break;
287 }
288 }
289 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
290 G4double weightedelasticXS;
291 G4double weightedcaptureXS;
292 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
293 weightedelasticXS = (*G4ParticleHPManager::GetInstance()->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
294 weightedcaptureXS = (*G4ParticleHPManager::GetInstance()->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
295 } else {
296 weightedelasticXS = this->GetDopplerBroadenedElasticXS( dp, indexEl, isotopeJ );
297 weightedcaptureXS = this->GetDopplerBroadenedCaptureXS( dp, indexEl, isotopeJ );
298 }
299 xsela_cache[id] = weightedelasticXS / frac;
300 xscap_cache[id] = weightedcaptureXS / frac;
301 if ( Z < 88 ) {
302 xsfiss_cache[id] = 0.0;
303 } else {
304 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
305 G4double weightedfissionXS = (*G4ParticleHPManager::GetInstance()->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
306 xsfiss_cache[id] = weightedfissionXS / frac;
307 } else {
308 G4double weightedfissionXS = this->GetDopplerBroadenedFissionXS( dp, indexEl, isotopeJ );
309 xsfiss_cache[id] = weightedfissionXS / frac;
310 }
311 }
312 } else if ( lssf_flag == 0 ) {
313 G4int indexE = theEnergies->GetEnergyIndex(kineticEnergy);
314 std::vector< G4double >* theProbability1 = theProbabilities->at(indexE - 1);
315 std::vector< G4double >* theProbability2 = theProbabilities->at(indexE);
316 G4double ene1 = theEnergies->GetEnergy(indexE - 1);
317 G4double ene2 = theEnergies->GetEnergy(indexE);
318 G4double rand = G4UniformRand();
319 random_number_cache[id] = rand;
320 G4int indexP1;
321 G4int indexP2;
322 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
323 if ( rand <= theProbability1->at(indexP1) ) break;
324 }
325 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
326 if ( rand <= theProbability2->at(indexP2) ) break;
327 }
328 G4double ela1 = theElasticData->at(indexE - 1)->at(indexP1);
329 G4double ela2 = theElasticData->at(indexE)->at(indexP2);
330 G4double cap1 = theCaptureData->at(indexE - 1)->at(indexP1);
331 G4double cap2 = theCaptureData->at(indexE)->at(indexP2);
332 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * barn;
333 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * barn;
334 if ( Z < 88 ) {
335 xsfiss_cache[id] = 0.0;
336 } else {
337 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
338 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
339 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * barn;
340 }
341 } else if ( lssf_flag == 1 ) {
342 G4int indexE = theEnergies->GetEnergyIndex( kineticEnergy );
343 std::vector< G4double >* theProbability1 = theProbabilities->at(indexE - 1);
344 std::vector< G4double >* theProbability2 = theProbabilities->at(indexE);
345 G4double ene1 = theEnergies->GetEnergy(indexE - 1);
346 G4double ene2 = theEnergies->GetEnergy(indexE);
347 G4double rand = G4UniformRand();
348 random_number_cache[id] = rand;
349 G4int indexP1;
350 G4int indexP2;
351 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
352 if ( rand <= theProbability1->at(indexP1) ) break;
353 }
354 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
355 if ( rand <= theProbability2->at(indexP2) ) break;
356 }
357 G4int indexEl = (G4int)ele->GetIndex();
358 G4int isotopeJ = -1; // index of isotope in the given element
359 for ( G4int j = 0; j < (G4int)ele->GetNumberOfIsotopes(); j++ ) {
360 if ( A == (G4int)ele->GetIsotope(j)->GetN() ) {
361 isotopeJ = j;
362 break;
363 }
364 }
365 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
366 G4double weightedelasticXS;
367 G4double weightedcaptureXS;
368 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
369 weightedelasticXS = (*G4ParticleHPManager::GetInstance()->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
370 weightedcaptureXS = (*G4ParticleHPManager::GetInstance()->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
371 } else {
372 weightedelasticXS = this->GetDopplerBroadenedElasticXS( dp, indexEl, isotopeJ );
373 weightedcaptureXS = this->GetDopplerBroadenedCaptureXS( dp, indexEl, isotopeJ );
374 }
375 G4double ela1 = theElasticData->at(indexE - 1)->at(indexP1);
376 G4double ela2 = theElasticData->at(indexE)->at(indexP2);
377 G4double cap1 = theCaptureData->at(indexE - 1)->at(indexP1);
378 G4double cap2 = theCaptureData->at(indexE)->at(indexP2);
379 G4double elasticXS = weightedelasticXS / frac;
380 G4double captureXS = weightedcaptureXS / frac;
381 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * elasticXS;
382 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * captureXS;
383 if ( Z < 88 ) {
384 xsfiss_cache[id] = 0.0;
385 } else {
386 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
387 G4double weightedfissionXS = (*G4ParticleHPManager::GetInstance()->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
388 G4double fissionXS = weightedfissionXS / frac;
389 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
390 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
391 xsfiss_cache[id] = theInt.Lin(kineticEnergy, ene1, ene2, fiss1, fiss2) * fissionXS;
392 } else {
393 G4double weightedfissionXS = this->GetDopplerBroadenedFissionXS( dp, indexEl, isotopeJ );
394 G4double fissionXS = weightedfissionXS / frac;
395 G4double fiss1 = theFissionData->at(indexE - 1)->at(indexP1);
396 G4double fiss2 = theFissionData->at(indexE)->at(indexP2);
397 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
398 }
399 }
400 }
401 if ( MTnumber == 2 ) { // elastic cross section
402 return xsela_cache[id];
403 } else if ( MTnumber == 102 ) { // radiative capture cross section
404 return xscap_cache[id];
405 } else if ( MTnumber == 18 ) { // fission cross section
406 return xsfiss_cache[id];
407 } else {
408 G4cout << "Reaction was not found, returns 0." << G4endl;
409 return 0;
410 }
411 }
412