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_CALENDF.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 CALENDF. 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_CALENDF.hh" 50 #include "G4SystemOfUnits.hh" 51 #include "Randomize.hh" 52 #include "G4ParticleHPVector.hh" 53 #include "G4DynamicParticle.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 #include <string> 62 #include <sstream> 63 64 ///-------------------------------------------------------------------------------------- 65 G4ParticleHPIsoProbabilityTable_CALENDF::G4ParticleHPIsoProbabilityTable_CALENDF() : theInelasticData(nullptr) {} 66 67 ///-------------------------------------------------------------------------------------- 68 G4ParticleHPIsoProbabilityTable_CALENDF::~G4ParticleHPIsoProbabilityTable_CALENDF() { 69 for ( auto it = theInelasticData->cbegin(); it != theInelasticData->cend(); ++it ) { 70 delete* it; 71 } 72 delete theInelasticData; 73 } 74 75 ///-------------------------------------------------------------------------------------- 76 void G4ParticleHPIsoProbabilityTable_CALENDF::Init( G4int theZ, G4int theA, G4int them, G4double theT, const G4String& dirName ) { 77 Z = theZ; 78 A = theA; 79 m = them; 80 T = theT; 81 G4cout << "The CALENDF probability tables are being initialized for Z=" << Z << " A=" << A << " and T=" << T << " K." << G4endl; 82 filename = std::to_string(Z) + "_" + std::to_string(A); 83 if ( m != 0 ) filename += "_m" + std::to_string(m); 84 G4String fullPathFileName = dirName + filename + "." + std::to_string( (G4int)(T) ) + ".pt"; 85 std::istringstream theData( std::ios::in ); 86 G4ParticleHPManager::GetInstance()->GetDataStream( fullPathFileName, theData ); 87 if ( theData.good() ) { 88 G4double emin; 89 G4double emax; 90 theData >> emin >> emax; 91 Emin = emin * eV; 92 Emax = emax * eV; 93 theData >> nEnergies; 94 theEnergies = new G4ParticleHPVector( nEnergies ); 95 theProbabilities = new std::vector< std::vector< G4double >* >; 96 theElasticData = new std::vector< std::vector< G4double >* >; 97 theCaptureData = new std::vector< std::vector< G4double >* >; 98 theFissionData = new std::vector< std::vector< G4double >* >; 99 theInelasticData = new std::vector< std::vector< G4double >* >; 100 G4double tableOrder; 101 G4double probability, total, elastic, capture, fission, inelastic; 102 for ( G4int i = 0; i < nEnergies; i++ ) { 103 theData >> emin >> emax >> tableOrder; 104 theEnergies->SetData( i, emax * eV, tableOrder ); 105 std::vector< G4double >* vecprob = new std::vector< G4double >; 106 std::vector< G4double >* vecela = new std::vector< G4double >; 107 std::vector< G4double >* veccap = new std::vector< G4double >; 108 std::vector< G4double >* vecfiss = new std::vector< G4double >; 109 std::vector< G4double >* vecinela = new std::vector< G4double >; 110 for ( G4int j = 0; j < tableOrder; j++ ) { 111 theData >> probability >> total >> elastic >> capture >> fission >> inelastic; 112 vecprob->push_back( probability ); 113 vecela->push_back( elastic * barn ); 114 veccap->push_back( capture * barn ); 115 vecfiss->push_back( fission * barn ); 116 vecinela->push_back( inelastic * barn ); 117 } 118 theProbabilities->push_back( vecprob ); 119 theElasticData->push_back( vecela ); 120 theCaptureData->push_back( veccap ); 121 theFissionData->push_back( vecfiss ); 122 theInelasticData->push_back( vecinela ); 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_CALENDF::GetCorrelatedIsoCrossSectionPT( const G4DynamicParticle* dp, 132 G4int MTnumber, const G4Element* ele, G4double &kineticEnergy, G4double &random_number, std::thread::id &id ) { 133 // if this function is called again for different reaction or the particle came to different region 134 // with the same temperature and unchanged temperature 135 if ( kineticEnergy == energy_cache[id] ) { 136 if ( MTnumber == 2 ) { // elastic cross section 137 return xsela_cache[id]; 138 } else if ( MTnumber == 102 ) { // radiative capture cross section 139 return xscap_cache[id]; 140 } else if ( MTnumber == 18 ) { // fission cross section 141 return xsfiss_cache[id]; 142 } else if ( MTnumber == 3 ) { // inelastic cross section 143 return xsinela_cache[id]; 144 } 145 } 146 energy_cache[id] = kineticEnergy; 147 if ( kineticEnergy < Emin || kineticEnergy > Emax ) { 148 // if the kinetic energy is outside of the URR limits for the given isotope, it finds the smooth cross section 149 G4int indexEl = (G4int)ele->GetIndex(); 150 G4int isotopeJ = -1; // index of isotope in the given element 151 G4int n_isotopes = (G4int)ele->GetNumberOfIsotopes(); 152 for ( G4int j = 0; j < n_isotopes; j++ ) { 153 if ( A == (G4int)( ele->GetIsotope(j)->GetN() ) ) { 154 isotopeJ = j; 155 break; 156 } 157 } 158 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ]; 159 G4double weightedelasticXS; 160 G4double weightedcaptureXS; 161 G4double weightedinelasticXS; 162 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) { 163 weightedelasticXS = (*G4ParticleHPManager::GetInstance()->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ ); 164 weightedcaptureXS = (*G4ParticleHPManager::GetInstance()->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ ); 165 weightedinelasticXS = ((*G4ParticleHPManager::GetInstance()->GetInelasticFinalStates(dp->GetDefinition()))[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ )) * barn; 166 } else { 167 weightedelasticXS = this->GetDopplerBroadenedElasticXS( dp, indexEl, isotopeJ ); 168 weightedcaptureXS = this->GetDopplerBroadenedCaptureXS( dp, indexEl, isotopeJ ); 169 weightedinelasticXS = this->GetDopplerBroadenedInelasticXS( dp, indexEl, isotopeJ ); 170 } 171 xsela_cache[id] = weightedelasticXS / frac; 172 xscap_cache[id] = weightedcaptureXS / frac; 173 xsinela_cache[id] = weightedinelasticXS / frac; 174 if ( Z < 88 ) { 175 xsfiss_cache[id] = 0.0; 176 } else { 177 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) { 178 G4double weightedfissionXS = (*G4ParticleHPManager::GetInstance()->GetFissionFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ ); 179 xsfiss_cache[id] = weightedfissionXS / frac; 180 } else { 181 G4double weightedfissionXS = this->GetDopplerBroadenedFissionXS( dp, indexEl, isotopeJ ); 182 xsfiss_cache[id] = weightedfissionXS / frac; 183 } 184 } 185 } else { 186 G4int indexE = theEnergies->GetEnergyIndex( kineticEnergy ); 187 G4int order = (G4int)( theEnergies->GetY(indexE) + 0.5 ); 188 std::vector< G4double >* theProbability = theProbabilities->at( indexE ); 189 G4double rand = random_number; 190 G4int indexP; 191 for ( indexP = 0; indexP < order; indexP++ ) { 192 if ( rand <= theProbability->at(indexP) ) break; 193 } 194 xsela_cache[id] = theElasticData->at(indexE)->at(indexP); 195 xscap_cache[id] = theCaptureData->at(indexE)->at(indexP); 196 xsinela_cache[id] = theInelasticData->at(indexE)->at(indexP); 197 if ( Z < 88 ) xsfiss_cache[id] = 0.0; 198 else xsfiss_cache[id] = theFissionData->at(indexE)->at(indexP); 199 } 200 if ( MTnumber == 2 ) { // elastic cross section 201 return xsela_cache[id]; 202 } else if ( MTnumber == 102 ) { // radiative capture cross section 203 return xscap_cache[id]; 204 } else if ( MTnumber == 18 ) { // fission cross section 205 return xsfiss_cache[id]; 206 } else if ( MTnumber == 3 ) { // inelastic cross section 207 return xsinela_cache[id]; 208 } else { 209 G4cout << "Reaction was not found, returns 0." << G4endl; 210 return 0; 211 } 212 } 213 214 ///-------------------------------------------------------------------------------------- 215 G4double G4ParticleHPIsoProbabilityTable_CALENDF::GetIsoCrossSectionPT( const G4DynamicParticle* dp, G4int MTnumber, 216 const G4Element* ele,G4double &kineticEnergy , std::map< std::thread::id, G4double > &random_number_cache, std::thread::id &id ) { 217 energy_cache[id] = kineticEnergy; 218 if ( kineticEnergy < Emin || kineticEnergy > Emax ) { 219 // if the kinetic energy is outside of the URR limits for the given isotope, it finds the smooth cross section 220 G4int indexEl = (G4int)ele->GetIndex(); 221 G4int isotopeJ = -1; //index of isotope in the given element 222 G4int n_isotopes = (G4int)ele->GetNumberOfIsotopes(); 223 for ( G4int j = 0; j < n_isotopes; j++ ) { 224 if ( A == (G4int)ele->GetIsotope(j)->GetN() ) { 225 isotopeJ = j; 226 break; 227 } 228 } 229 G4double frac = ele->GetRelativeAbundanceVector()[isotopeJ]; 230 G4double weightedelasticXS; 231 G4double weightedcaptureXS; 232 G4double weightedinelasticXS; 233 if (G4ParticleHPManager::GetInstance()->GetNeglectDoppler()) { 234 weightedelasticXS = (*G4ParticleHPManager::GetInstance()->GetElasticFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ ); 235 weightedcaptureXS = (*G4ParticleHPManager::GetInstance()->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ ); 236 weightedinelasticXS = ((*G4ParticleHPManager::GetInstance()->GetInelasticFinalStates(dp->GetDefinition()))[indexEl]->GetWeightedXsec( kineticEnergy, isotopeJ )) * barn; 237 } else { 238 weightedelasticXS = this->GetDopplerBroadenedElasticXS( dp, indexEl, isotopeJ ); 239 weightedcaptureXS = this->GetDopplerBroadenedCaptureXS( dp, indexEl, isotopeJ ); 240 weightedinelasticXS = this->GetDopplerBroadenedInelasticXS( dp, indexEl, isotopeJ ); 241 } 242 xsela_cache[id] = weightedelasticXS / frac; 243 xscap_cache[id] = weightedcaptureXS / frac; 244 xsinela_cache[id] = weightedinelasticXS / frac; 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 xsfiss_cache[id] = weightedfissionXS / frac; 251 } else { 252 G4double weightedfissionXS = this->GetDopplerBroadenedFissionXS( dp, indexEl, isotopeJ ); 253 xsfiss_cache[id] = weightedfissionXS / frac; 254 } 255 } 256 } else { 257 G4int indexE = theEnergies->GetEnergyIndex( kineticEnergy ); 258 G4int order = (G4int)( theEnergies->GetY(indexE) + 0.5 ); 259 std::vector< G4double >* theProbability = theProbabilities->at(indexE); 260 G4double rand = G4UniformRand(); 261 random_number_cache[id] = rand; 262 G4int indexP; 263 for ( indexP = 0; indexP < order; indexP++ ) { 264 if ( rand <= theProbability->at(indexP) ) break; 265 } 266 xsela_cache[id] = theElasticData->at(indexE)->at(indexP); 267 xscap_cache[id] = theCaptureData->at(indexE)->at(indexP); 268 xsinela_cache[id] = theInelasticData->at(indexE)->at(indexP); 269 if ( Z < 88 ) xsfiss_cache[id] = 0.0; 270 else xsfiss_cache[id] = theFissionData->at(indexE)->at(indexP); 271 } 272 if (MTnumber == 2) { // elastic cross section 273 return xsela_cache[id]; 274 } else if (MTnumber == 102) { // radiative capture cross section 275 return xscap_cache[id]; 276 } else if (MTnumber == 18) { // fission cross section 277 return xsfiss_cache[id]; 278 } else if (MTnumber == 3) { // inelastic cross section 279 return xsinela_cache[id]; 280 } else { 281 G4cout << "Reaction was not found, returns 0." << G4endl; 282 return 0; 283 } 284 } 285