Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 //-------------------------------------------- 28 // 29 // Geant4 source file 30 // 31 // File name: G4ParticleHPIsoProbabilityT 32 // 33 // Authors: Marek Zmeskal (CTU, Czech Tec 34 // Loic Thulliez (CEA France) 35 // 36 // Creation date: 4 June 2024 37 // 38 // Description: Class for the probability 39 // and for the given tempera 40 // It reads the files with p 41 // finds the correct cross-s 42 // 43 // Modifications: 44 // 45 // ------------------------------------------- 46 // 47 // 48 49 #include "G4ParticleHPIsoProbabilityTable_NJOY 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::G4Partic 67 tableOrder = 0; 68 lssf_flag = -1; 69 } 70 71 ///------------------------------------------- 72 G4ParticleHPIsoProbabilityTable_NJOY::~G4Parti 73 74 ///------------------------------------------- 75 void G4ParticleHPIsoProbabilityTable_NJOY::Ini 76 Z = theZ; 77 A = theA; 78 m = them; 79 T = theT; 80 G4cout << "The NJOY probability tables are b 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 + filena 89 std::istringstream theData( std::ios::in ); 90 G4ParticleHPManager::GetInstance()->GetDataS 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 99 theEnergies = new G4ParticleHPVector(nEner 100 theProbabilities = new std::vector< std::v 101 theElasticData = new std::vector< std::vec 102 theCaptureData = new std::vector< std::vec 103 theFissionData = new std::vector< std::vec 104 G4double probability, total, elastic, capt 105 for ( G4int i = 0; i < nEnergies; i++ ) { 106 theData >> energy; 107 theEnergies->SetEnergy( i, energy * eV ) 108 std::vector< G4double >* vecprob = new s 109 std::vector< G4double >* vecela = new s 110 std::vector< G4double >* veccap = new s 111 std::vector< G4double >* vecfiss = new s 112 for ( G4int j = 0; j < tableOrder; j++ ) 113 theData >> probability >> total >> ela 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 su 125 } else { 126 G4cout << "No probability tables found for 127 } 128 } 129 130 ///------------------------------------------- 131 G4double G4ParticleHPIsoProbabilityTable_NJOY: 132 const G4Element* ele, G4double& kineticEnerg 133 if ( kineticEnergy == energy_cache[id] ) { 134 if ( MTnumber == 2 ) { // elasti 135 return xsela_cache[id]; 136 } else if ( MTnumber == 102 ) { // radiat 137 return xscap_cache[id]; 138 } else if ( MTnumber == 18 ) { // fissio 139 return xsfiss_cache[id]; 140 } 141 } 142 energy_cache[id] = kineticEnergy; 143 if ( kineticEnergy < Emin || kineticEnergy 144 // if the kinetic energy is outside of the 145 G4int indexEl = (G4int)ele->GetIndex(); 146 G4int isotopeJ = -1; // index of isotope 147 for ( G4int j = 0; j < (G4int)ele->GetNumb 148 if ( A == (G4int)ele->GetIsotope(j)->Get 149 isotopeJ = j; 150 break; 151 } 152 } 153 G4double frac = ele->GetRelativeAbundanceV 154 G4double weightedelasticXS; 155 G4double weightedcaptureXS; 156 if ( G4ParticleHPManager::GetInstance()->G 157 weightedelasticXS = (*G4ParticleHPManage 158 weightedcaptureXS = (*G4ParticleHPManage 159 } else { 160 weightedelasticXS = this->GetDopplerBroa 161 weightedcaptureXS = this->GetDopplerBroa 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()- 169 G4double weightedfissionXS = (*G4Parti 170 xsfiss_cache[id] = weightedfissionXS / 171 } else { 172 G4double weightedfissionXS = this->Get 173 xsfiss_cache[id] = weightedfissionXS / 174 } 175 } 176 } else if ( lssf_flag == 0 ) { 177 G4int indexE = theEnergies->GetEnergyIndex 178 std::vector< G4double >* theProbability1 = 179 std::vector< G4double >* theProbability2 = 180 G4double ene1 = theEnergies->GetEnergy(ind 181 G4double ene2 = theEnergies->GetEnergy(ind 182 G4double rand = random_number; 183 G4int indexP1; 184 G4int indexP2; 185 for ( indexP1 = 0; indexP1 < tableOrder; i 186 if ( rand <= theProbability1->at(indexP1 187 } 188 for ( indexP2 = 0; indexP2 < tableOrder; i 189 if ( rand <= theProbability2->at(indexP2 190 } 191 G4double ela1 = theElasticData->at(indexE 192 G4double ela2 = theElasticData->at(indexE) 193 G4double cap1 = theCaptureData->at(indexE 194 G4double cap2 = theCaptureData->at(indexE) 195 xsela_cache[id] = theInt.Lin( kineticEnerg 196 xscap_cache[id] = theInt.Lin( kineticEnerg 197 if ( Z < 88 ) { 198 xsfiss_cache[id] = 0.0; 199 } else { 200 G4double fiss1 = theFissionData->at(inde 201 G4double fiss2 = theFissionData->at(inde 202 xsfiss_cache[id] = theInt.Lin( kineticEn 203 } 204 } else if ( lssf_flag == 1 ) { 205 G4int indexE = theEnergies->GetEnergyIndex 206 std::vector< G4double >* theProbability1 = 207 std::vector< G4double >* theProbability2 = 208 G4double ene1 = theEnergies->GetEnergy(ind 209 G4double ene2 = theEnergies->GetEnergy(ind 210 G4double rand = random_number; 211 G4int indexP1; 212 G4int indexP2; 213 for ( indexP1 = 0; indexP1 < tableOrder; i 214 if ( rand <= theProbability1->at(indexP1 215 } 216 for ( indexP2 = 0; indexP2 < tableOrder; i 217 if ( rand <= theProbability2->at(indexP2 218 } 219 G4int indexEl = (G4int)ele->GetIndex(); 220 G4int isotopeJ = -1; // index of isotope 221 for ( G4int j = 0; j < (G4int)ele->GetNumb 222 if ( A == (G4int)ele->GetIsotope(j)->Get 223 isotopeJ = j; 224 break; 225 } 226 } 227 G4double frac = ele->GetRelativeAbundanceV 228 G4double weightedelasticXS; 229 G4double weightedcaptureXS; 230 if ( G4ParticleHPManager::GetInstance()->G 231 weightedelasticXS = (*G4ParticleHPManage 232 weightedcaptureXS = (*G4ParticleHPManage 233 } else { 234 weightedelasticXS = this->GetDopplerBroa 235 weightedcaptureXS = this->GetDopplerBroa 236 } 237 G4double ela1 = theElasticData->at(indexE 238 G4double ela2 = theElasticData->at(indexE) 239 G4double cap1 = theCaptureData->at(indexE 240 G4double cap2 = theCaptureData->at(indexE) 241 G4double elasticXS = weightedelasticXS / f 242 G4double captureXS = weightedcaptureXS / f 243 xsela_cache[id] = theInt.Lin( kineticEnerg 244 xscap_cache[id] = theInt.Lin( kineticEnerg 245 if ( Z < 88 ) { 246 xsfiss_cache[id] = 0.0; 247 } else { 248 if ( G4ParticleHPManager::GetInstance()- 249 G4double weightedfissionXS = (*G4Parti 250 G4double fissionXS = weightedfissionXS 251 G4double fiss1 = theFissionData->at(in 252 G4double fiss2 = theFissionData->at(in 253 xsfiss_cache[id] = theInt.Lin( kinetic 254 } else { 255 G4double weightedfissionXS = this->Get 256 G4double fissionXS = weightedfissionXS 257 G4double fiss1 = theFissionData->at(in 258 G4double fiss2 = theFissionData->at(in 259 xsfiss_cache[id] = theInt.Lin( kinetic 260 } 261 } 262 } 263 if ( MTnumber == 2 ) { // elastic 264 return xsela_cache[id]; 265 } else if ( MTnumber == 102 ) { // radiativ 266 return xscap_cache[id]; 267 } else if ( MTnumber == 18 ) { // fission 268 return xsfiss_cache[id]; 269 } else { 270 G4cout << "Reaction was not found, returns 271 return 0; 272 } 273 } 274 275 ///------------------------------------------- 276 G4double G4ParticleHPIsoProbabilityTable_NJOY: 277 const G4Element* ele, G4double& kineticEnerg 278 energy_cache[id] = kineticEnergy; 279 if ( kineticEnergy < Emin || kineticEnergy 280 // if the kinetic energy is outside of the 281 G4int indexEl = (G4int)ele->GetIndex(); 282 G4int isotopeJ = -1; // index of isotope 283 for ( G4int j = 0; j < (G4int)ele->GetNumb 284 if ( A == (G4int)ele->GetIsotope(j)->Get 285 isotopeJ = j; 286 break; 287 } 288 } 289 G4double frac = ele->GetRelativeAbundanceV 290 G4double weightedelasticXS; 291 G4double weightedcaptureXS; 292 if ( G4ParticleHPManager::GetInstance()->G 293 weightedelasticXS = (*G4ParticleHPManage 294 weightedcaptureXS = (*G4ParticleHPManage 295 } else { 296 weightedelasticXS = this->GetDopplerBroa 297 weightedcaptureXS = this->GetDopplerBroa 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()- 305 G4double weightedfissionXS = (*G4Parti 306 xsfiss_cache[id] = weightedfissionXS / 307 } else { 308 G4double weightedfissionXS = this->Get 309 xsfiss_cache[id] = weightedfissionXS / 310 } 311 } 312 } else if ( lssf_flag == 0 ) { 313 G4int indexE = theEnergies->GetEnergyIndex 314 std::vector< G4double >* theProbability1 = 315 std::vector< G4double >* theProbability2 = 316 G4double ene1 = theEnergies->GetEnergy(ind 317 G4double ene2 = theEnergies->GetEnergy(ind 318 G4double rand = G4UniformRand(); 319 random_number_cache[id] = rand; 320 G4int indexP1; 321 G4int indexP2; 322 for ( indexP1 = 0; indexP1 < tableOrder; i 323 if ( rand <= theProbability1->at(indexP1 324 } 325 for ( indexP2 = 0; indexP2 < tableOrder; i 326 if ( rand <= theProbability2->at(indexP2 327 } 328 G4double ela1 = theElasticData->at(indexE 329 G4double ela2 = theElasticData->at(indexE) 330 G4double cap1 = theCaptureData->at(indexE 331 G4double cap2 = theCaptureData->at(indexE) 332 xsela_cache[id] = theInt.Lin( kineticEnerg 333 xscap_cache[id] = theInt.Lin( kineticEnerg 334 if ( Z < 88 ) { 335 xsfiss_cache[id] = 0.0; 336 } else { 337 G4double fiss1 = theFissionData->at(inde 338 G4double fiss2 = theFissionData->at(inde 339 xsfiss_cache[id] = theInt.Lin( kineticEn 340 } 341 } else if ( lssf_flag == 1 ) { 342 G4int indexE = theEnergies->GetEnergyIndex 343 std::vector< G4double >* theProbability1 = 344 std::vector< G4double >* theProbability2 = 345 G4double ene1 = theEnergies->GetEnergy(ind 346 G4double ene2 = theEnergies->GetEnergy(ind 347 G4double rand = G4UniformRand(); 348 random_number_cache[id] = rand; 349 G4int indexP1; 350 G4int indexP2; 351 for ( indexP1 = 0; indexP1 < tableOrder; i 352 if ( rand <= theProbability1->at(indexP1 353 } 354 for ( indexP2 = 0; indexP2 < tableOrder; i 355 if ( rand <= theProbability2->at(indexP2 356 } 357 G4int indexEl = (G4int)ele->GetIndex(); 358 G4int isotopeJ = -1; // index of isotope 359 for ( G4int j = 0; j < (G4int)ele->GetNumb 360 if ( A == (G4int)ele->GetIsotope(j)->Get 361 isotopeJ = j; 362 break; 363 } 364 } 365 G4double frac = ele->GetRelativeAbundanceV 366 G4double weightedelasticXS; 367 G4double weightedcaptureXS; 368 if ( G4ParticleHPManager::GetInstance()->G 369 weightedelasticXS = (*G4ParticleHPManage 370 weightedcaptureXS = (*G4ParticleHPManage 371 } else { 372 weightedelasticXS = this->GetDopplerBroa 373 weightedcaptureXS = this->GetDopplerBroa 374 } 375 G4double ela1 = theElasticData->at(indexE 376 G4double ela2 = theElasticData->at(indexE) 377 G4double cap1 = theCaptureData->at(indexE 378 G4double cap2 = theCaptureData->at(indexE) 379 G4double elasticXS = weightedelasticXS / f 380 G4double captureXS = weightedcaptureXS / f 381 xsela_cache[id] = theInt.Lin( kineticEnerg 382 xscap_cache[id] = theInt.Lin( kineticEnerg 383 if ( Z < 88 ) { 384 xsfiss_cache[id] = 0.0; 385 } else { 386 if ( G4ParticleHPManager::GetInstance()- 387 G4double weightedfissionXS = (*G4Par 388 G4double fissionXS = weightedfission 389 G4double fiss1 = theFissionData->at( 390 G4double fiss2 = theFissionData->at( 391 xsfiss_cache[id] = theInt.Lin(kineti 392 } else { 393 G4double weightedfissionXS = this->G 394 G4double fissionXS = weightedfission 395 G4double fiss1 = theFissionData->at( 396 G4double fiss2 = theFissionData->at( 397 xsfiss_cache[id] = theInt.Lin( kinet 398 } 399 } 400 } 401 if ( MTnumber == 2 ) { // elastic 402 return xsela_cache[id]; 403 } else if ( MTnumber == 102 ) { // radiativ 404 return xscap_cache[id]; 405 } else if ( MTnumber == 18 ) { // fission 406 return xsfiss_cache[id]; 407 } else { 408 G4cout << "Reaction was not found, returns 409 return 0; 410 } 411 } 412