00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00064
00065 #include "G4PhysicalConstants.hh"
00066 #include "G4SystemOfUnits.hh"
00067 #include "G4NuclearLevelManager.hh"
00068 #include "G4NuclearLevelStore.hh"
00069 #include "G4NuclearDecayChannel.hh"
00070 #include "G4DynamicParticle.hh"
00071 #include "G4DecayProducts.hh"
00072 #include "G4DecayTable.hh"
00073 #include "G4PhysicsLogVector.hh"
00074 #include "G4ParticleChangeForRadDecay.hh"
00075 #include "G4IonTable.hh"
00076 #include "G4HadTmpUtil.hh"
00077
00078 #include "G4BetaFermiFunction.hh"
00079 #include "G4PhotonEvaporation.hh"
00080
00081 #include "G4VAtomDeexcitation.hh"
00082 #include "G4AtomicShells.hh"
00083 #include "G4LossTableManager.hh"
00084
00085
00086
00087
00088
00089 const G4double G4NuclearDecayChannel:: pTolerance = 0.001;
00090 const G4double G4NuclearDecayChannel:: levelTolerance = 2.0*keV;
00091
00092
00093
00094
00095
00096 G4NuclearDecayChannel::
00097 G4NuclearDecayChannel(const G4RadioactiveDecayMode& theMode,
00098 G4int Verbose,
00099 const G4ParticleDefinition* theParentNucleus,
00100 G4double theBR,
00101 G4double theQtransition,
00102 G4int A,
00103 G4int Z,
00104 G4double theDaughterExcitation)
00105 :G4GeneralPhaseSpaceDecay(Verbose), decayMode(theMode), dynamicDaughter(0),
00106 RandomEnergy(0)
00107 {
00108 #ifdef G4VERBOSE
00109 if (GetVerboseLevel()>1)
00110 {G4cout <<"G4NuclearDecayChannel constructor for " <<G4int(theMode) <<G4endl;}
00111 #endif
00112 SetParent(theParentNucleus);
00113 FillParent();
00114 parent_mass = theParentNucleus->GetPDGMass();
00115 SetBR (theBR);
00116 SetNumberOfDaughters (1);
00117 FillDaughterNucleus (0, A, Z, theDaughterExcitation);
00118 Qtransition = theQtransition;
00119 halflifethreshold = -1.*second;
00120 applyICM = true;
00121 applyARM = true;
00122 }
00123
00124
00125
00126
00127 G4NuclearDecayChannel::
00128 G4NuclearDecayChannel(const G4RadioactiveDecayMode& theMode,
00129 G4int Verbose,
00130 const G4ParticleDefinition *theParentNucleus,
00131 G4double theBR,
00132 G4double theQtransition,
00133 G4int A,
00134 G4int Z,
00135 G4double theDaughterExcitation,
00136 const G4String theDaughterName1)
00137 :G4GeneralPhaseSpaceDecay(Verbose), decayMode(theMode), dynamicDaughter(0),
00138 RandomEnergy(0)
00139 {
00140 #ifdef G4VERBOSE
00141 if (GetVerboseLevel()>1)
00142 {G4cout <<"G4NuclearDecayChannel constructor for " <<G4int(theMode) <<G4endl;}
00143 #endif
00144 SetParent (theParentNucleus);
00145 FillParent();
00146 parent_mass = theParentNucleus->GetPDGMass();
00147 SetBR (theBR);
00148 SetNumberOfDaughters (2);
00149 SetDaughter(0, theDaughterName1);
00150 FillDaughterNucleus (1, A, Z, theDaughterExcitation);
00151 Qtransition = theQtransition;
00152 halflifethreshold = -1.*second;
00153 applyICM = true;
00154 applyARM = true;
00155 }
00156
00157
00158
00159
00160 G4NuclearDecayChannel::
00161 G4NuclearDecayChannel(const G4RadioactiveDecayMode &theMode,
00162 G4int Verbose,
00163 const G4ParticleDefinition *theParentNucleus,
00164 G4double theBR,
00165 G4double ,
00166 G4bool ,
00167 CLHEP::RandGeneral* randBeta,
00168 G4double theQtransition,
00169 G4int A,
00170 G4int Z,
00171 G4double theDaughterExcitation,
00172 const G4String theDaughterName1,
00173 const G4String theDaughterName2)
00174 :G4GeneralPhaseSpaceDecay(Verbose), decayMode(theMode), dynamicDaughter(0)
00175 {
00176 #ifdef G4VERBOSE
00177 if (GetVerboseLevel()>1)
00178 {G4cout <<"G4NuclearDecayChannel constructor for " <<G4int(theMode) <<G4endl;}
00179 #endif
00180 SetParent (theParentNucleus);
00181 FillParent();
00182 parent_mass = theParentNucleus->GetPDGMass();
00183 SetBR (theBR);
00184 SetNumberOfDaughters (3);
00185 SetDaughter(0, theDaughterName1);
00186 SetDaughter(2, theDaughterName2);
00187 FillDaughterNucleus(1, A, Z, theDaughterExcitation);
00188 RandomEnergy = randBeta;
00189 Qtransition = theQtransition;
00190 halflifethreshold = -1*second;
00191 applyICM = true;
00192 applyARM = true;
00193 }
00194
00195
00196 void G4NuclearDecayChannel::FillDaughterNucleus(G4int index, G4int A, G4int Z,
00197 G4double theDaughterExcitation)
00198 {
00199
00200
00201 if (A<1 || Z<0 || theDaughterExcitation <0.0)
00202 {
00203
00204
00205
00206
00207 G4ExceptionDescription ed;
00208 ed << "Inappropriate values of daughter A, Z or excitation: "
00209 << A << " , " << Z << " , " << theDaughterExcitation*MeV << " MeV "
00210 << G4endl;
00211 G4Exception("G4NuclearDecayChannel::FillDaughterNucleus()", "HAD_RDM_006",
00212 FatalException, ed);
00213
00214
00215 }
00216
00217
00218
00219 daughterA = A;
00220 daughterZ = Z;
00221 if (Z == 1 && A == 1) {
00222 daughterNucleus = G4Proton::Definition();
00223 } else if (Z == 0 && A == 1) {
00224 daughterNucleus = G4Neutron::Definition();
00225 } else {
00226 G4IonTable *theIonTable =
00227 (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
00228 daughterNucleus = theIonTable->GetIon(daughterZ, daughterA, theDaughterExcitation*MeV);
00229 }
00230 daughterExcitation = theDaughterExcitation;
00231 SetDaughter(index, daughterNucleus);
00232 }
00233
00234
00235 G4DecayProducts* G4NuclearDecayChannel::DecayIt(G4double theParentMass)
00236 {
00237
00238
00239
00240 if (parent == 0) FillParent();
00241 if (daughters == 0) FillDaughters();
00242
00243
00244
00245
00246 theParentMass = 0.0;
00247 for( G4int index=0; index < numberOfDaughters; index++)
00248 {theParentMass += daughters[index]->GetPDGMass();}
00249 theParentMass += Qtransition ;
00250
00251 if (decayMode == 2) theParentMass -= 2*0.511 * MeV;
00252
00253 #ifdef G4VERBOSE
00254 if (GetVerboseLevel()>1) {
00255 G4cout << "G4NuclearDecayChannel::DecayIt "<< G4endl;
00256 G4cout << "the decay mass = " << theParentMass << G4endl;
00257 }
00258 #endif
00259
00260 SetParentMass (theParentMass);
00261
00262
00263 G4DecayProducts* products = 0;
00264
00265
00266
00267 switch (numberOfDaughters) {
00268 case 0:
00269 G4cerr << "G4NuclearDecayChannel::DecayIt ";
00270 G4cerr << " daughters not defined " <<G4endl;
00271 break;
00272 case 1:
00273 products = OneBodyDecayIt();
00274 break;
00275 case 2:
00276 products = TwoBodyDecayIt();
00277 break;
00278 case 3:
00279 products = BetaDecayIt();
00280 break;
00281 default:
00282
00283
00284
00285 {
00286 G4ExceptionDescription ed;
00287 ed << " More than 3 daughters in decay: N = " << numberOfDaughters
00288 << G4endl;
00289 G4Exception("G4NuclearDecayChannel::DecayIt()", "HAD_RDM_007",
00290 FatalException, ed);
00291 }
00292 }
00293
00294 if (products == 0) {
00295 G4ExceptionDescription ed;
00296 ed << " Parent nucleus " << *parent_name << " was not decayed " << G4endl;
00297 G4Exception("G4NuclearDecayChannel::DecayIt()", "HAD_RDM_008",
00298 JustWarning, ed);
00299 DumpInfo();
00300 } else {
00301
00302
00303
00304
00305
00306 G4int shellIndex = -1;
00307
00308 if (daughterExcitation > 0.0) {
00309
00310
00311
00312 dynamicDaughter = products->PopProducts();
00313 G4LorentzVector daughterMomentum = dynamicDaughter->Get4Momentum();
00314 G4ThreeVector const daughterMomentum1(static_cast<const G4LorentzVector> (daughterMomentum));
00315
00316
00317
00318 G4Fragment nucleus(daughterA, daughterZ, daughterMomentum);
00319 G4PhotonEvaporation* deexcitation = new G4PhotonEvaporation;
00320 deexcitation->SetVerboseLevel(GetVerboseLevel());
00321
00322
00323 deexcitation->SetICM(applyICM);
00324
00325
00326
00327
00328 deexcitation->SetMaxHalfLife(halflifethreshold);
00329
00330
00331 if (decayMode == 0) {
00332 deexcitation->RDMForced(true);
00333 } else {
00334 deexcitation->RDMForced(false);
00335 }
00336
00338
00339
00340
00341
00342
00344
00345
00346 G4FragmentVector* gammas = 0;
00347 if (applyICM) {
00348 gammas = deexcitation->BreakUp(nucleus);
00349 } else {
00350 gammas = deexcitation->BreakItUp(nucleus);
00351 }
00352
00353
00354
00355 G4int nGammas=gammas->size()-1;
00356
00357
00358
00359
00360 for (G4int ig = 0; ig < nGammas; ig++) {
00361 G4DynamicParticle* theGammaRay =
00362 new G4DynamicParticle(gammas->operator[](ig)->GetParticleDefinition(),
00363 gammas->operator[](ig)->GetMomentum());
00364 theGammaRay -> SetProperTime(gammas->operator[](ig)->GetCreationTime());
00365 products->PushProducts (theGammaRay);
00366 }
00367
00368
00369 G4double finalDaughterExcitation = gammas->operator[](nGammas)->GetExcitationEnergy();
00370
00371 if (finalDaughterExcitation <= 1.0*keV) finalDaughterExcitation = 0;
00372
00373 if (dynamicDaughter) delete dynamicDaughter;
00374
00375 G4IonTable* theIonTable =
00376 (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
00377 dynamicDaughter = new G4DynamicParticle(
00378 theIonTable->GetIon(daughterZ,daughterA,finalDaughterExcitation),
00379 daughterMomentum1);
00380 products->PushProducts (dynamicDaughter);
00381
00382
00383 shellIndex = deexcitation->GetVacantShellNumber();
00384
00385
00386 while (!gammas->empty()) {
00387 delete *(gammas->end()-1);
00388 gammas->pop_back();
00389 }
00390 delete gammas;
00391 delete deexcitation;
00392 }
00393
00394
00395 G4int eShell = -1;
00396 if (decayMode == 3 || decayMode == 4 || decayMode == 5) {
00397 switch (decayMode)
00398 {
00399 case KshellEC:
00400
00401 {
00402 eShell = 0;
00403 }
00404 break;
00405 case LshellEC:
00406
00407 {
00408 eShell = G4int(G4UniformRand()*3)+1;
00409 }
00410 break;
00411 case MshellEC:
00412
00413 {
00414
00415
00416 eShell = G4int(G4UniformRand()*3)+4;
00417 }
00418 break;
00419 case ERROR:
00420 default:
00421 G4Exception("G4NuclearDecayChannel::DecayIt()", "HAD_RDM_009",
00422 FatalException, "Incorrect decay mode selection");
00423 }
00424 } else {
00425
00426
00427 eShell = shellIndex;
00428 }
00429
00430
00431 if (applyARM && eShell != -1) {
00432 G4int aZ = daughterZ;
00433
00434
00435
00436
00437 G4VAtomDeexcitation* atomDeex =
00438 G4LossTableManager::Instance()->AtomDeexcitation();
00439 if (atomDeex) {
00440 if(atomDeex->IsFluoActive() && aZ > 5 && aZ < 100) {
00441 if (eShell >= G4AtomicShells::GetNumberOfShells(aZ)){
00442 eShell = G4AtomicShells::GetNumberOfShells(aZ)-1;
00443 }
00444 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(eShell);
00445 const G4AtomicShell* shell = atomDeex->GetAtomicShell(aZ, as);
00446 std::vector<G4DynamicParticle*> armProducts;
00447 const G4double deexLimit = 0.1*keV;
00448 atomDeex->GenerateParticles(&armProducts, shell, aZ, deexLimit, deexLimit);
00449 size_t narm = armProducts.size();
00450 if (narm > 0) {
00451
00452
00453 dynamicDaughter = products->PopProducts();
00454 G4ThreeVector bst = dynamicDaughter->Get4Momentum().boostVector();
00455 for (size_t i = 0; i<narm; ++i) {
00456 G4DynamicParticle* dp = armProducts[i];
00457 G4LorentzVector lv = dp->Get4Momentum().boost(bst);
00458 dp->Set4Momentum(lv);
00459 products->PushProducts(dp);
00460 }
00461 products->PushProducts(dynamicDaughter);
00462 }
00463 }
00464 }
00465 }
00466 }
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520 return products;
00521 }
00522
00523
00524 G4DecayProducts* G4NuclearDecayChannel::BetaDecayIt()
00525 {
00526 if (GetVerboseLevel()>1) G4cout << "G4Decay::BetaDecayIt()"<<G4endl;
00527
00528
00529 G4double daughtermass[3];
00530 G4double sumofdaughtermass = 0.0;
00531 G4double pmass = GetParentMass();
00532 for (G4int index=0; index<3; index++)
00533 {
00534 daughtermass[index] = daughters[index]->GetPDGMass();
00535 sumofdaughtermass += daughtermass[index];
00536 }
00537
00538
00539 G4ParticleMomentum dummy;
00540 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00541
00542
00543 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00544 delete parentparticle;
00545
00546 G4double Q = pmass - sumofdaughtermass;
00547
00548
00549 G4double daughtermomentum[3];
00550 G4double daughterenergy[3];
00551
00552 daughterenergy[0] = RandomEnergy->shoot() * Q;
00553 daughtermomentum[0] = std::sqrt(daughterenergy[0]*daughterenergy[0] +
00554 2.0*daughterenergy[0] * daughtermass[0]);
00555
00556
00557 G4double rd = 2*G4UniformRand()-1;
00558
00559 G4double Mme=pmass-daughtermass[0];
00560 G4double K=0.5-daughtermass[1]*daughtermass[1]/(2*Mme*Mme-4*pmass*daughterenergy[0]);
00561
00562 daughterenergy[2]=K*(Mme-daughterenergy[0]+rd*daughtermomentum[0]);
00563 daughtermomentum[2] = daughterenergy[2] ;
00564
00565
00566 daughterenergy[1] = Q-daughterenergy[0]-daughterenergy[2];
00567 G4double recoilmomentumsquared = daughterenergy[1]*daughterenergy[1] +
00568 2.0*daughterenergy[1] * daughtermass[1];
00569 if (recoilmomentumsquared < 0.0) recoilmomentumsquared = 0.0;
00570 daughtermomentum[1] = std::sqrt(recoilmomentumsquared);
00571
00572
00573 if (GetVerboseLevel()>1) {
00574 G4cout <<" daughter 0:" <<daughtermomentum[0]/GeV <<"[GeV/c]" <<G4endl;
00575 G4cout <<" daughter 1:" <<daughtermomentum[1]/GeV <<"[GeV/c]" <<G4endl;
00576 G4cout <<" daughter 2:" <<daughtermomentum[2]/GeV <<"[GeV/c]" <<G4endl;
00577 }
00578
00579 G4double costheta, sintheta, phi, sinphi, cosphi;
00580 G4double costhetan, sinthetan, phin, sinphin, cosphin;
00581 costheta = 2.*G4UniformRand()-1.0;
00582 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00583 phi = twopi*G4UniformRand()*rad;
00584 sinphi = std::sin(phi);
00585 cosphi = std::cos(phi);
00586 G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta);
00587 G4DynamicParticle * daughterparticle
00588 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
00589 products->PushProducts(daughterparticle);
00590
00591 costhetan = (daughtermomentum[1]*daughtermomentum[1]-
00592 daughtermomentum[2]*daughtermomentum[2]-
00593 daughtermomentum[0]*daughtermomentum[0])/
00594 (2.0*daughtermomentum[2]*daughtermomentum[0]);
00595
00596
00597 if (costhetan > 1.) costhetan = 1.;
00598 if (costhetan < -1.) costhetan = -1.;
00599 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00600 phin = twopi*G4UniformRand()*rad;
00601 sinphin = std::sin(phin);
00602 cosphin = std::cos(phin);
00603 G4ParticleMomentum direction2;
00604 direction2.setX(sinthetan*cosphin*costheta*cosphi -
00605 sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
00606 direction2.setY(sinthetan*cosphin*costheta*sinphi +
00607 sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
00608 direction2.setZ(-sinthetan*cosphin*sintheta + costhetan*costheta);
00609 daughterparticle = new G4DynamicParticle(daughters[2],
00610 direction2*(daughtermomentum[2]/direction2.mag()));
00611 products->PushProducts(daughterparticle);
00612
00613 daughterparticle =
00614 new G4DynamicParticle(daughters[1],
00615 (direction0*daughtermomentum[0] +
00616 direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0));
00617 products->PushProducts(daughterparticle);
00618
00619 if (GetVerboseLevel()>1) {
00620 G4cout << "G4NuclearDecayChannel::BetaDecayIt ";
00621 G4cout << " create decay products in rest frame " <<G4endl;
00622 products->DumpInfo();
00623 }
00624 return products;
00625 }