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 Class 30 // 31 // File name: G4MuonicAtomDecay 32 // 33 // 20170522 K L Genser first implementation ba 34 // V.Ivantchenko & G4HadronicProcess & G4Decay 35 // 36 // Class Description: 37 // 38 // MuonicAtom Process where Muon either decays 39 // 40 // Modifications: 41 // 42 // 43 //-------------------------------------------- 44 45 #include "G4MuonicAtomDecay.hh" 46 #include "G4HadronicProcessStore.hh" 47 #include "G4HadronicProcessType.hh" 48 #include "G4Nucleus.hh" 49 #include "G4ProcessManager.hh" 50 #include "G4HadFinalState.hh" 51 #include "G4HadProjectile.hh" 52 #include "G4HadSecondary.hh" 53 #include "G4ForceCondition.hh" 54 #include "G4MuonicAtom.hh" 55 #include "G4MuonicAtomHelper.hh" 56 #include "G4VDecayChannel.hh" 57 #include "G4DecayTable.hh" 58 #include "G4DecayProducts.hh" 59 #include "G4CascadeInterface.hh" 60 #include "G4MuMinusCapturePrecompound.hh" 61 #include "G4RandomDirection.hh" 62 #include "G4PhysicalConstants.hh" 63 #include "G4SystemOfUnits.hh" 64 65 //....oooOO0OOooo........oooOO0OOooo........oo 66 67 G4MuonicAtomDecay::G4MuonicAtomDecay(G4Hadroni 68 const G4S 69 : G4VRestDiscreteProcess(name, fDecay), 70 fMuMass(G4MuonMinus::MuonMinus()->GetPDGMa 71 cmptr(hiptr), 72 verboseLevel(0) 73 { 74 // This is not a hadronic process; assume it 75 enableAtRestDoIt = true; 76 enablePostStepDoIt = true; // it is a streac 77 theProcessSubType = 221; // (see G4DecayProc 78 if (!cmptr) { 79 // cmptr = new G4CascadeInterface(); // Be 80 cmptr = new G4MuMinusCapturePrecompound(); 81 } 82 } 83 84 //....oooOO0OOooo........oooOO0OOooo........oo 85 86 G4MuonicAtomDecay::~G4MuonicAtomDecay() 87 {} 88 89 //....oooOO0OOooo........oooOO0OOooo........oo 90 91 G4bool G4MuonicAtomDecay::IsApplicable(const G 92 { 93 return ( a.GetParticleType() == "MuonicAtom" 94 } 95 96 //....oooOO0OOooo........oooOO0OOooo........oo 97 98 // void 99 // G4MuonicAtomDecay::PreparePhysicsTable(cons 100 // { 101 // G4HadronicProcessStore::Instance()->Regis 102 // } 103 104 // //....oooOO0OOooo........oooOO0OOooo....... 105 106 // void G4MuonicAtomDecay::BuildPhysicsTable(c 107 // { 108 // G4HadronicProcessStore::Instance()->Print 109 // } 110 111 //....oooOO0OOooo........oooOO0OOooo........oo 112 113 G4double G4MuonicAtomDecay::AtRestGetPhysicalI 114 const G4Track& aTrack, G4ForceCondition* c 115 { 116 *condition = NotForced; 117 // check if this is the beginning of tracki 118 if (theNumberOfInteractionLengthLeft < 0.) { 119 ResetNumberOfInteractionLengthLeft(); 120 } 121 return theNumberOfInteractionLengthLeft*GetM 122 } 123 124 //....oooOO0OOooo........oooOO0OOooo........oo 125 126 G4double G4MuonicAtomDecay::PostStepGetPhysica 127 const G4Track&, G4double, G4ForceCondition 128 { 129 *condition = NotForced; 130 return DBL_MAX; // this will need to be chan 131 } 132 133 //....oooOO0OOooo........oooOO0OOooo........oo 134 135 G4double G4MuonicAtomDecay::GetMeanLifeTime(co 136 G4 137 { 138 const G4DynamicParticle* aParticle = aTrack 139 G4ParticleDefinition* aParticleDef = aParti 140 G4MuonicAtom* muatom = static_cast<G4Muonic 141 G4double meanlife = muatom->GetPDGLifeTime( 142 #ifdef G4VERBOSE 143 if (GetVerboseLevel()>1) { 144 G4cout << "mean life time: "<< meanlife/n 145 } 146 #endif 147 return meanlife; 148 } 149 150 151 G4VParticleChange* G4MuonicAtomDecay::DecayIt( 152 153 { 154 155 // mainly based on G4HadronStoppingProcess & 156 // if primary is not Alive then do nothing 157 theTotalResult.Clear(); // G4ParticleChang 158 theTotalResult.Initialize(aTrack); 159 theTotalResult.ProposeWeight(aTrack.GetWeigh 160 if(aTrack.GetTrackStatus() != fAlive && 161 aTrack.GetTrackStatus() != fStopButAlive) 162 return &theTotalResult; 163 } 164 165 const G4DynamicParticle* aParticle = aTrack. 166 const G4ParticleDefinition* aParticleDef = a 167 G4MuonicAtom const* muatom = static_cast<G4M 168 G4Ions const* baseion = muatom->GetBaseIon() 169 G4int Z = baseion->GetAtomicNumber(); 170 G4double Zd = Z; 171 G4double KEnergy = G4MuonicAtomHelper::GetKS 172 G4HadProjectile thePro(aTrack); // G4HadProj 173 thePro.SetBoundEnergy(KEnergy); 174 175 G4ForceCondition* condition = nullptr; // it 176 G4double meanlife = GetMeanLifeTime(aTrack, 177 178 G4HadFinalState* result = nullptr; // resul 179 // G4int nSecondaries = 0; 180 // save track time and start from zero time 181 // G4double time0 = aTrack.GetGlobalTime(); 182 // see G4Decay DecayIt 183 // see time0 down below 184 thePro.SetGlobalTime(0.0); 185 186 // do we need G4double fRemainderLifeTime; ? 187 188 G4double maDTime = theNumberOfInteractionLe 189 #ifdef G4VERBOSE 190 if (GetVerboseLevel()>1) { 191 G4cout << "G4MuonicAtomDecay::DecayIt time 192 } 193 #endif 194 195 // decide on DIO or Capture 196 197 G4double lambdac = 1./muatom->GetDIOLifeTime 198 G4double lambdad = 1./muatom->GetNCLifeTime( 199 G4double lambda = lambdac + lambdad; 200 201 if ( G4UniformRand()*lambda < lambdac) { 202 // if ( false ) { // force NC for testing 203 204 // DIO 205 // result = dmptr->ApplyYourself(thePro, * 206 // using G4PhaseSpaceDecayChannel 207 208 #ifdef G4VERBOSE 209 if (GetVerboseLevel()>0) { 210 G4cout << "G4MuonicAtomDecay::DecayIt: selec 211 } 212 #endif 213 214 // decay table; we use it only for the DIO 215 // code mostly copied from G4Decay 216 217 G4DecayProducts* products = nullptr; 218 G4DecayTable *decaytable = aParticleDef- 219 G4VDecayChannel* decaychannel = nullptr; 220 G4double massParent = aParticle->GetMass() 221 decaychannel = decaytable->SelectADecayCha 222 if (decaychannel == nullptr) { 223 // decay channel not found 224 G4ExceptionDescription ed; 225 ed << "Can not determine decay channel f 226 << aParticleDef->GetParticleName() << 227 << " mass of dynamic particle: " << 228 << " dacay table has " << decaytable 229 G4double checkedmass=massParent; 230 if (massParent < 0.) { 231 checkedmass=aParticleDef->GetPDGMass() 232 ed << "Using PDG mass ("<<checkedmass/ 233 } 234 for (G4int ic =0;ic <decaytable->entries 235 G4VDecayChannel * dc= decaytable->GetD 236 ed << ic << ": BR " << dc->GetBR() << 237 << dc->IsOKWithParentMass(checkedma 238 << ", --> "; 239 G4int ndaughters=dc->GetNumberOfDaught 240 for (G4int id=0;id<ndaughters;++id) { 241 if (id>0) ed << " + "; // seperato 242 ed << dc->GetDaughterName(id); 243 } 244 ed << G4endl; 245 } 246 G4Exception("G4MuonicAtomDecay::DecayIt" 247 return &theTotalResult; 248 } else { 249 // execute DecayIt() 250 #ifdef G4VERBOSE 251 G4int temp = decaychannel->GetVerboseLev 252 if (GetVerboseLevel()>1) { 253 G4cout << "G4MuonicAtomDecay::DecayIt : sel 254 << decaychannel <<G4endl; 255 decaychannel->SetVerboseLevel(GetVerboseLeve 256 } 257 #endif 258 products = decaychannel->DecayIt(aPartic 259 if(products == nullptr) { 260 G4ExceptionDescription ed; 261 ed << "No products are generated for " 262 << aParticleDef->GetParticleName(); 263 G4Exception("G4MuonicAtomDecay::DecayIt","DE 264 return &theTotalResult; 265 } 266 #ifdef G4VERBOSE 267 if (GetVerboseLevel()>1) { 268 decaychannel->SetVerboseLevel(temp); 269 } 270 if (GetVerboseLevel()>2) { 271 if (! products->IsChecked() ) products->Dump 272 } 273 #endif 274 } 275 276 // get parent particle information ....... 277 G4double ParentEnergy = aParticle->GetT 278 G4double ParentMass = aParticle->GetM 279 if (ParentEnergy < ParentMass) { 280 if (GetVerboseLevel()>0) { 281 G4cout << "G4MuonicAtomDecay::DecayIt 282 G4cout << " Particle: " << aParticle-> 283 G4cout << " Energy:" << ParentEnerg 284 G4cout << " Mass:" << ParentMass/Me 285 G4cout << G4endl; 286 } 287 G4Exception( "G4MuonicAtomDecay::DecayIt 288 "DECAY102",JustWarning, 289 "Total Energy is less than 290 ParentEnergy = ParentMass; 291 } 292 293 G4ThreeVector ParentDirection(aParticle->G 294 295 //boost all decay products to laboratory f 296 G4double energyDeposit = 0.0; 297 G4double finalGlobalTime = aTrack.GetGloba 298 G4double finalLocalTime = aTrack.GetLocalT 299 if (aTrack.GetTrackStatus() == fStopButAli 300 // AtRest case 301 finalGlobalTime += maDTime; 302 finalLocalTime += maDTime; 303 energyDeposit += aParticle->GetKineticEn 304 } else { 305 // PostStep case 306 products->Boost( ParentEnergy, ParentDir 307 } 308 309 //add products in theTotalResult 310 G4int numberOfSecondaries = products->entr 311 theTotalResult.SetNumberOfSecondaries(numb 312 313 #ifdef G4VERBOSE 314 if (GetVerboseLevel()>1) { 315 G4cout << "G4MuonicAtomDecay::DecayIt : 316 G4cout << " Time: " << finalGlobalTime/n 317 G4cout << " X:" << (aTrack.GetPosition() 318 G4cout << " Y:" << (aTrack.GetPosition() 319 G4cout << " Z:" << (aTrack.GetPosition() 320 G4cout << G4endl; 321 G4cout << "G4MuonicAtomDecay::DecayIt : 322 products->DumpInfo(); 323 } 324 #endif 325 G4int index; 326 G4ThreeVector currentPosition; 327 const G4TouchableHandle thand = aTrack.Get 328 for (index=0; index < numberOfSecondaries; 329 { 330 // get current position of the track 331 currentPosition = aTrack.GetPosition() 332 // create a new track object 333 G4Track* secondary = new G4Track( prod 334 fina 335 curr 336 // switch on good for tracking flag 337 secondary->SetGoodForTrackingFlag(); 338 secondary->SetTouchableHandle(thand); 339 // add the secondary track in the List 340 theTotalResult.AddSecondary(secondary) 341 } 342 delete products; 343 344 // Kill the parent particle 345 theTotalResult.ProposeTrackStatus( fStopAn 346 theTotalResult.ProposeLocalEnergyDeposit(e 347 theTotalResult.ProposeLocalTime( finalLoca 348 349 // Clear NumberOfInteractionLengthLeft 350 ClearNumberOfInteractionLengthLeft(); 351 352 return &theTotalResult ; 353 354 } else { //either or 355 356 // nuclearCapture 357 358 // model 359 // need to be able to choose between preco 360 // hardcoded in the constructor for now 361 362 #ifdef G4VERBOSE 363 if (GetVerboseLevel()>0) { 364 G4cout << "G4MuonicAtomDecay::DecayIt: selec 365 } 366 #endif 367 368 G4int A = baseion->GetAtomicMass(); 369 // G4Nucleus* nucleus = GetTargetNucleu 370 G4Nucleus nucleus; 371 nucleus.SetParameters(A, Z); 372 373 // we define a local projectile here which 374 // we shall assume it is at rest; fixme 375 376 // G4HadProjectile, here the muon 377 G4HadProjectile theMuPro(G4DynamicParticle 378 379 theMuPro.SetBoundEnergy(KEnergy); 380 theMuPro.SetGlobalTime(0.0); 381 382 G4int reentryCount = 0; // this may be in 383 do { 384 // sample final state 385 // nuclear interaction should keep G4Had 386 // model should define time of each seco 387 try { 388 result = cmptr->ApplyYourself(theMuPro 389 ++reentryCount; 390 } 391 catch(G4HadronicException & aR) { 392 G4ExceptionDescription ed; 393 ed << "Call for " << cmptr->GetModelNa 394 ed << " Z= " 395 << nucleus.GetZ_asInt() 396 << " A= " << nucleus.GetA_asInt() 397 DumpState(aTrack,"ApplyYourself",ed); 398 ed << " ApplyYourself failed" << G4end 399 G4Exception("G4MuonicAtomDecay::DecayI 400 FatalException, ed); 401 } 402 403 // Check the result for catastrophic ene 404 // result = CheckResult(theMuPro, nucleu 405 406 if(reentryCount>100) { 407 G4ExceptionDescription ed; 408 ed << "Call for " << cmptr->GetModelNa 409 ed << " Z= " 410 << nucleus.GetZ_asInt() 411 << " A= " << nucleus.GetA_asInt() 412 DumpState(aTrack,"ApplyYourself",ed); 413 ed << " ApplyYourself does not complet 414 G4Exception("G4MuonicAtomDecay::DecayI 415 FatalException, ed); 416 } 417 // Loop checking, 06-Aug-2015, Vladimir 418 } while(result == nullptr); 419 420 // add delay time of capture (inter + intr 421 G4int nsec = (G4int)result->GetNumberOfSec 422 for(G4int i=0; i<nsec; ++i) { 423 G4HadSecondary* sec = result->GetSeconda 424 G4double ctime = sec->GetTime(); 425 sec->SetTime(maDTime + ctime); // we add 426 #ifdef G4VERBOSE 427 if (GetVerboseLevel()>1) { 428 G4cout << "G4MuonicAtomDecay::DecayIt 429 << (maDTime + ctime)/ns << "[ns 430 } 431 #endif 432 } 433 434 FillResult(result,aTrack); 435 436 ClearNumberOfInteractionLengthLeft(); 437 return &theTotalResult; 438 } 439 } 440 441 //....oooOO0OOooo........oooOO0OOooo........oo 442 443 void G4MuonicAtomDecay::ProcessDescription(std 444 { 445 outFile << "MuonicAtom process where Muon de 446 } 447 448 //....oooOO0OOooo........oooOO0OOooo........oo 449 450 void G4MuonicAtomDecay::FillResult(G4HadFinalS 451 { 452 // based on G4HadronicProcess::FillResult 453 454 theTotalResult.ProposeLocalEnergyDeposit(aR- 455 456 G4double rotation = CLHEP::twopi*G4UniformRa 457 G4ThreeVector it(0., 0., 1.); 458 459 G4double efinal = aR->GetEnergyChange(); 460 if(efinal < 0.0) { efinal = 0.0; } 461 462 // check status of primary 463 if(aR->GetStatusChange() == stopAndKill) { 464 theTotalResult.ProposeTrackStatus(fStopAnd 465 theTotalResult.ProposeEnergy( 0.0 ); 466 467 // check its final energy 468 } else if(0.0 == efinal) { 469 theTotalResult.ProposeEnergy( 0.0 ); 470 if(aT.GetParticleDefinition()->GetProcessM 471 ->GetAtRestProcessVector()->size() > 0) 472 { theTotalResult.ProposeTrackStatus(f 473 else { theTotalResult.ProposeTrackStatus(f 474 475 // primary is not killed apply rotation an 476 } else { 477 theTotalResult.ProposeTrackStatus(fAlive); 478 G4double mass = aT.GetParticleDefinition() 479 G4double newE = efinal + mass; 480 G4double newP = std::sqrt(efinal*(efinal + 481 G4ThreeVector newPV = newP*aR->GetMomentum 482 G4LorentzVector newP4(newE, newPV); 483 newP4.rotate(rotation, it); 484 newP4 *= aR->GetTrafoToLab(); 485 theTotalResult.ProposeMomentumDirection(ne 486 newE = newP4.e() - mass; 487 #ifdef G4VERBOSE 488 if (GetVerboseLevel()>1 && newE <= 0.0) { 489 G4ExceptionDescription ed; 490 DumpState(aT,"Primary has zero energy af 491 G4Exception("G4MuonicAtomDecay::FillResu 492 } 493 #endif 494 if(newE < 0.0) { newE = 0.0; } 495 theTotalResult.ProposeEnergy( newE ); 496 } 497 //G4cout << "FillResult: Efinal= " << efinal 498 // << theTotalResult.GetTrackStatus() 499 // << " fKill= " << fStopAndKill << G4end 500 501 // check secondaries: apply rotation and Lor 502 G4int nSec = (G4int)aR->GetNumberOfSecondari 503 theTotalResult.SetNumberOfSecondaries(nSec); 504 G4double weight = aT.GetWeight(); 505 506 if (nSec > 0) { 507 G4double time0 = aT.GetGlobalTime(); 508 for (G4int i = 0; i < nSec; ++i) { 509 G4LorentzVector theM = aR->GetSecondary( 510 theM.rotate(rotation, it); 511 theM *= aR->GetTrafoToLab(); 512 aR->GetSecondary(i)->GetParticle()->Set4 513 514 // time of interaction starts from zero 515 G4double time = aR->GetSecondary(i)->Get 516 if (time < 0.0) { time = 0.0; } 517 518 // take into account global time 519 time += time0; 520 521 G4Track* track = new G4Track(aR->GetSeco 522 time, aT.Ge 523 track->SetCreatorModelID(aR->GetSecondar 524 G4double newWeight = weight*aR->GetSecon 525 // G4cout << "#### ParticleDebug " 526 // <<GetProcessName()<<" " 527 //<<aR->GetSecondary(i)->GetParticle()->GetD 528 // <<aScaleFactor<<" " 529 // <<XBiasSurvivalProbability()<<" " 530 // <<XBiasSecondaryWeight()<<" " 531 // <<aT.GetWeight()<<" " 532 // <<aR->GetSecondary(i)->GetWeight()<<" " 533 // <<aR->GetSecondary(i)->GetParticle()->Get 534 // <<G4endl; 535 track->SetWeight(newWeight); 536 track->SetTouchableHandle(aT.GetTouchabl 537 theTotalResult.AddSecondary(track); 538 #ifdef G4VERBOSE 539 if (GetVerboseLevel()>1) { 540 G4double e = track->GetKineticEnergy() 541 if (e <= 0.0) { 542 G4ExceptionDescription ed; 543 DumpState(aT,"Secondary has zero ene 544 ed << "Secondary " << track->GetDefi 545 << G4endl; 546 G4Exception("G4MuonicAtomDecay::Fill 547 JustWarning,ed); 548 } 549 } 550 #endif 551 } 552 } 553 aR->Clear(); 554 } 555 556 //....oooOO0OOooo........oooOO0OOooo........oo 557 558 void G4MuonicAtomDecay::DumpState(const G4Trac 559 const G4String& method, 560 G4ExceptionDescription& ed) 561 { 562 ed << "Unrecoverable error in the method " < 563 << GetProcessName() << G4endl; 564 ed << "TrackID= "<< aTrack.GetTrackID() << " 565 << aTrack.GetParentID() 566 << " " << aTrack.GetParticleDefinition() 567 << G4endl; 568 ed << "Ekin(GeV)= " << aTrack.GetKineticEner 569 << "; direction= " << aTrack.GetMomentum 570 ed << "Position(mm)= " << aTrack.GetPosition 571 572 if (aTrack.GetMaterial()) { 573 ed << " material " << aTrack.GetMaterial( 574 } 575 ed << G4endl; 576 577 if (aTrack.GetVolume()) { 578 ed << "PhysicalVolume <" << aTrack.GetVol 579 << ">" << G4endl; 580 } 581 } 582 583 //....oooOO0OOooo........oooOO0OOooo........oo 584 585 G4double G4MuonicAtomDecay::GetMeanFreePath(co 586 { 587 // based on G4Decay::GetMeanFreePath check; 588 589 // get particle 590 const G4DynamicParticle* aParticle = aTrack. 591 const G4ParticleDefinition* aParticleDef = a 592 G4double aMass = aParticle->GetMass(); 593 G4double aLife = aParticleDef->GetPDGLifeTim 594 595 // returns the mean free path in GEANT4 inte 596 G4double pathlength; 597 G4double aCtau = c_light * aLife; 598 599 // check if the particle is stable? 600 if (aParticleDef->GetPDGStable()) { 601 pathlength = DBL_MAX; 602 603 //check if the particle has very short lif 604 } else if (aCtau < DBL_MIN) { 605 pathlength = DBL_MIN; 606 607 } else { 608 //calculate the mean free path 609 // by using normalized kinetic energy (= E 610 G4double rKineticEnergy = aParticle->Get 611 const G4double HighestValue = 20.0; // 612 if ( rKineticEnergy > HighestValue) { 613 // gamma >> 1 614 pathlength = ( rKineticEnergy + 1.0)* aC 615 } else if ( rKineticEnergy < DBL_MIN ) { 616 // too slow particle 617 #ifdef G4VERBOSE 618 if (GetVerboseLevel()>1) { 619 G4cout << "G4MuonicAtomDecay::GetMeanF 620 G4cout << aParticleDef->GetParticleNam 621 G4cout << "KineticEnergy:" << aParticl 622 } 623 #endif 624 pathlength = DBL_MIN; 625 } else { 626 // beta <1 627 pathlength = (aParticle->GetTotalMomentu 628 } 629 } 630 return pathlength; 631 } 632