Changeset 962 for trunk/source/processes/hadronic/models/radioactive_decay/src/G4NuclearDecayChannel.cc
- Timestamp:
- Apr 6, 2009, 12:30:29 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/radioactive_decay/src/G4NuclearDecayChannel.cc
r819 r962 66 66 #include "G4BetaFermiFunction.hh" 67 67 #include "G4PhotonEvaporation.hh" 68 #include "G4AtomicTransitionManager.hh" 69 #include "G4AtomicShell.hh" 68 70 #include "G4AtomicDeexcitation.hh" 69 70 71 71 72 const G4double G4NuclearDecayChannel:: pTolerance = 0.001; … … 241 242 if (decayMode == 2) theParentMass -= 2*0.511 * MeV; 242 243 // 244 #ifdef G4VERBOSE 243 245 if (GetVerboseLevel()>1) { 244 G4cout << "G4NuclearDecayChannel::DecayIt " ;246 G4cout << "G4NuclearDecayChannel::DecayIt "<< G4endl; 245 247 G4cout << "the decay mass = " << theParentMass << G4endl; 246 248 } 249 #endif 247 250 248 251 SetParentMass (theParentMass); … … 261 264 { 262 265 case 0: 263 if (GetVerboseLevel()>0) 264 { 265 G4cout << "G4NuclearDecayChannel::DecayIt "; 266 G4cout << " daughters not defined " <<G4endl; 267 } 266 G4cerr << "G4NuclearDecayChannel::DecayIt "; 267 G4cerr << " daughters not defined " <<G4endl; 268 268 break; 269 269 case 1: … … 281 281 G4Exception(__FILE__, G4inttostring(__LINE__), FatalException, "G4NuclearDecayChannel::DecayIt"); 282 282 } 283 if ( (products == 0) && (GetVerboseLevel()>0)) {283 if (products == 0) { 284 284 G4cerr << "G4NuclearDecayChannel::DecayIt "; 285 285 G4cerr << *parent_name << " can not decay " << G4endl; 286 286 DumpInfo(); 287 287 } 288 289 // 290 // now we have to take care of the EC product which have go through the ARM 288 // 289 // If the decay is to an excited state of the daughter nuclide, we need 290 // to apply the photo-evaporation process. 291 // 292 // needed to hold the shell idex after ICM 293 G4int shellIndex = -1; 294 // 295 if (daughterExcitation > 0.0) 296 { 297 // 298 // Pop the daughter nucleus off the product vector - we need to retain 299 // the momentum of this particle. 300 // 301 dynamicDaughter = products->PopProducts(); 302 G4LorentzVector daughterMomentum = dynamicDaughter->Get4Momentum(); 303 G4ThreeVector const daughterMomentum1(static_cast<const G4LorentzVector> (daughterMomentum)); 304 // 305 // 306 // Now define a G4Fragment with the correct A, Z and excitation, and declare and 307 // initialise a G4PhotonEvaporation object. 308 // 309 G4Fragment nucleus(daughterA, daughterZ, daughterMomentum); 310 G4PhotonEvaporation* deexcitation = new G4PhotonEvaporation; 311 deexcitation->SetVerboseLevel(GetVerboseLevel()); 312 // switch on/off internal electron conversion 313 deexcitation->SetICM(true); 314 // set the maximum life-time for a level that will be treated. Level with life-time longer than this 315 // will be outputed as meta-stable isotope 316 // 317 deexcitation->SetMaxHalfLife(1e-6*second); 318 // but in IT mode, we need to force the transition 319 if (decayMode == 0) { 320 deexcitation->RDMForced(true); 321 } else { 322 deexcitation->RDMForced(false); 323 } 324 // 325 // Get the gammas by deexciting the nucleus. 326 // 327 G4FragmentVector* gammas = deexcitation->BreakItUp(nucleus); 328 // in the case of BreakItUp(nucleus), the returned G4FragmentVector contains the residual nuclide 329 // as its last entry. 330 G4int nGammas=gammas->size()-1; 331 // 332 // Go through each gamma/e- and add it to the decay product. The angular distribution 333 // of the gammas is isotropic, and the residual nucleus is assumed not to have suffered 334 // any recoil as a result of this de-excitation. 335 // 336 for (G4int ig=0; ig<nGammas; ig++) 337 { 338 G4DynamicParticle *theGammaRay = new 339 G4DynamicParticle (gammas->operator[](ig)->GetParticleDefinition(), 340 gammas->operator[](ig)->GetMomentum()); 341 theGammaRay -> SetProperTime(gammas->operator[](ig)->GetCreationTime()); 342 products->PushProducts (theGammaRay); 343 } 344 // 345 // now the nucleus 346 G4double finalDaughterExcitation = gammas->operator[](nGammas)->GetExcitationEnergy(); 347 // f.lei (03/01/03) this is needed to fix the crach in test18 348 if (finalDaughterExcitation <= 1.0*keV) finalDaughterExcitation = 0 ; 349 350 // f.lei (07/03/05) added the delete to fix bug#711 351 if (dynamicDaughter) delete dynamicDaughter; 352 353 G4IonTable *theIonTable = (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable()); 354 dynamicDaughter = new G4DynamicParticle 355 (theIonTable->GetIon(daughterZ,daughterA,finalDaughterExcitation), 356 daughterMomentum1); 357 products->PushProducts (dynamicDaughter); 358 // retrive the ICM shell index 359 shellIndex = deexcitation->GetVacantShellNumber(); 360 361 // 362 // Delete/reset variables associated with the gammas. 363 // 364 while (!gammas->empty()) { 365 delete *(gammas->end()-1); 366 gammas->pop_back(); 367 } 368 // gammas->clearAndDestroy(); 369 delete gammas; 370 delete deexcitation; 371 } 372 // 373 // now we have to take care of the EC product which have to go through the ARM 374 // 375 G4int eShell = -1; 291 376 if (decayMode == 3 || decayMode == 4 || decayMode == 5) { 292 G4int eShell = 0;293 377 switch (decayMode) 294 378 { … … 296 380 // 297 381 { 298 eShell = 1;382 eShell = 0; // --> 0 from 1 (f.lei 30/4/2008) 299 383 } 300 384 break; … … 313 397 case ERROR: 314 398 default: 315 G4cout << " There is an error in decay mode selection! exit RDM now" << G4endl; 316 exit(0); 399 G4Exception("G4NuclearDecayChannel::DecayIt()", "601", 400 FatalException, "Error in decay mode selection"); 317 401 } 402 } 403 // now deal with the IT case where ICM may have been applied 404 // 405 if (decayMode == 0) { 406 eShell = shellIndex; 407 } 408 // now apply ARM if there is a vaccancy 409 // 410 if (eShell != -1) { 318 411 G4int aZ = daughterZ; 319 if (aZ > 5 && aZ < 101) { // only applies to 5< Z <101 320 G4AtomicDeexcitation* atomDeex = new G4AtomicDeexcitation(); 321 //no Auger electron generation 322 // atomDeex->ActivateAugerElectronProduction(0); 323 std::vector<G4DynamicParticle*>* armProducts = atomDeex->GenerateParticles(aZ,eShell); 324 325 // pop up the daughter before insertion 412 if (aZ > 5 && aZ < 100) { // only applies to 5< Z <100 413 // Retrieve the corresponding identifier and binding energy of the selected shell 414 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance(); 415 const G4AtomicShell* shell = transitionManager->Shell(aZ, eShell); 416 G4double bindingEnergy = shell->BindingEnergy(); 417 G4int shellId = shell->ShellId(); 418 419 G4AtomicDeexcitation* atomDeex = new G4AtomicDeexcitation(); 420 //the default is no Auger electron generation. 421 // Switch it on/off here! 422 atomDeex->ActivateAugerElectronProduction(true); 423 std::vector<G4DynamicParticle*>* armProducts = atomDeex->GenerateParticles(aZ,shellId); 424 425 // pop up the daughter before insertion; 426 // f.lei (30/04/2008) check if the total kinetic energy is less than 427 // the shell binding energy; if true add the difference to the daughter to conserve the energy 326 428 dynamicDaughter = products->PopProducts(); 327 for (size_t i = 0; i < armProducts->size(); i++) 429 G4double tARMEnergy = 0.0; 430 for (size_t i = 0; i < armProducts->size(); i++) { 328 431 products->PushProducts ((*armProducts)[i]); 432 tARMEnergy += (*armProducts)[i]->GetKineticEnergy(); 433 } 434 if ((bindingEnergy - tARMEnergy) > 0.1*keV){ 435 G4double dEnergy = dynamicDaughter->GetKineticEnergy() + (bindingEnergy - tARMEnergy); 436 dynamicDaughter->SetKineticEnergy(dEnergy); 437 } 438 products->PushProducts(dynamicDaughter); 439 440 #ifdef G4VERBOSE 441 if (GetVerboseLevel()>0) 442 { 443 G4cout <<"G4NuclearDecayChannel::Selected shell number for ARM = " <<shellId <<G4endl; 444 G4cout <<"G4NuclearDecayChannel::ARM products = " <<armProducts->size()<<G4endl; 445 G4cout <<" The binding energy = " << bindingEnergy << G4endl; 446 G4cout <<" Total ARM particle kinetic energy = " << tARMEnergy << G4endl; 447 } 448 #endif 449 329 450 delete armProducts; 330 451 delete atomDeex; 331 products->PushProducts (dynamicDaughter); 332 } 333 } 334 335 // 336 // If the decay is to an excited state of the daughter nuclide, we need 337 // to apply the photo-evaporation process. 338 // 339 if (daughterExcitation > 0.0) 340 { 341 // 342 // 343 // Pop the daughter nucleus off the product vector - we need to retain 344 // the momentum of this particle. 345 // 346 dynamicDaughter = products->PopProducts(); 347 G4LorentzVector daughterMomentum = dynamicDaughter->Get4Momentum(); 348 G4ThreeVector const daughterMomentum1(static_cast<const G4LorentzVector> (daughterMomentum)); 349 // 350 // 351 // Now define a G4Fragment with the correct A, Z and excitation, and declare and 352 // initialise a G4DiscreteGammaDeexcitation object. 353 // 354 355 // daughterMomentum.setT(daughterMomentum.t()+G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( daughterZ, daughterA )+daughterExcitation); 356 357 // daughterMomentum.setT(daughterMomentum.t()+daughterExcitation); 358 G4Fragment nucleus(daughterA, daughterZ, daughterMomentum); 359 //G4LorentzVector p4(0.,0.,0.,G4NucleiProperties::GetNuclearMass(daughterA,daughterZ) 360 // +daughterExcitation); 361 //G4Fragment nucleus(daughterA, daughterZ, p4); 362 // nucleus.SetExcitationEnergy(daughterExcitation); 363 364 // G4VGammaDeexcitation* deexcitation = new G4DiscreteGammaDeexcitation; 365 G4PhotonEvaporation* deexcitation = new G4PhotonEvaporation; 366 deexcitation->SetVerboseLevel(GetVerboseLevel()); 367 // deexcitation->Initialize(nucleus); 368 deexcitation->SetICM(true); 369 if (decayMode == 0) { 370 deexcitation->RDMForced(true); 371 } else { 372 deexcitation->RDMForced(false); 373 } 374 // ARM in G4 is applied but no auger electrons! 375 deexcitation->SetARM(true); 376 // not applied 377 //deexcitation->SetARM(false); 378 // 379 deexcitation->SetMaxHalfLife(1e-6*second); 380 // 381 // Get the gammas by deexciting the nucleus. 382 // 383 G4FragmentVector* gammas = deexcitation->BreakItUp(nucleus); 384 // in the case of BreakItUp(nucleus), the returned G4FragmentVector contains the residual nuclide 385 // as its last entry. 386 G4int nGammas=gammas->size()-1; 387 // 388 // 389 // Go through each gamma/e- and add it to the decay product. The angular distribution 390 // of the gammas is isotropic, and the residual nucleus is assumed not to suffer 391 // any recoil as a result of this de-excitation. 392 // 393 for (G4int ig=0; ig<nGammas; ig++) 394 { 395 // G4double costheta = 2.0*G4UniformRand() - 1.0; 396 // G4double sintheta = std::sqrt((1.0 - costheta) * (1.0+costheta)); 397 // G4double phi = twopi * G4UniformRand(); 398 // G4ParticleMomentum gDirection 399 // (sintheta*std::cos(phi),sintheta*std::sin(phi),costheta); 400 //G4double gEnergy = gammas->operator[](ig)->GetMomentum().e() 401 // - gammas->operator[](ig)->GetParticleDefinition()->GetPDGMass() ; 402 G4DynamicParticle *theGammaRay = new 403 G4DynamicParticle (gammas->operator[](ig)->GetParticleDefinition(), 404 gammas->operator[](ig)->GetMomentum()); 405 theGammaRay -> SetProperTime(gammas->operator[](ig)->GetCreationTime()); 406 products->PushProducts (theGammaRay); 407 } 408 // 409 // now the nucleus 410 G4double finalDaughterExcitation = gammas->operator[](nGammas)->GetExcitationEnergy(); 411 // f.lei (03/01/03) this is needed to fix the crach in test18 412 if (finalDaughterExcitation <= 1.0*keV) finalDaughterExcitation = 0 ; 413 G4IonTable *theIonTable = (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable()); 414 // f.lei (07/03/05) added the delete to fix bug#711 415 if (dynamicDaughter) delete dynamicDaughter; 416 dynamicDaughter = new G4DynamicParticle 417 (theIonTable->GetIon(daughterZ,daughterA,finalDaughterExcitation), 418 daughterMomentum1); 419 products->PushProducts (dynamicDaughter); 420 // 421 // Delete/reset variables associated with the gammas. 422 // 423 // if (nGammas != 0) gammas->clearAndDestroy(); 424 while (!gammas->empty()) { 425 delete *(gammas->end()-1); 426 gammas->pop_back(); 427 } 428 // gammas->clearAndDestroy(); 429 delete gammas; 430 delete deexcitation; 431 } 452 } 453 } 454 432 455 return products; 433 456 }
Note: See TracChangeset
for help on using the changeset viewer.