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 // particle_hp -- source file 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron trans 29 // 30 // 080801 Give a warning message for irregular 31 // Introduce theNDLDataA,Z which has A 32 // 081024 G4NucleiPropertiesTable:: to G4Nucle 33 // 101111 Add Special treatment for Be9(n,2n)B 34 // 35 // P. Arce, June-2014 Conversion neutron_hp to 36 // 37 // June-2019 - E. Mendoza --> Added protection 38 // is not applied when data is in MF=6 format 39 // (add Q value info to G4ParticleHPNBodyPhase 40 41 #include "G4ParticleHPInelasticBaseFS.hh" 42 43 #include "G4Alpha.hh" 44 #include "G4Electron.hh" 45 #include "G4He3.hh" 46 #include "G4IonTable.hh" 47 #include "G4NucleiProperties.hh" 48 #include "G4Nucleus.hh" 49 #include "G4ParticleHPDataUsed.hh" 50 #include "G4ParticleHPManager.hh" 51 52 G4ParticleHPInelasticBaseFS::G4ParticleHPInela 53 { 54 hasXsec = true; 55 theXsection = new G4ParticleHPVector; 56 } 57 58 G4ParticleHPInelasticBaseFS::~G4ParticleHPInel 59 { 60 delete theXsection; 61 delete theEnergyDistribution; 62 delete theFinalStatePhotons; 63 delete theEnergyAngData; 64 delete theAngularDistribution; 65 } 66 67 void G4ParticleHPInelasticBaseFS::InitGammas(G 68 { 69 G4int Z = G4lrint(ZR); 70 G4int A = G4lrint(AR); 71 std::ostringstream ost; 72 ost << gammaPath << "z" << Z << ".a" << A; 73 G4String aName = ost.str(); 74 std::ifstream from(aName, std::ios::in); 75 76 if (!from) return; // no data found for thi 77 std::ifstream theGammaData(aName, std::ios:: 78 79 theNuclearMassDifference = G4NucleiPropertie 80 - G4NucleiProperties::GetBindingEnergy(the 81 theGammas.Init(theGammaData); 82 } 83 84 void G4ParticleHPInelasticBaseFS::Init(G4doubl 85 const G 86 const G 87 { 88 gammaPath = fManager->GetNeutronHPPath() + " 89 G4String tString = dirName; 90 SetA_Z(A, Z, M); 91 G4bool dbool = true; 92 G4ParticleHPDataUsed aFile = 93 theNames.GetName(theBaseA, theBaseZ, M, tS 94 SetAZMs(aFile); 95 G4String filename = aFile.GetName(); 96 #ifdef G4VERBOSE 97 if (fManager->GetDEBUG()) 98 G4cout << " G4ParticleHPInelasticBaseFS::I 99 #endif 100 101 if (!dbool || (theBaseZ <= 2 && (theNDLDataZ 102 { 103 #ifdef G4PHPDEBUG 104 if (fManager->GetDEBUG()) 105 G4cout << "Skipped = " << filename << " 106 #endif 107 hasAnyData = false; 108 hasFSData = false; 109 hasXsec = false; 110 return; 111 } 112 113 std::istringstream theData(std::ios::in); 114 fManager->GetDataStream(filename, theData); 115 116 if (!theData) //"!" is a operator of ios 117 { 118 hasAnyData = false; 119 hasFSData = false; 120 hasXsec = false; 121 return; // no data for exactly this isoto 122 } 123 // here we go 124 G4int infoType, dataType, dummy = INT_MAX; 125 hasFSData = false; 126 while (theData >> infoType) // Loop checkin 127 { 128 theData >> dataType; 129 130 if (dummy == INT_MAX) theData >> Qvalue >> 131 Qvalue *= CLHEP::eV; 132 // In G4NDL4.5 this value is the MT number 133 // in others is que Q-value in eV 134 135 if (dataType == 3) { 136 G4int total; 137 theData >> total; 138 theXsection->Init(theData, total, CLHEP: 139 } 140 else if (dataType == 4) { 141 theAngularDistribution = new G4ParticleH 142 theAngularDistribution->Init(theData); 143 hasFSData = true; 144 } 145 else if (dataType == 5) { 146 theEnergyDistribution = new G4ParticleHP 147 theEnergyDistribution->Init(theData); 148 hasFSData = true; 149 } 150 else if (dataType == 6) { 151 theEnergyAngData = new G4ParticleHPEnAng 152 theEnergyAngData->Init(theData); 153 hasFSData = true; 154 } 155 else if (dataType == 12) { 156 theFinalStatePhotons = new G4ParticleHPP 157 theFinalStatePhotons->InitMean(theData); 158 hasFSData = true; 159 } 160 else if (dataType == 13) { 161 theFinalStatePhotons = new G4ParticleHPP 162 theFinalStatePhotons->InitPartials(theDa 163 hasFSData = true; 164 } 165 else if (dataType == 14) { 166 theFinalStatePhotons->InitAngular(theDat 167 hasFSData = true; 168 } 169 else if (dataType == 15) { 170 theFinalStatePhotons->InitEnergies(theDa 171 hasFSData = true; 172 } 173 else { 174 G4ExceptionDescription ed; 175 ed << "Z=" << theBaseZ << " A=" << theBa 176 << " projectile: " << theProjectile->GetPar 177 G4Exception("G4ParticleHPInelasticBaseFS 178 ed, "Data-type unknown"); 179 } 180 } 181 } 182 183 void G4ParticleHPInelasticBaseFS::BaseApply(co 184 G4 185 G4 186 { 187 // G4cout << "G4ParticleHPInelasticBaseFS::B 188 // prepare neutron 189 if (theResult.Get() == nullptr) theResult.Pu 190 theResult.Get()->Clear(); 191 G4double eKinetic = theTrack.GetKineticEnerg 192 G4ReactionProduct incidReactionProduct(theTr 193 incidReactionProduct.SetMomentum(theTrack.Ge 194 incidReactionProduct.SetKineticEnergy(eKinet 195 196 // prepare target 197 G4double targetMass = G4NucleiProperties::Ge 198 /CLHEP::neutron_mass_c2; 199 /* 200 G4cout << " Target Z=" << theBaseZ << " A=" 201 << " projectile: " << theTrack.GetDefinitio 202 << " Ekin(MeV)=" << theTrack.GetKineticEner 203 << " Mtar/Mn=" << targetMass 204 << " hasFS=" << HasFSData() << G4endl; 205 */ 206 // give priority to ENDF vales for target ma 207 if (theEnergyAngData != nullptr) { 208 G4double x = theEnergyAngData->GetTargetMa 209 if(0.0 < x) { targetMass = x/CLHEP::neutro 210 } 211 if (theAngularDistribution != nullptr) { 212 G4double x = theAngularDistribution->GetTa 213 if(0.0 < x) { targetMass = x/CLHEP::neutr 214 } 215 216 // 110512 TKDB ENDF-VII.0 21Sc45 has trouble 217 // recorded. TKDB targetMass = 0; ENDF-VII.0 218 219 G4Nucleus aNucleus; 220 G4ReactionProduct theTarget; 221 222 G4ThreeVector neuVelo = 223 incidReactionProduct.GetMomentum() / CLHEP 224 225 theTarget = 226 aNucleus.GetBiasedThermalNucleus(targetMas 227 theTrack.GetMaterial()->GetTemper 228 theTarget.SetDefinition(ionTable->GetIon(the 229 230 // prepare energy in target rest frame 231 G4ReactionProduct boosted; 232 boosted.Lorentz(incidReactionProduct, theTar 233 eKinetic = boosted.GetKineticEnergy(); 234 235 G4double orgMomentum = boosted.GetMomentum() 236 237 // Take N-body phase-space distribution, if 238 if (!HasFSData()) // adding the residual is 239 { 240 G4ParticleHPNBodyPhaseSpace thePhaseSpaceD 241 G4double aPhaseMass = 0; 242 G4int ii; 243 for (ii = 0; ii < nDef; ++ii) { 244 aPhaseMass += theDefs[ii]->GetPDGMass(); 245 } 246 247 //---------------------------------------- 248 if (std::abs(Qvalue) < CLHEP::keV) { 249 // Not in the G4NDL lib or not calculate 250 251 // Calculate residual: 252 G4int ResidualA = theBaseA; 253 G4int ResidualZ = theBaseZ; 254 for (ii = 0; ii < nDef; ++ii) { 255 ResidualZ -= theDefs[ii]->GetAtomicNum 256 ResidualA -= theDefs[ii]->GetBaryonNum 257 } 258 259 if (ResidualA > 0 && ResidualZ > 0) { 260 G4ParticleDefinition* resid = ionTable 261 Qvalue = 262 incidReactionProduct.GetMass() + the 263 } 264 265 if (std::abs(Qvalue) > 400 * CLHEP::MeV) 266 // Then Q value is probably too large 267 Qvalue = 1.1 * CLHEP::keV; 268 } 269 } 270 //---------------------------------------- 271 272 thePhaseSpaceDistribution.Init(aPhaseMass, 273 thePhaseSpaceDistribution.SetProjectileRP( 274 thePhaseSpaceDistribution.SetTarget(&theTa 275 thePhaseSpaceDistribution.SetQValue(Qvalue 276 277 for (ii = 0; ii < nDef; ++ii) { 278 G4double massCode = 1000. * std::abs(the 279 massCode += theDefs[ii]->GetBaryonNumber 280 G4double dummy = 0; 281 G4ReactionProduct* aSec = thePhaseSpaceD 282 aSec->Lorentz(*aSec, -1. * theTarget); 283 auto aPart = new G4DynamicParticle(); 284 aPart->SetDefinition(aSec->GetDefinition 285 aPart->SetMomentum(aSec->GetMomentum()); 286 delete aSec; 287 theResult.Get()->AddSecondary(aPart, sec 288 #ifdef G4VERBOSE 289 if (fManager->GetDEBUG()) 290 G4cout << this << " G4ParticleHPInelas 291 << aPart->GetParticleDefinition 292 << " E= " << aPart->GetKineticE 293 << theResult.Get()->GetNumberOf 294 #endif 295 } 296 297 theResult.Get()->SetStatusChange(stopAndKi 298 // Final momentum check should be done bef 299 const G4ParticleDefinition* targ_pd = ionT 300 G4double mass = targ_pd->GetPDGMass(); 301 G4LorentzVector targ_4p_lab(theTarget.GetM 302 std::sqrt(mass*mass + theTarget.GetMo 303 G4LorentzVector proj_4p_lab = theTrack.Get 304 G4LorentzVector init_4p_lab = proj_4p_lab 305 adjust_final_state(init_4p_lab); 306 return; 307 } 308 309 // set target and neutron in the relevant ex 310 if (theAngularDistribution != nullptr) { 311 theAngularDistribution->SetTarget(theTarge 312 theAngularDistribution->SetProjectileRP(in 313 } 314 else if (theEnergyAngData != nullptr) { 315 theEnergyAngData->SetTarget(theTarget); 316 theEnergyAngData->SetProjectileRP(incidRea 317 } 318 319 G4ReactionProductVector* tmpHadrons = nullpt 320 G4int dummy; 321 std::size_t i; 322 323 if (theEnergyAngData != nullptr) { 324 tmpHadrons = theEnergyAngData->Sample(eKin 325 auto hproj = theTrack.GetDefinition(); 326 G4int projA = hproj->GetBaryonNumber(); 327 G4int projZ = G4lrint(hproj->GetPDGCharge( 328 329 if (!fManager->GetDoNotAdjustFinalState()) 330 // Adjust A and Z in the case of miss mu 331 if (tmpHadrons != nullptr) { 332 G4int sumA = 0; 333 G4int sumZ = 0; 334 G4int maxA = 0; 335 G4int jAtMaxA = 0; 336 for (G4int j = 0; j != (G4int)tmpHadro 337 auto had = ((*tmpHadrons)[j])->GetDefiniti 338 if (had->GetBaryonNumber() > maxA) { 339 maxA = had->GetBaryonNumber(); 340 jAtMaxA = j; 341 } 342 sumA += had->GetBaryonNumber(); 343 sumZ += G4lrint(had->GetPDGCharge()/ 344 } 345 G4int dA = theBaseA + projA - sumA; 346 G4int dZ = theBaseZ + projZ - sumZ; 347 if (dA < 0 || dZ < 0) { 348 auto p = ((*tmpHadrons)[jAtMaxA])->G 349 G4int newA = p->GetBaryonNumber() + 350 G4int newZ = G4lrint(p->GetPDGCharge 351 if (newA > newZ && newZ > 0) { 352 G4ParticleDefinition* pd = ionTabl 353 ((*tmpHadrons)[jAtMaxA])->SetDefin 354 } 355 } 356 } 357 } 358 } 359 else if (theAngularDistribution != nullptr) 360 361 auto Done = new G4bool[nDef]; 362 G4int i0; 363 for (i0 = 0; i0 < nDef; ++i0) 364 Done[i0] = false; 365 366 tmpHadrons = new G4ReactionProductVector; 367 G4ReactionProduct* aHadron; 368 G4double localMass = G4NucleiProperties::G 369 G4ThreeVector bufferedDirection(0, 0, 0); 370 for (i0 = 0; i0 < nDef; ++i0) { 371 if (!Done[i0]) { 372 aHadron = new G4ReactionProduct; 373 if (theEnergyDistribution != nullptr) 374 aHadron->SetDefinition(theDefs[i0]); 375 aHadron->SetKineticEnergy(theEnergyD 376 } 377 else if (nDef == 1) { 378 aHadron->SetDefinition(theDefs[i0]); 379 aHadron->SetKineticEnergy(eKinetic); 380 } 381 else if (nDef == 2) { 382 aHadron->SetDefinition(theDefs[i0]); 383 aHadron->SetKineticEnergy(50 * CLHEP 384 } 385 else { 386 throw G4HadronicException( 387 __FILE__, __LINE__, 388 "No energy distribution to sample 389 } 390 theAngularDistribution->SampleAndUpdat 391 if (theEnergyDistribution == nullptr & 392 if (i0 == 0) { 393 G4double mass1 = theDefs[0]->GetPD 394 G4double mass2 = theDefs[1]->GetPD 395 G4double massn = theProjectile->Ge 396 G4int z1 = theBaseZ - 397 G4lrint((theDefs[0]->GetPDGCharge() + 398 G4int a1 = theBaseA - theDefs[0]-> 399 G4double concreteMass = G4NucleiPr 400 G4double availableEnergy = eKineti 401 - concreteMass; 402 // available kinetic energy in CMS 403 G4double emin = availableEnergy + 404 - std::sqrt((mass1 + mass2) * (m 405 G4double p1 = std::sqrt(2. * mass2 406 bufferedDirection = p1 * aHadron-> 407 #ifdef G4PHPDEBUG 408 if (fManager->GetDEBUG()) // @@@@ 409 { 410 G4cout << "G4ParticleHPInelastic 411 << " " << a1 << " " << theBaseA << " 412 << " " << emin << G4endl; 413 } 414 #endif 415 } 416 else { 417 bufferedDirection = -bufferedDirec 418 } 419 // boost from cms to lab 420 aHadron->SetTotalEnergy( 421 std::sqrt(aHadron->GetMass() * aHa 422 aHadron->SetMomentum(bufferedDirecti 423 aHadron->Lorentz(*aHadron, -1. * (th 424 #ifdef G4PHPDEBUG 425 if (fManager->GetDEBUG()) { 426 G4cout << " G4ParticleHPInelasticB 427 << " P=" << aHadron->GetMomentum() 428 << " bufDir^2=" << bufferedDirection.ma 429 << G4endl; 430 } 431 #endif 432 } 433 tmpHadrons->push_back(aHadron); 434 #ifdef G4PHPDEBUG 435 if (fManager->GetDEBUG()) 436 G4cout << " G4ParticleHPInelasticBas 437 << aHadron->GetDefinition()-> 438 << " E= " << aHadron->GetKine 439 #endif 440 } 441 } 442 delete[] Done; 443 } 444 else { 445 throw G4HadronicException(__FILE__, __LINE 446 } 447 448 G4ReactionProductVector* thePhotons = nullpt 449 if (theFinalStatePhotons != nullptr) { 450 // the photon distributions are in the Nuc 451 G4ReactionProduct boosted_tmp; 452 boosted_tmp.Lorentz(incidReactionProduct, 453 G4double anEnergy = boosted_tmp.GetKinetic 454 thePhotons = theFinalStatePhotons->GetPhot 455 if (thePhotons != nullptr) { 456 for (i = 0; i < thePhotons->size(); ++i) 457 // back to lab 458 thePhotons->operator[](i)->Lorentz(*(t 459 } 460 } 461 } 462 else if (theEnergyAngData != nullptr) { 463 // PA130927: do not create photons to adju 464 G4bool bAdjustPhotons{true}; 465 #ifdef PHP_AS_HP 466 bAdjustPhotons = true; 467 #else 468 if (fManager->GetDoNotAdjustFinalState()) 469 #endif 470 471 if (bAdjustPhotons) { 472 G4double theGammaEnergy = theEnergyAngDa 473 G4double anEnergy = boosted.GetKineticEn 474 theGammaEnergy = anEnergy - theGammaEner 475 theGammaEnergy += theNuclearMassDifferen 476 G4double eBindProducts = 0; 477 G4double eBindN = 0; 478 G4double eBindP = 0; 479 G4double eBindD = G4NucleiProperties::Ge 480 G4double eBindT = G4NucleiProperties::Ge 481 G4double eBindHe3 = G4NucleiProperties:: 482 G4double eBindA = G4NucleiProperties::Ge 483 G4int ia = 0; 484 for (auto const & had : *tmpHadrons) { 485 if (had->GetDefinition() == G4Neutron: 486 eBindProducts += eBindN; 487 } 488 else if (had->GetDefinition() == G4Pro 489 eBindProducts += eBindP; 490 } 491 else if (had->GetDefinition() == G4Deu 492 eBindProducts += eBindD; 493 } 494 else if (had->GetDefinition() == G4Tri 495 eBindProducts += eBindT; 496 } 497 else if (had->GetDefinition() == G4He3 498 eBindProducts += eBindHe3; 499 } 500 else if (had->GetDefinition() == G4Alp 501 eBindProducts += eBindA; 502 ia++; 503 } 504 } 505 theGammaEnergy += eBindProducts; 506 507 #ifdef G4PHPDEBUG 508 if (fManager->GetDEBUG()) 509 G4cout << " G4ParticleHPInelasticBaseF 510 << " eBindProducts " << eBindPr 511 #endif 512 513 // Special treatment for Be9 + n -> 2n + 514 if (theBaseZ == 4 && theBaseA == 9) { 515 // This only valid for G4NDL3.13,,, 516 if (ia == 2 && std::abs(G4NucleiProper 517 G4NucleiProperties::GetBindingEnergy(9 518 theNuclearMassDifference) < CLHEP::keV 519 { 520 theGammaEnergy -= (2 * eBindA); 521 } 522 } 523 524 if (theGammaEnergy > 0.0) { 525 for (G4int iLevel = theGammas.GetNumberOfLev 526 G4double e = theGammas.GetLevelEnerg 527 if (e < theGammaEnergy) { 528 thePhotons = theGammas.GetDecayGam 529 theGammaEnergy -= e; 530 break; 531 } 532 } 533 } 534 } 535 } 536 // fill the result 537 std::size_t nSecondaries = tmpHadrons->size( 538 std::size_t nPhotons = 0; 539 if (thePhotons != nullptr) { 540 nPhotons = thePhotons->size(); 541 } 542 nSecondaries += nPhotons; 543 544 G4DynamicParticle* theSec; 545 #ifdef G4PHPDEBUG 546 if (fManager->GetDEBUG()) 547 G4cout << " G4ParticleHPInelasticBaseFS::B 548 << G4endl; 549 #endif 550 551 for (i = 0; i < nSecondaries - nPhotons; ++i 552 theSec = new G4DynamicParticle; 553 theSec->SetDefinition(tmpHadrons->operator 554 theSec->SetMomentum(tmpHadrons->operator[] 555 theResult.Get()->AddSecondary(theSec, secI 556 #ifdef G4PHPDEBUG 557 if (fManager->GetDEBUG()) 558 G4cout << " G4ParticleHPInelasticBaseFS: 559 << theSec->GetParticleDefinition( 560 << " E=" << theSec->GetKineticEne 561 << theResult.Get()->GetNumberOfSe 562 #endif 563 delete (*tmpHadrons)[i]; 564 } 565 #ifdef G4PHPDEBUG 566 if (fManager->GetDEBUG()) 567 G4cout << " G4ParticleHPInelasticBaseFS::B 568 #endif 569 if (thePhotons != nullptr) { 570 for (i = 0; i < nPhotons; ++i) { 571 theSec = new G4DynamicParticle; 572 theSec->SetDefinition((*thePhotons)[i]-> 573 theSec->SetMomentum((*thePhotons)[i]->Ge 574 theResult.Get()->AddSecondary(theSec, se 575 #ifdef G4PHPDEBUG 576 if (fManager->GetDEBUG()) 577 G4cout << " G4ParticleHPInelasticBaseF 578 << theSec->GetParticleDefinitio 579 << " E=" << theSec->GetKineticE 580 << theResult.Get()->GetNumberOf 581 #endif 582 delete thePhotons->operator[](i); 583 } 584 } 585 586 // some garbage collection 587 delete thePhotons; 588 delete tmpHadrons; 589 590 G4ParticleDefinition* targ_pd = ionTable->Ge 591 G4LorentzVector targ_4p_lab( 592 theTarget.GetMomentum(), 593 std::sqrt(targ_pd->GetPDGMass() * targ_pd- 594 G4LorentzVector proj_4p_lab = theTrack.Get4M 595 G4LorentzVector init_4p_lab = proj_4p_lab + 596 597 // if data in MF=6 format (no correlated par 598 // severe errors: 599 if (theEnergyAngData == nullptr) { 600 adjust_final_state(init_4p_lab); 601 } 602 603 // clean up the primary neutron 604 theResult.Get()->SetStatusChange(stopAndKill 605 } 606