| 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