80 G4cout <<
"The NJOY probability tables are being initialized for Z=" <<
Z <<
" A=" <<
A <<
" and T=" <<
T <<
" K." <<
G4endl;
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 ) {
145 if ( kineticEnergy < Emin || kineticEnergy >
Emax ) {
149 for (
G4int j = 0; j < (
G4int)ele->GetNumberOfIsotopes(); j++ ) {
150 if (
A == (
G4int)ele->GetIsotope(j)->GetN() ) {
155 if (isotopeJ == -1) {
return 0.0; }
156 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
159 if ( manHP->GetNeglectDoppler() ) {
160 weightedelasticXS = (*manHP->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
161 weightedcaptureXS = (*manHP->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
171 if ( manHP->GetNeglectDoppler() ) {
172 G4double weightedfissionXS = (*manHP->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
179 }
else if ( lssf_flag == 0 ) {
181 std::vector< G4double >* theProbability1 =
theProbabilities->at(indexE - 1);
188 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
189 if ( rand <= theProbability1->at(indexP1) )
break;
191 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
192 if ( rand <= theProbability2->at(indexP2) )
break;
198 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * barn;
199 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * barn;
205 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * barn;
207 }
else if ( lssf_flag == 1 ) {
209 std::vector< G4double >* theProbability1 =
theProbabilities->at(indexE - 1);
216 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
217 if ( rand <= theProbability1->at(indexP1) )
break;
219 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
220 if ( rand <= theProbability2->at(indexP2) )
break;
224 for (
G4int j = 0; j < (
G4int)ele->GetNumberOfIsotopes(); j++ ) {
225 if (
A == (
G4int)ele->GetIsotope(j)->GetN() ) {
230 if (isotopeJ == -1) {
return 0.0; }
231 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
234 if ( manHP->GetNeglectDoppler() ) {
235 weightedelasticXS = (*manHP->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
236 weightedcaptureXS = (*manHP->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
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;
252 if ( manHP->GetNeglectDoppler() ) {
253 G4double weightedfissionXS = (*manHP->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ );
254 G4double fissionXS = weightedfissionXS / frac;
257 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
260 G4double fissionXS = weightedfissionXS / frac;
263 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
267 if ( MTnumber == 2 ) {
269 }
else if ( MTnumber == 102 ) {
271 }
else if ( MTnumber == 18 ) {
281 const G4Element* ele,
G4double& kineticEnergy, std::map< std::thread::id, G4double >& random_number_cache, std::thread::id&
id ) {
283 if ( kineticEnergy < Emin || kineticEnergy >
Emax ) {
287 for (
G4int j = 0; j < (
G4int)ele->GetNumberOfIsotopes(); j++ ) {
288 if (
A == (
G4int)ele->GetIsotope(j)->GetN() ) {
293 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
316 }
else if ( lssf_flag == 0 ) {
318 std::vector< G4double >* theProbability1 =
theProbabilities->at(indexE - 1);
323 random_number_cache[id] = rand;
326 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
327 if ( rand <= theProbability1->at(indexP1) )
break;
329 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
330 if ( rand <= theProbability2->at(indexP2) )
break;
336 xsela_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, ela1, ela2 ) * barn;
337 xscap_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, cap1, cap2 ) * barn;
343 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * barn;
345 }
else if ( lssf_flag == 1 ) {
347 std::vector< G4double >* theProbability1 =
theProbabilities->at(indexE - 1);
352 random_number_cache[id] = rand;
355 for ( indexP1 = 0; indexP1 < tableOrder; indexP1++ ) {
356 if ( rand <= theProbability1->at(indexP1) )
break;
358 for ( indexP2 = 0; indexP2 < tableOrder; indexP2++ ) {
359 if ( rand <= theProbability2->at(indexP2) )
break;
363 for (
G4int j = 0; j < (
G4int)ele->GetNumberOfIsotopes(); j++ ) {
364 if (
A == (
G4int)ele->GetIsotope(j)->GetN() ) {
369 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ];
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;
392 G4double fissionXS = weightedfissionXS / frac;
395 xsfiss_cache[id] = theInt.Lin(kineticEnergy, ene1, ene2, fiss1, fiss2) * fissionXS;
398 G4double fissionXS = weightedfissionXS / frac;
401 xsfiss_cache[id] = theInt.Lin( kineticEnergy, ene1, ene2, fiss1, fiss2 ) * fissionXS;
405 if ( MTnumber == 2 ) {
407 }
else if ( MTnumber == 102 ) {
409 }
else if ( MTnumber == 18 ) {
412 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)