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 // GEANT4 Class source file 29 // 30 // G4RadioactiveDecay 31 // 32 // Author: D.H. Wright (SLAC) 33 // Date: 29 August 2017 34 // 35 // Based on the code of F. Lei and P.R. Trusc 36 // 37 ////////////////////////////////////////////// 38 39 #include "G4RadioactiveDecay.hh" 40 #include "G4RadioactivationMessenger.hh" 41 42 #include "G4SystemOfUnits.hh" 43 #include "G4DynamicParticle.hh" 44 #include "G4DecayProducts.hh" 45 #include "G4DecayTable.hh" 46 #include "G4ParticleChangeForRadDecay.hh" 47 #include "G4ITDecay.hh" 48 #include "G4BetaDecayType.hh" 49 #include "G4BetaMinusDecay.hh" 50 #include "G4BetaPlusDecay.hh" 51 #include "G4ECDecay.hh" 52 #include "G4AlphaDecay.hh" 53 #include "G4TritonDecay.hh" 54 #include "G4ProtonDecay.hh" 55 #include "G4NeutronDecay.hh" 56 #include "G4SFDecay.hh" 57 #include "G4VDecayChannel.hh" 58 #include "G4NuclearDecay.hh" 59 #include "G4RadioactiveDecayMode.hh" 60 #include "G4Fragment.hh" 61 #include "G4Ions.hh" 62 #include "G4IonTable.hh" 63 #include "G4BetaDecayType.hh" 64 #include "Randomize.hh" 65 #include "G4LogicalVolumeStore.hh" 66 #include "G4NuclearLevelData.hh" 67 #include "G4DeexPrecoParameters.hh" 68 #include "G4LevelManager.hh" 69 #include "G4ThreeVector.hh" 70 #include "G4Electron.hh" 71 #include "G4Positron.hh" 72 #include "G4Neutron.hh" 73 #include "G4Gamma.hh" 74 #include "G4Alpha.hh" 75 #include "G4Triton.hh" 76 #include "G4Proton.hh" 77 78 #include "G4HadronicProcessType.hh" 79 #include "G4HadronicProcessStore.hh" 80 #include "G4HadronicException.hh" 81 #include "G4LossTableManager.hh" 82 #include "G4VAtomDeexcitation.hh" 83 #include "G4UAtomicDeexcitation.hh" 84 #include "G4PhotonEvaporation.hh" 85 86 #include <vector> 87 #include <sstream> 88 #include <algorithm> 89 #include <fstream> 90 91 using namespace CLHEP; 92 93 G4RadioactiveDecay::G4RadioactiveDecay(const G 94 const G 95 : G4VRadioactiveDecay(processName, timeThres 96 { 97 #ifdef G4VERBOSE 98 if (GetVerboseLevel() > 1) { 99 G4cout << "G4RadioactiveDecay constructor: 100 << G4endl; 101 } 102 #endif 103 104 theRadioactivationMessenger = new G4Radioact 105 106 // Apply default values. 107 NSourceBin = 1; 108 SBin[0] = 0.* s; 109 SBin[1] = 1.* s; // Convert to ns 110 SProfile[0] = 1.; 111 SProfile[1] = 0.; 112 NDecayBin = 1; 113 DBin[0] = 0. * s ; 114 DBin[1] = 1. * s; 115 DProfile[0] = 1.; 116 DProfile[1] = 0.; 117 decayWindows[0] = 0; 118 G4RadioactivityTable* rTable = new G4Radioac 119 theRadioactivityTables.push_back(rTable); 120 NSplit = 1; 121 AnalogueMC = true; 122 BRBias = true; 123 halflifethreshold = 1000.*nanosecond; 124 } 125 126 void G4RadioactiveDecay::ProcessDescription(st 127 { 128 outFile << "The G4RadioactiveDecay process p 129 << "nuclides (G4GenericIon) in biase 130 << "duplication, branching ratio bia 131 << "and detector time convolution. 132 << "activation physics.\n" 133 << "The required half-lives and deca 134 << "the RadioactiveDecay database wh 135 } 136 137 138 G4RadioactiveDecay::~G4RadioactiveDecay() 139 { 140 delete theRadioactivationMessenger; 141 } 142 143 G4bool 144 G4RadioactiveDecay::IsRateTableReady(const G4P 145 { 146 // Check whether the radioactive decay rates 147 // been calculated. 148 G4String aParticleName = aParticle.GetPartic 149 for (std::size_t i = 0; i < theParentChainTa 150 if (theParentChainTable[i].GetIonName() == 151 } 152 return false; 153 } 154 155 void 156 G4RadioactiveDecay::GetChainsFromParent(const 157 { 158 // Retrieve the decay rate table for the spe 159 G4String aParticleName = aParticle.GetPartic 160 161 for (std::size_t i = 0; i < theParentChainTa 162 if (theParentChainTable[i].GetIonName() == 163 theDecayRateVector = theParentChainTable 164 } 165 } 166 #ifdef G4VERBOSE 167 if (GetVerboseLevel() > 1) { 168 G4cout << "The DecayRate Table for " << aP 169 << G4endl; 170 } 171 #endif 172 } 173 174 // ConvolveSourceTimeProfile performs the conv 175 // function with a single exponential characte 176 // decay chain. The time profile is treated a 177 // convolution integral can be done bin-by-bin 178 // This implements Eq. 4.13 of DERA technical 179 180 G4double 181 G4RadioactiveDecay::ConvolveSourceTimeProfile( 182 { 183 G4double convolvedTime = 0.0; 184 G4int nbin; 185 if ( t > SBin[NSourceBin]) { 186 nbin = NSourceBin; 187 } else { 188 nbin = 0; 189 190 G4int loop = 0; 191 while (t > SBin[nbin]) { // Loop checking 192 loop++; 193 if (loop > 1000) { 194 G4Exception("G4RadioactiveDecay::Convo 195 "HAD_RDM_100", JustWarning 196 break; 197 } 198 ++nbin; 199 } 200 --nbin; 201 } 202 203 // Use expm1 wherever possible to avoid larg 204 // 1 - exp(x) for small x 205 G4double earg = 0.0; 206 if (nbin > 0) { 207 for (G4int i = 0; i < nbin; ++i) { 208 earg = (SBin[i+1] - SBin[i])/tau; 209 if (earg < 100.) { 210 convolvedTime += SProfile[i] * std::ex 211 std::expm1(earg); 212 } else { 213 convolvedTime += SProfile[i] * 214 (std::exp(-(t-SBin[i+1])/tau)-std::e 215 } 216 } 217 } 218 convolvedTime -= SProfile[nbin] * std::expm1 219 // tau divided out of final result to provid 220 221 if (convolvedTime < 0.) { 222 G4cout << " Convolved time =: " << convolv 223 G4cout << " t = " << t << " tau = " << tau 224 G4cout << SBin[nbin] << " " << SBin[0] << 225 convolvedTime = 0.; 226 } 227 #ifdef G4VERBOSE 228 if (GetVerboseLevel() > 2) 229 G4cout << " Convolved time: " << convolved 230 #endif 231 return convolvedTime; 232 } 233 234 235 ////////////////////////////////////////////// 236 // 237 // GetDecayTime 238 // Randomly select a decay time for the dec 239 // supplied decay time bias scheme. 240 // 241 ////////////////////////////////////////////// 242 243 G4double G4RadioactiveDecay::GetDecayTime() 244 { 245 G4double decaytime = 0.; 246 G4double rand = G4UniformRand(); 247 G4int i = 0; 248 249 G4int loop = 0; 250 while (DProfile[i] < rand) { /* Loop checki 251 // Entries in DProfile[i] are all between 252 // Comparison with rand chooses which time 253 ++i; 254 loop++; 255 if (loop > 100000) { 256 G4Exception("G4RadioactiveDecay::GetDeca 257 JustWarning, "While loop cou 258 break; 259 } 260 } 261 262 rand = G4UniformRand(); 263 decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i 264 #ifdef G4VERBOSE 265 if (GetVerboseLevel() > 2) 266 G4cout <<" Decay time: " <<decaytime/s <<" 267 #endif 268 return decaytime; 269 } 270 271 272 G4int G4RadioactiveDecay::GetDecayTimeBin(cons 273 { 274 G4int i = 0; 275 276 G4int loop = 0; 277 while (aDecayTime > DBin[i] ) { /* Loop ch 278 ++i; 279 loop++; 280 if (loop > 100000) { 281 G4Exception("G4RadioactiveDecay::GetDeca 282 JustWarning, "While loop cou 283 break; 284 } 285 } 286 287 return i; 288 } 289 290 ////////////////////////////////////////////// 291 // 292 // GetMeanLifeTime (required by the base clas 293 // 294 ////////////////////////////////////////////// 295 296 G4double G4RadioactiveDecay::GetMeanLifeTime(c 297 G4 298 { 299 // For variance reduction time is set to 0 s 300 // to decay immediately. 301 // In analogue mode it returns the particle' 302 G4double meanlife = 0.; 303 if (AnalogueMC) meanlife = G4VRadioactiveDec 304 return meanlife; 305 } 306 307 308 void 309 G4RadioactiveDecay::SetDecayRate(G4int theZ, G 310 G4int theG, std::vector<G4double>& the 311 std::vector<G4double>& theTaos) 312 // Why not make this a method of G4Radioactiv 313 { 314 //fill the decay rate vector 315 ratesToDaughter.SetZ(theZ); 316 ratesToDaughter.SetA(theA); 317 ratesToDaughter.SetE(theE); 318 ratesToDaughter.SetGeneration(theG); 319 ratesToDaughter.SetDecayRateC(theCoefficient 320 ratesToDaughter.SetTaos(theTaos); 321 } 322 323 324 void G4RadioactiveDecay:: 325 CalculateChainsFromParent(const G4ParticleDefi 326 { 327 // Use extended Bateman equation to calculat 328 // progeny of theParentNucleus. The coeffic 329 // calculated using the method of P. Truscot 330 // DERA Technical Note DERA/CIS/CIS2/7/36/4/ 331 // Coefficients are then added to the decay 332 333 // Create and initialise variables used in t 334 theDecayRateVector.clear(); 335 336 G4int nGeneration = 0; 337 338 std::vector<G4double> taos; 339 340 // Dimensionless A coefficients of Eqs. 4.24 341 std::vector<G4double> Acoeffs; 342 343 // According to Eq. 4.26 the first coefficie 344 Acoeffs.push_back(-1.); 345 346 const G4Ions* ion = static_cast<const G4Ions 347 G4int A = ion->GetAtomicMass(); 348 G4int Z = ion->GetAtomicNumber(); 349 G4double E = ion->GetExcitationEnergy(); 350 G4double tao = ion->GetPDGLifeTime(); 351 if (tao < 0.) tao = 1e-100; 352 taos.push_back(tao); 353 G4int nEntry = 0; 354 355 // Fill the decay rate container (G4Radioact 356 // isotope data 357 SetDecayRate(Z,A,E,nGeneration,Acoeffs,taos) 358 359 // store the decay rate in decay rate vector 360 theDecayRateVector.push_back(ratesToDaughter 361 ++nEntry; 362 363 // Now start treating the secondary generati 364 G4bool stable = false; 365 G4int j; 366 G4VDecayChannel* theChannel = 0; 367 G4NuclearDecay* theNuclearDecayChannel = 0; 368 369 G4ITDecay* theITChannel = 0; 370 G4BetaMinusDecay* theBetaMinusChannel = 0; 371 G4BetaPlusDecay* theBetaPlusChannel = 0; 372 G4AlphaDecay* theAlphaChannel = 0; 373 G4ProtonDecay* theProtonChannel = 0; 374 G4TritonDecay* theTritonChannel = 0; 375 G4NeutronDecay* theNeutronChannel = 0; 376 G4SFDecay* theFissionChannel = 0; 377 378 G4RadioactiveDecayMode theDecayMode; 379 G4double theBR = 0.0; 380 G4int AP = 0; 381 G4int ZP = 0; 382 G4int AD = 0; 383 G4int ZD = 0; 384 G4double EP = 0.; 385 std::vector<G4double> TP; 386 std::vector<G4double> RP; // A coefficient 387 G4ParticleDefinition *theDaughterNucleus; 388 G4double daughterExcitation; 389 G4double nearestEnergy = 0.0; 390 G4int nearestLevelIndex = 0; 391 G4ParticleDefinition *aParentNucleus; 392 G4IonTable* theIonTable; 393 G4DecayTable* parentDecayTable; 394 G4double theRate; 395 G4double TaoPlus; 396 G4int nS = 0; // Running index of fir 397 G4int nT = nEntry; // Total number of deca 398 const G4int nMode = G4RadioactiveDecayModeSi 399 G4double brs[nMode]; 400 // 401 theIonTable = G4ParticleTable::GetParticleTa 402 403 G4int loop = 0; 404 while (!stable) { /* Loop checking, 01.09. 405 loop++; 406 if (loop > 10000) { 407 G4Exception("G4RadioactiveDecay::Calcula 408 JustWarning, "While loop cou 409 break; 410 } 411 nGeneration++; 412 for (j = nS; j < nT; ++j) { 413 // First time through, get data for pare 414 ZP = theDecayRateVector[j].GetZ(); 415 AP = theDecayRateVector[j].GetA(); 416 EP = theDecayRateVector[j].GetE(); 417 RP = theDecayRateVector[j].GetDecayRateC 418 TP = theDecayRateVector[j].GetTaos(); 419 if (GetVerboseLevel() > 1) { 420 G4cout << "G4RadioactiveDecay::Calcula 421 << ZP << ", " << AP << ", " << 422 << ") are being calculated, ge 423 << G4endl; 424 } 425 426 aParentNucleus = theIonTable->GetIon(ZP, 427 parentDecayTable = GetDecayTable(aParent 428 if (nullptr == parentDecayTable) { conti 429 430 G4DecayTable* summedDecayTable = new G4D 431 // This instance of G4DecayTable is for 432 // channels. It will contain one decay 433 // (alpha, beta, etc.); its branching ra 434 // branching ratios for that type of dec 435 // halflife of a particular channel is l 436 // that channel will be inserted specifi 437 // ratio will not be included in the abo 438 // This instance is not used to perform 439 440 for (G4int k = 0; k < nMode; ++k) brs[k] 441 442 // Go through the decay table and sum al 443 for (G4int i = 0; i < parentDecayTable-> 444 theChannel = parentDecayTable->GetDeca 445 theNuclearDecayChannel = static_cast<G 446 theDecayMode = theNuclearDecayChannel- 447 daughterExcitation = theNuclearDecayCh 448 theDaughterNucleus = theNuclearDecayCh 449 AD = ((const G4Ions*)(theDaughterNucle 450 ZD = ((const G4Ions*)(theDaughterNucle 451 const G4LevelManager* levelManager = 452 G4NuclearLevelData::GetInstance()->G 453 454 // Check each nuclide to see if it is 455 // If so, add it to the decay chain by 456 // summedDecayTable. If not, just add 457 if (levelManager->NumberOfTransitions( 458 nearestEnergy = levelManager->Neares 459 if ((std::abs(daughterExcitation - n 460 (std::abs(daughterExcitation - nearest 461 // Level half-life is in ns and th 462 // by default, user can set it via 463 nearestLevelIndex = (G4int)levelMa 464 if (levelManager->LifeTime(nearest 465 // save the metastable decay cha 466 summedDecayTable->Insert(theChan 467 } else { 468 brs[theDecayMode] += theChannel- 469 } 470 } else { 471 brs[theDecayMode] += theChannel->G 472 } 473 } else { 474 brs[theDecayMode] += theChannel->Get 475 } 476 477 } // Combine decay channels (loop i) 478 479 brs[BetaPlus] = brs[BetaPlus]+brs[Kshell 480 brs[KshellEC] = brs[LshellEC] = brs[Mshe 481 for (G4int i = 0; i < nMode; ++i) { 482 if (brs[i] > 0.) { 483 switch (i) { 484 case IT: 485 // Decay mode is isomeric transiti 486 theITChannel = new G4ITDecay(aPare 487 488 summedDecayTable->Insert(theITChan 489 break; 490 491 case BetaMinus: 492 // Decay mode is beta- 493 theBetaMinusChannel = new G4BetaMi 494 495 496 summedDecayTable->Insert(theBetaMi 497 break; 498 499 case BetaPlus: 500 // Decay mode is beta+ + EC. 501 theBetaPlusChannel = new G4BetaPlu 502 503 504 summedDecayTable->Insert(theBetaPl 505 break; 506 507 case Alpha: 508 // Decay mode is alpha. 509 theAlphaChannel = new G4AlphaDecay 510 511 summedDecayTable->Insert(theAlphaC 512 break; 513 514 case Proton: 515 // Decay mode is proton. 516 theProtonChannel = new G4ProtonDec 517 518 summedDecayTable->Insert(theProton 519 break; 520 521 case Neutron: 522 // Decay mode is neutron. 523 theNeutronChannel = new G4NeutronD 524 525 summedDecayTable->Insert(theNeutro 526 break; 527 528 case SpFission: 529 // Decay mode is spontaneous fissi 530 theFissionChannel = new G4SFDecay( 531 532 summedDecayTable->Insert(theFissio 533 break; 534 535 case BDProton: 536 // Not yet implemented 537 break; 538 539 case BDNeutron: 540 // Not yet implemented 541 break; 542 543 case Beta2Minus: 544 // Not yet implemented 545 break; 546 547 case Beta2Plus: 548 // Not yet implemented 549 break; 550 551 case Proton2: 552 // Not yet implemented 553 break; 554 555 case Neutron2: 556 // Not yet implemented 557 break; 558 559 case Triton: 560 // Decay mode is Triton. 561 theTritonChannel = new G4TritonDec 562 563 summedDecayTable->Insert(theTriton 564 break; 565 566 default: 567 break; 568 } 569 } 570 } 571 572 // loop over all branches in summedDecay 573 // 574 for (G4int i = 0; i < summedDecayTable-> 575 theChannel = summedDecayTable->GetDeca 576 theNuclearDecayChannel = static_cast<G 577 theBR = theChannel->GetBR(); 578 theDaughterNucleus = theNuclearDecayCh 579 580 // First check if the decay of the ori 581 // if true create a new ground-state n 582 if (theNuclearDecayChannel->GetDecayMo 583 584 A = ((const G4Ions*)(theDaughterNucl 585 Z = ((const G4Ions*)(theDaughterNucl 586 theDaughterNucleus=theIonTable->GetI 587 } 588 if (IsApplicable(*theDaughterNucleus) 589 aParentNucleus != theDaughterNucle 590 // need to make sure daughter has de 591 parentDecayTable = GetDecayTable(the 592 if (nullptr != parentDecayTable && p 593 A = ((const G4Ions*)(theDaughterNu 594 Z = ((const G4Ions*)(theDaughterNu 595 E = ((const G4Ions*)(theDaughterNu 596 597 TaoPlus = theDaughterNucleus->GetP 598 if (TaoPlus <= 0.) TaoPlus = 1e-1 599 600 // first set the taos, one simply 601 taos.clear(); 602 taos = TP; // load lifetimes of 603 std::size_t k; 604 //check that TaoPlus differs from 605 //for (k = 0; k < TP.size(); ++k){ 606 //if (std::abs((TaoPlus-TP[k])/TP[ 607 //} 608 taos.push_back(TaoPlus); // add d 609 // now calculate the coefficiencie 610 // 611 // they are in two parts, first th 612 // Eq 4.24 of the TN 613 Acoeffs.clear(); 614 long double ta1,ta2; 615 ta2 = (long double)TaoPlus; 616 for (k = 0; k < RP.size(); ++k){ 617 ta1 = (long double)TP[k]; // 618 if (ta1 == ta2) { 619 theRate = 1.e100; 620 } else { 621 theRate = ta1/(ta1-ta2); 622 } 623 theRate = theRate * theBR * RP[k 624 Acoeffs.push_back(theRate); 625 } 626 627 // the second part: the n:n coeffi 628 // Eq 4.25 of the TN. Note Yn+1 i 629 // as treated at line 1013 630 theRate = 0.; 631 long double aRate, aRate1; 632 aRate1 = 0.L; 633 for (k = 0; k < RP.size(); ++k){ 634 ta1 = (long double)TP[k]; 635 if (ta1 == ta2 ) { 636 aRate = 1.e100; 637 } else { 638 aRate = ta2/(ta1-ta2); 639 } 640 aRate = aRate * (long double)(th 641 aRate1 += aRate; 642 } 643 theRate = -aRate1; 644 Acoeffs.push_back(theRate); 645 SetDecayRate (Z,A,E,nGeneration,Ac 646 theDecayRateVector.push_back(rates 647 nEntry++; 648 } // there are entries in the table 649 } // nuclide is OK to decay 650 } // end of loop (i) over decay table br 651 652 delete summedDecayTable; 653 654 } // Getting contents of decay rate vector 655 nS = nT; 656 nT = nEntry; 657 if (nS == nT) stable = true; 658 } // while nuclide is not stable 659 660 // end of while loop 661 // the calculation completed here 662 663 664 // fill the first part of the decay rate tab 665 // which is the name of the original particl 666 chainsFromParent.SetIonName(theParentNucleus 667 668 // now fill the decay table with the newly c 669 chainsFromParent.SetItsRates(theDecayRateVec 670 671 // finally add the decayratetable to the tab 672 theParentChainTable.push_back(chainsFromPare 673 } 674 675 ////////////////////////////////////////////// 676 // 677 // SetSourceTimeProfile 678 // read in the source time profile functio 679 // 680 ////////////////////////////////////////////// 681 682 void G4RadioactiveDecay::SetSourceTimeProfile( 683 { 684 std::ifstream infile ( filename, std::ios::i 685 if (!infile) { 686 G4ExceptionDescription ed; 687 ed << " Could not open file " << filename 688 G4Exception("G4RadioactiveDecay::SetSource 689 FatalException, ed); 690 } 691 692 G4double bin, flux; 693 NSourceBin = -1; 694 695 G4int loop = 0; 696 while (infile >> bin >> flux) { /* Loop che 697 loop++; 698 if (loop > 10000) { 699 G4Exception("G4RadioactiveDecay::SetSour 700 JustWarning, "While loop cou 701 break; 702 } 703 704 NSourceBin++; 705 if (NSourceBin > 99) { 706 G4Exception("G4RadioactiveDecay::SetSour 707 FatalException, "Input sourc 708 709 } else { 710 SBin[NSourceBin] = bin * s; // Conver 711 SProfile[NSourceBin] = flux; // Dimens 712 } 713 } 714 715 AnalogueMC = false; 716 infile.close(); 717 718 #ifdef G4VERBOSE 719 if (GetVerboseLevel() > 2) 720 G4cout <<" Source Timeprofile Nbin = " << 721 #endif 722 } 723 724 ////////////////////////////////////////////// 725 // 726 // SetDecayBiasProfile 727 // read in the decay bias scheme function ( 728 // 729 ////////////////////////////////////////////// 730 731 void G4RadioactiveDecay::SetDecayBias(const G4 732 { 733 std::ifstream infile(filename, std::ios::in) 734 if (!infile) G4Exception("G4RadioactiveDecay 735 FatalException, "Un 736 737 G4double bin, flux; 738 G4int dWindows = 0; 739 G4int i ; 740 741 theRadioactivityTables.clear(); 742 743 NDecayBin = -1; 744 745 G4int loop = 0; 746 while (infile >> bin >> flux ) { /* Loop ch 747 NDecayBin++; 748 loop++; 749 if (loop > 10000) { 750 G4Exception("G4RadioactiveDecay::SetDeca 751 JustWarning, "While loop cou 752 break; 753 } 754 755 if (NDecayBin > 99) { 756 G4Exception("G4RadioactiveDecay::SetDeca 757 FatalException, "Input bias 758 } else { 759 DBin[NDecayBin] = bin * s; // Conver 760 DProfile[NDecayBin] = flux; // Dimens 761 if (flux > 0.) { 762 decayWindows[NDecayBin] = dWindows; 763 dWindows++; 764 G4RadioactivityTable *rTable = new G4R 765 theRadioactivityTables.push_back(rTabl 766 } 767 } 768 } 769 for ( i = 1; i<= NDecayBin; ++i) DProfile[i] 770 for ( i = 0; i<= NDecayBin; ++i) DProfile[i] 771 // Normalize so e 772 // converted to accumulated probabilities 773 774 AnalogueMC = false; 775 infile.close(); 776 777 #ifdef G4VERBOSE 778 if (GetVerboseLevel() > 2) 779 G4cout <<" Decay Bias Profile Nbin = " << 780 #endif 781 } 782 783 ////////////////////////////////////////////// 784 // 785 // DecayIt 786 // 787 ////////////////////////////////////////////// 788 789 G4VParticleChange* 790 G4RadioactiveDecay::DecayIt(const G4Track& the 791 { 792 // Initialize G4ParticleChange object, get p 793 fParticleChangeForRadDecay.Initialize(theTra 794 fParticleChangeForRadDecay.ProposeWeight(the 795 const G4DynamicParticle* theParticle = theTr 796 const G4ParticleDefinition* theParticleDef = 797 798 // First check whether RDM applies to the cu 799 if (!isAllVolumesMode) { 800 if (!std::binary_search(ValidVolumes.begin 801 theTrack.GetVolume()->Get 802 #ifdef G4VERBOSE 803 if (GetVerboseLevel()>0) { 804 G4cout <<"G4RadioactiveDecay::DecayIt 805 << theTrack.GetVolume()->GetLog 806 << " is not selected for the RD 807 G4cout << " There are " << ValidVolume 808 G4cout << " The Valid volumes are " << 809 for (std::size_t i = 0; i< ValidVolume 810 G4cout << Va 811 } 812 #endif 813 fParticleChangeForRadDecay.SetNumberOfSe 814 815 // Kill the parent particle. 816 fParticleChangeForRadDecay.ProposeTrackS 817 fParticleChangeForRadDecay.ProposeLocalE 818 ClearNumberOfInteractionLengthLeft(); 819 return &fParticleChangeForRadDecay; 820 } 821 } 822 823 // Now check if particle is valid for RDM 824 if (!(IsApplicable(*theParticleDef) ) ) { 825 // Particle is not an ion or is outside th 826 #ifdef G4VERBOSE 827 if (GetVerboseLevel() > 1) { 828 G4cout << "G4RadioactiveDecay::DecayIt : 829 << theParticleDef->GetParticleNam 830 << " is not an ion or is outside 831 << " Set particle change accordin 832 << G4endl; 833 } 834 #endif 835 fParticleChangeForRadDecay.SetNumberOfSeco 836 837 // Kill the parent particle 838 fParticleChangeForRadDecay.ProposeTrackSta 839 fParticleChangeForRadDecay.ProposeLocalEne 840 ClearNumberOfInteractionLengthLeft(); 841 return &fParticleChangeForRadDecay; 842 } 843 844 G4DecayTable* theDecayTable = GetDecayTable( 845 846 if (theDecayTable == nullptr || theDecayTabl 847 // No data in the decay table. Set partic 848 // to indicate this. 849 #ifdef G4VERBOSE 850 if (GetVerboseLevel() > 1) { 851 G4cout << "G4RadioactiveDecay::DecayIt : 852 << "decay table not defined for " 853 << theParticleDef->GetParticleNam 854 << ". Set particle change accordi 855 << G4endl; 856 } 857 #endif 858 fParticleChangeForRadDecay.SetNumberOfSeco 859 860 // Kill the parent particle. 861 fParticleChangeForRadDecay.ProposeTrackSta 862 fParticleChangeForRadDecay.ProposeLocalEne 863 ClearNumberOfInteractionLengthLeft(); 864 return &fParticleChangeForRadDecay; 865 866 } else { 867 // Data found. Try to decay nucleus 868 if (AnalogueMC) { 869 G4VRadioactiveDecay::DecayAnalog(theTrac 870 871 } else { 872 // Proceed with decay using variance red 873 G4double energyDeposit = 0.0; 874 G4double finalGlobalTime = theTrack.GetG 875 G4double finalLocalTime = theTrack.GetLo 876 G4int index; 877 G4ThreeVector currentPosition; 878 currentPosition = theTrack.GetPosition() 879 880 G4IonTable* theIonTable; 881 G4ParticleDefinition* parentNucleus; 882 883 // Get decay chains for the given nuclid 884 if (!IsRateTableReady(*theParticleDef)) 885 CalculateChainsFromParent(*theParticleDef); 886 GetChainsFromParent(*theParticleDef); 887 888 // Declare some of the variables require 889 G4int PZ; 890 G4int PA; 891 G4double PE; 892 G4String keyName; 893 std::vector<G4double> PT; 894 std::vector<G4double> PR; 895 G4double taotime; 896 long double decayRate; 897 898 std::size_t i; 899 G4int numberOfSecondaries; 900 G4int totalNumberOfSecondaries = 0; 901 G4double currentTime = 0.; 902 G4int ndecaych; 903 G4DynamicParticle* asecondaryparticle; 904 std::vector<G4DynamicParticle*> secondar 905 std::vector<G4double> pw; 906 std::vector<G4double> ptime; 907 pw.clear(); 908 ptime.clear(); 909 910 // Now apply the nucleus splitting 911 for (G4int n = 0; n < NSplit; ++n) { 912 // Get the decay time following the de 913 // supplied by user 914 G4double theDecayTime = GetDecayTime() 915 G4int nbin = GetDecayTimeBin(theDecayT 916 917 // calculate the first part of the wei 918 G4double weight1 = 1.; 919 if (nbin == 1) { 920 weight1 = 1./DProfile[nbin-1] 921 *(DBin[nbin]-DBin[nbin-1]) 922 } else if (nbin > 1) { 923 // Go from nbin to nbin-2 because fl 924 weight1 = 1./(DProfile[nbin]-DProfil 925 *(DBin[nbin]-DBin[nbin-1]) 926 // weight1 = (probability of choosin 927 } 928 // it should be calculated in seconds 929 weight1 /= s ; 930 931 // loop over all the possible secondar 932 // the first one is itself. 933 for (i = 0; i < theDecayRateVector.siz 934 PZ = theDecayRateVector[i].GetZ(); 935 PA = theDecayRateVector[i].GetA(); 936 PE = theDecayRateVector[i].GetE(); 937 PT = theDecayRateVector[i].GetTaos() 938 PR = theDecayRateVector[i].GetDecayR 939 940 // The array of arrays theDecayRateV 941 // chains of a given parent nucleus 942 // nuclide (Z,A,E). 943 // 944 // theDecayRateVector[0] contains th 945 // nucleus 946 // PZ = ZP 947 // PA = AP 948 // PE = EP 949 // PT[] = {TP} 950 // PR[] = {RP} 951 // 952 // theDecayRateVector[1] contains th 953 // generation daughter (Z1,A1,E1). 954 // PZ = Z1 955 // PA = A1 956 // PE = E1 957 // PT[] = {TP, T1} 958 // PR[] = {RP, R1} 959 // 960 // theDecayRateVector[2] contains th 961 // generation daughter (Z1,A1,E1) an 962 // generation daughter to the second 963 // PZ = Z2 964 // PA = A2 965 // PE = E2 966 // PT[] = {TP, T1, T2} 967 // PR[] = {RP, R1, R2} 968 // 969 // theDecayRateVector[3] may contain 970 // PZ = Z2a 971 // PA = A2a 972 // PE = E2a 973 // PT[] = {TP, T1, T2a} 974 // PR[] = {RP, R1, R2a} 975 // 976 // and so on. 977 978 // Calculate the decay rate of the i 979 // radioactivity of isotope (PZ,PA,P 980 // it will be used to calculate the 981 // decay products of this isotope 982 983 // For each nuclide, calculate all t 984 // the parent nuclide 985 decayRate = 0.L; 986 for (G4int j = 0; j < G4int(PT.size( 987 taotime = ConvolveSourceTimeProfil 988 decayRate -= PR[j] * (long double) 989 // Eq.4.23 of of the TN 990 // note the negative here is requi 991 // equation is defined to be negat 992 // i.e. decay away, but we need po 993 994 // G4cout << j << "\t"<< PT[j]/s < 995 } 996 997 // At this point any negative decay 998 // (order 10**-30) that negative val 999 // errors. Set them to zero. 1000 if (decayRate < 0.0) decayRate = 0. 1001 1002 // G4cout <<theDecayTime/s <<"\t"< 1003 // G4cout << theTrack.GetWeight() 1004 1005 // Add isotope to the radioactivity 1006 // One table for each observation t 1007 // SetDecayBias(G4String filename) 1008 1009 theRadioactivityTables[decayWindows 1010 ->AddIsotope(PZ,PA,PE,weight1 1011 1012 // Now calculate the statistical we 1013 // One needs to fold the source bia 1014 // also need to include the track w 1015 G4double weight = weight1*decayRate 1016 1017 // decay the isotope 1018 theIonTable = (G4IonTable *)(G4Part 1019 parentNucleus = theIonTable->GetIon 1020 1021 // Create a temprary products buffe 1022 // Its contents to be transfered to 1023 G4DecayProducts* tempprods = nullpt 1024 1025 // Decide whether to apply branchin 1026 if (BRBias) { 1027 G4DecayTable* decayTable = GetDec 1028 G4VDecayChannel* theDecayChannel = null 1029 if (nullptr != decayTable) { 1030 ndecaych = G4int(decayTable->entries( 1031 theDecayChannel = decayTable->G 1032 } 1033 1034 if (theDecayChannel == nullptr) { 1035 // Decay channel not found. 1036 1037 if (GetVerboseLevel() > 0) { 1038 G4cout << " G4RadioactiveDeca 1039 G4cout << " for this nucleus; 1040 G4cout << G4endl; 1041 if (nullptr != decayTable) { 1042 } 1043 // DHW 6 Dec 2010 - do decay as if no 1044 tempprods = DoDecay(*parentNucl 1045 } else { 1046 // A decay channel has been ide 1047 G4double tempmass = parentNucle 1048 tempprods = theDecayChannel->De 1049 weight *= (theDecayChannel->Get 1050 } 1051 } else { 1052 tempprods = DoDecay(*parentNucleu 1053 } 1054 1055 // save the secondaries for buffers 1056 numberOfSecondaries = tempprods->en 1057 currentTime = finalGlobalTime + the 1058 for (index = 0; index < numberOfSec 1059 asecondaryparticle = tempprods->P 1060 if (asecondaryparticle->GetDefini 1061 pw.push_back(weight); 1062 ptime.push_back(currentTime); 1063 secondaryparticles.push_back(as 1064 } 1065 // Generate gammas and Xrays from 1066 else if (((const G4Ions*)(asecond 1067 ->GetExcitationEnergy() 1068 G4ParticleDefinition* apartDef 1069 AddDeexcitationSpectrumForBiasM 1070 1071 } 1072 } 1073 1074 delete tempprods; 1075 } // end of i loop 1076 } // end of n loop 1077 1078 // now deal with the secondaries in the 1079 // and submmit them back to the trackin 1080 totalNumberOfSecondaries = (G4int)pw.si 1081 fParticleChangeForRadDecay.SetNumberOfS 1082 for (index=0; index < totalNumberOfSeco 1083 G4Track* secondary = new G4Track(seco 1084 ptim 1085 secondary->SetGoodForTrackingFlag(); 1086 secondary->SetTouchableHandle(theTrac 1087 secondary->SetWeight(pw[index]); 1088 fParticleChangeForRadDecay.AddSeconda 1089 } 1090 1091 // Kill the parent particle 1092 fParticleChangeForRadDecay.ProposeTrack 1093 fParticleChangeForRadDecay.ProposeLocal 1094 fParticleChangeForRadDecay.ProposeLocal 1095 // Reset NumberOfInteractionLengthLeft. 1096 ClearNumberOfInteractionLengthLeft(); 1097 } // end VR decay 1098 1099 return &fParticleChangeForRadDecay; 1100 } // end of data found branch 1101 } 1102 1103 1104 // Add gamma, X-ray, conversion and auger ele 1105 void 1106 G4RadioactiveDecay::AddDeexcitationSpectrumFo 1107 G4dou 1108 std:: 1109 std:: 1110 std:: 1111 { 1112 G4double elevel=((const G4Ions*)(apartDef)) 1113 G4double life_time=apartDef->GetPDGLifeTime 1114 G4ITDecay* anITChannel = 0; 1115 1116 while (life_time < halflifethreshold && ele 1117 decayIT->SetupDecay(apartDef); 1118 G4DecayProducts* pevap_products = decayIT 1119 G4int nb_pevapSecondaries = pevap_product 1120 1121 G4DynamicParticle* a_pevap_secondary = 0; 1122 G4ParticleDefinition* secDef = 0; 1123 for (G4int ind = 0; ind < nb_pevapSeconda 1124 a_pevap_secondary= pevap_products->PopP 1125 secDef = a_pevap_secondary->GetDefiniti 1126 1127 if (secDef->GetBaryonNumber() > 4) { 1128 elevel = ((const G4Ions*)(secDef))->G 1129 life_time = secDef->GetPDGLifeTime(); 1130 apartDef = secDef; 1131 if (secDef->GetPDGStable() ) { 1132 weights_v.push_back(weight); 1133 times_v.push_back(currentTime); 1134 secondaries_v.push_back(a_pevap_sec 1135 } 1136 } else { 1137 weights_v.push_back(weight); 1138 times_v.push_back(currentTime); 1139 secondaries_v.push_back(a_pevap_secon 1140 } 1141 } 1142 1143 delete anITChannel; 1144 delete pevap_products; 1145 } 1146 } 1147 1148