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);
85 std::string strm = std::to_string(
m);
89 std::istringstream theData( std::ios::in );
91 if ( theData.good() ) {
95 theData >> emin >> emax;
98 theData >>
nEnergies >> tableOrder >> lssf_flag;
104 G4double probability, total, elastic, capture, fission;
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 );
124 G4cout <<
"Probability tables found and succesfully read from " <<
Emin / keV <<
" keV to " <<
Emax / keV <<
" keV." <<
G4endl;
126 G4cout <<
"No probability tables found for this isotope and temperature, smooth cross section will be used instead." <<
G4endl;
134 if ( MTnumber == 2 ) {
136 }
else if ( MTnumber == 102 ) {
138 }
else if ( MTnumber == 18 ) {
143 if ( kineticEnergy < Emin || kineticEnergy >
Emax ) {
147 for (
G4int j = 0; j < (
G4int)ele->GetNumberOfIsotopes(); j++ ) {
148 if (
A == (
G4int)ele->GetIsotope(j)->GetN() ) {
153 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
176 }
else if ( lssf_flag == 0 ) {
178 std::vector< G4double >* theProbability1 =
theProbabilities->at(indexE - 1);
185 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
186 if ( rand <= theProbability1->at(indexP1) )
break;
188 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
189 if ( rand <= theProbability2->at(indexP2) )
break;
195 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * barn;
196 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * barn;
202 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * barn;
204 }
else if ( lssf_flag == 1 ) {
206 std::vector< G4double >* theProbability1 =
theProbabilities->at(indexE - 1);
213 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
214 if ( rand <= theProbability1->at(indexP1) )
break;
216 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
217 if ( rand <= theProbability2->at(indexP2) )
break;
221 for (
G4int j = 0; j < (
G4int)ele->GetNumberOfIsotopes(); j++ ) {
222 if (
A == (
G4int)ele->GetIsotope(j)->GetN() ) {
227 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
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;
250 G4double fissionXS = weightedfissionXS / frac;
253 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
256 G4double fissionXS = weightedfissionXS / frac;
259 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
263 if ( MTnumber == 2 ) {
265 }
else if ( MTnumber == 102 ) {
267 }
else if ( MTnumber == 18 ) {
270 G4cout <<
"Reaction was not found, returns 0." <<
G4endl;
277 const G4Element* ele,
G4double& kineticEnergy, std::map< std::thread::id, G4double >& random_number_cache, std::thread::id&
id ) {
279 if ( kineticEnergy < Emin || kineticEnergy >
Emax ) {
283 for (
G4int j = 0; j < (
G4int)ele->GetNumberOfIsotopes(); j++ ) {
284 if (
A == (
G4int)ele->GetIsotope(j)->GetN() ) {
289 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
312 }
else if ( lssf_flag == 0 ) {
314 std::vector< G4double >* theProbability1 =
theProbabilities->at(indexE - 1);
319 random_number_cache[id] = rand;
322 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
323 if ( rand <= theProbability1->at(indexP1) )
break;
325 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
326 if ( rand <= theProbability2->at(indexP2) )
break;
332 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * barn;
333 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * barn;
339 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * barn;
341 }
else if ( lssf_flag == 1 ) {
343 std::vector< G4double >* theProbability1 =
theProbabilities->at(indexE - 1);
348 random_number_cache[id] = rand;
351 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
352 if ( rand <= theProbability1->at(indexP1) )
break;
354 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
355 if ( rand <= theProbability2->at(indexP2) )
break;
359 for (
G4int j = 0; j < (
G4int)ele->GetNumberOfIsotopes(); j++ ) {
360 if (
A == (
G4int)ele->GetIsotope(j)->GetN() ) {
365 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
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;
388 G4double fissionXS = weightedfissionXS / frac;
391 xsfiss_cache[id] = theInt.Lin(kineticEnergy, ene1, ene2, fiss1, fiss2) * fissionXS;
394 G4double fissionXS = weightedfissionXS / frac;
397 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
401 if ( MTnumber == 2 ) {
403 }
else if ( MTnumber == 102 ) {
405 }
else if ( MTnumber == 18 ) {
408 G4cout <<
"Reaction was not found, returns 0." <<
G4endl;
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 GetDataStream(const G4String &, std::istringstream &iss)