| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //---------------------------------------------------------------------
28 //
29 // GEANT4 Class
30 //
31 // File name: G4MuonicAtomDecay
32 //
33 // 20170522 K L Genser first implementation based on code by
34 // V.Ivantchenko & G4HadronicProcess & G4Decay
35 //
36 // Class Description:
37 //
38 // MuonicAtom Process where Muon either decays in orbit or is captured by the nucleus
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........oooOO0OOooo........oooOO0OOooo....
66
67 G4MuonicAtomDecay::G4MuonicAtomDecay(G4HadronicInteraction* hiptr,
68 const G4String& name)
69 : G4VRestDiscreteProcess(name, fDecay),
70 fMuMass(G4MuonMinus::MuonMinus()->GetPDGMass()),
71 cmptr(hiptr),
72 verboseLevel(0)
73 {
74 // This is not a hadronic process; assume it is a kind of decay
75 enableAtRestDoIt = true;
76 enablePostStepDoIt = true; // it is a streach; fixme
77 theProcessSubType = 221; // (see G4DecayProcessType.hh) fixme
78 if (!cmptr) {
79 // cmptr = new G4CascadeInterface(); // Bertini - Pointer owned by InteractionRegistry
80 cmptr = new G4MuMinusCapturePrecompound(); // Precompound
81 }
82 }
83
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85
86 G4MuonicAtomDecay::~G4MuonicAtomDecay()
87 {}
88
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90
91 G4bool G4MuonicAtomDecay::IsApplicable(const G4ParticleDefinition& a)
92 {
93 return ( a.GetParticleType() == "MuonicAtom" );
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97
98 // void
99 // G4MuonicAtomDecay::PreparePhysicsTable(const G4ParticleDefinition& p)
100 // {
101 // G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this,&p);
102 // }
103
104 // //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105
106 // void G4MuonicAtomDecay::BuildPhysicsTable(const G4ParticleDefinition& p)
107 // {
108 // G4HadronicProcessStore::Instance()->PrintInfo(&p);
109 // }
110
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112
113 G4double G4MuonicAtomDecay::AtRestGetPhysicalInteractionLength(
114 const G4Track& aTrack, G4ForceCondition* condition)
115 {
116 *condition = NotForced;
117 // check if this is the beginning of tracking
118 if (theNumberOfInteractionLengthLeft < 0.) {
119 ResetNumberOfInteractionLengthLeft();
120 }
121 return theNumberOfInteractionLengthLeft*GetMeanLifeTime(aTrack, condition);
122 }
123
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125
126 G4double G4MuonicAtomDecay::PostStepGetPhysicalInteractionLength(
127 const G4Track&, G4double, G4ForceCondition* condition)
128 {
129 *condition = NotForced;
130 return DBL_MAX; // this will need to be changed in future; fixme
131 }
132
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134
135 G4double G4MuonicAtomDecay::GetMeanLifeTime(const G4Track& aTrack,
136 G4ForceCondition*)
137 {
138 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
139 G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
140 G4MuonicAtom* muatom = static_cast<G4MuonicAtom*>(aParticleDef);
141 G4double meanlife = muatom->GetPDGLifeTime();
142 #ifdef G4VERBOSE
143 if (GetVerboseLevel()>1) {
144 G4cout << "mean life time: "<< meanlife/ns << "[ns]" << G4endl;
145 }
146 #endif
147 return meanlife;
148 }
149
150
151 G4VParticleChange* G4MuonicAtomDecay::DecayIt(const G4Track& aTrack,
152 const G4Step&)
153 {
154
155 // mainly based on G4HadronStoppingProcess & G4Decay
156 // if primary is not Alive then do nothing
157 theTotalResult.Clear(); // G4ParticleChange*
158 theTotalResult.Initialize(aTrack);
159 theTotalResult.ProposeWeight(aTrack.GetWeight());
160 if(aTrack.GetTrackStatus() != fAlive &&
161 aTrack.GetTrackStatus() != fStopButAlive) {
162 return &theTotalResult;
163 }
164
165 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
166 const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
167 G4MuonicAtom const* muatom = static_cast<G4MuonicAtom const*>(aParticleDef);
168 G4Ions const* baseion = muatom->GetBaseIon();
169 G4int Z = baseion->GetAtomicNumber();
170 G4double Zd = Z;
171 G4double KEnergy = G4MuonicAtomHelper::GetKShellEnergy(Zd); // fixme check
172 G4HadProjectile thePro(aTrack); // G4HadProjectile, here the muonic atom
173 thePro.SetBoundEnergy(KEnergy);
174
175 G4ForceCondition* condition = nullptr; // it is unused in the following call anyway
176 G4double meanlife = GetMeanLifeTime(aTrack, condition);
177
178 G4HadFinalState* result = nullptr; // result before converting to G4VParticleChange*
179 // G4int nSecondaries = 0;
180 // save track time and start from zero time
181 // G4double time0 = aTrack.GetGlobalTime(); FillResult does it
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 = theNumberOfInteractionLengthLeft*meanlife; //fixme check
189 #ifdef G4VERBOSE
190 if (GetVerboseLevel()>1) {
191 G4cout << "G4MuonicAtomDecay::DecayIt time set to: "<< maDTime/ns << "[ns]" << G4endl;
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, *nucleus); // not quite the reaction;
206 // using G4PhaseSpaceDecayChannel
207
208 #ifdef G4VERBOSE
209 if (GetVerboseLevel()>0) {
210 G4cout << "G4MuonicAtomDecay::DecayIt: selected DIO mode" << G4endl;
211 }
212 #endif
213
214 // decay table; we use it only for the DIO which is more of a decay
215 // code mostly copied from G4Decay
216
217 G4DecayProducts* products = nullptr;
218 G4DecayTable *decaytable = aParticleDef->GetDecayTable();
219 G4VDecayChannel* decaychannel = nullptr;
220 G4double massParent = aParticle->GetMass();
221 decaychannel = decaytable->SelectADecayChannel(massParent);
222 if (decaychannel == nullptr) {
223 // decay channel not found
224 G4ExceptionDescription ed;
225 ed << "Can not determine decay channel for "
226 << aParticleDef->GetParticleName() << G4endl
227 << " mass of dynamic particle: " << massParent/GeV << " (GEV)" << G4endl
228 << " dacay table has " << decaytable->entries() << " entries" << G4endl;
229 G4double checkedmass=massParent;
230 if (massParent < 0.) {
231 checkedmass=aParticleDef->GetPDGMass();
232 ed << "Using PDG mass ("<<checkedmass/GeV << "(GeV)) in IsOKWithParentMass" << G4endl;
233 }
234 for (G4int ic =0;ic <decaytable->entries();++ic) {
235 G4VDecayChannel * dc= decaytable->GetDecayChannel(ic);
236 ed << ic << ": BR " << dc->GetBR() << ", IsOK? "
237 << dc->IsOKWithParentMass(checkedmass)
238 << ", --> ";
239 G4int ndaughters=dc->GetNumberOfDaughters();
240 for (G4int id=0;id<ndaughters;++id) {
241 if (id>0) ed << " + "; // seperator, except for first
242 ed << dc->GetDaughterName(id);
243 }
244 ed << G4endl;
245 }
246 G4Exception("G4MuonicAtomDecay::DecayIt", "DECAY003", FatalException,ed);
247 return &theTotalResult;
248 } else {
249 // execute DecayIt()
250 #ifdef G4VERBOSE
251 G4int temp = decaychannel->GetVerboseLevel();
252 if (GetVerboseLevel()>1) {
253 G4cout << "G4MuonicAtomDecay::DecayIt : selected decay channel addr:"
254 << decaychannel <<G4endl;
255 decaychannel->SetVerboseLevel(GetVerboseLevel());
256 }
257 #endif
258 products = decaychannel->DecayIt(aParticle->GetMass());
259 if(products == nullptr) {
260 G4ExceptionDescription ed;
261 ed << "No products are generated for "
262 << aParticleDef->GetParticleName();
263 G4Exception("G4MuonicAtomDecay::DecayIt","DECAY003",FatalException,ed);
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->DumpInfo();
272 }
273 #endif
274 }
275
276 // get parent particle information ...................................
277 G4double ParentEnergy = aParticle->GetTotalEnergy();
278 G4double ParentMass = aParticle->GetMass();
279 if (ParentEnergy < ParentMass) {
280 if (GetVerboseLevel()>0) {
281 G4cout << "G4MuonicAtomDecay::DecayIt : Total Energy is less than its mass" << G4endl;
282 G4cout << " Particle: " << aParticle->GetDefinition()->GetParticleName();
283 G4cout << " Energy:" << ParentEnergy/MeV << "[MeV]";
284 G4cout << " Mass:" << ParentMass/MeV << "[MeV]";
285 G4cout << G4endl;
286 }
287 G4Exception( "G4MuonicAtomDecay::DecayIt ",
288 "DECAY102",JustWarning,
289 "Total Energy is less than its mass");
290 ParentEnergy = ParentMass;
291 }
292
293 G4ThreeVector ParentDirection(aParticle->GetMomentumDirection());
294
295 //boost all decay products to laboratory frame
296 G4double energyDeposit = 0.0;
297 G4double finalGlobalTime = aTrack.GetGlobalTime();
298 G4double finalLocalTime = aTrack.GetLocalTime();
299 if (aTrack.GetTrackStatus() == fStopButAlive ){
300 // AtRest case
301 finalGlobalTime += maDTime;
302 finalLocalTime += maDTime;
303 energyDeposit += aParticle->GetKineticEnergy();
304 } else {
305 // PostStep case
306 products->Boost( ParentEnergy, ParentDirection);
307 }
308
309 //add products in theTotalResult
310 G4int numberOfSecondaries = products->entries();
311 theTotalResult.SetNumberOfSecondaries(numberOfSecondaries);
312
313 #ifdef G4VERBOSE
314 if (GetVerboseLevel()>1) {
315 G4cout << "G4MuonicAtomDecay::DecayIt : Decay vertex :";
316 G4cout << " Time: " << finalGlobalTime/ns << "[ns]";
317 G4cout << " X:" << (aTrack.GetPosition()).x() /cm << "[cm]";
318 G4cout << " Y:" << (aTrack.GetPosition()).y() /cm << "[cm]";
319 G4cout << " Z:" << (aTrack.GetPosition()).z() /cm << "[cm]";
320 G4cout << G4endl;
321 G4cout << "G4MuonicAtomDecay::DecayIt : decay products in Lab. Frame" << G4endl;
322 products->DumpInfo();
323 }
324 #endif
325 G4int index;
326 G4ThreeVector currentPosition;
327 const G4TouchableHandle thand = aTrack.GetTouchableHandle();
328 for (index=0; index < numberOfSecondaries; index++)
329 {
330 // get current position of the track
331 currentPosition = aTrack.GetPosition();
332 // create a new track object
333 G4Track* secondary = new G4Track( products->PopProducts(),
334 finalGlobalTime ,
335 currentPosition );
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( fStopAndKill ) ;
346 theTotalResult.ProposeLocalEnergyDeposit(energyDeposit);
347 theTotalResult.ProposeLocalTime( finalLocalTime );
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 or bertini; no good way to do it?
360 // hardcoded in the constructor for now
361
362 #ifdef G4VERBOSE
363 if (GetVerboseLevel()>0) {
364 G4cout << "G4MuonicAtomDecay::DecayIt: selected NC mode" << G4endl;
365 }
366 #endif
367
368 G4int A = baseion->GetAtomicMass();
369 // G4Nucleus* nucleus = GetTargetNucleusPointer(); // from G4HadronicProcess
370 G4Nucleus nucleus;
371 nucleus.SetParameters(A, Z);
372
373 // we define a local projectile here which will be the orbiting muon
374 // we shall assume it is at rest; fixme
375
376 // G4HadProjectile, here the muon
377 G4HadProjectile theMuPro(G4DynamicParticle(G4MuonMinus::MuonMinus(),
378 G4ThreeVector(0.,0.,0.)));
379 theMuPro.SetBoundEnergy(KEnergy);
380 theMuPro.SetGlobalTime(0.0);
381
382 G4int reentryCount = 0; // this may be in the model already; check fixme <---
383 do {
384 // sample final state
385 // nuclear interaction should keep G4HadFinalState object
386 // model should define time of each secondary particle
387 try {
388 result = cmptr->ApplyYourself(theMuPro, nucleus); // muon and muonic atom nucleus
389 ++reentryCount;
390 }
391 catch(G4HadronicException & aR) {
392 G4ExceptionDescription ed;
393 ed << "Call for " << cmptr->GetModelName() << G4endl;
394 ed << " Z= "
395 << nucleus.GetZ_asInt()
396 << " A= " << nucleus.GetA_asInt() << G4endl;
397 DumpState(aTrack,"ApplyYourself",ed);
398 ed << " ApplyYourself failed" << G4endl;
399 G4Exception("G4MuonicAtomDecay::DecayIt", "HAD_MAD_101",
400 FatalException, ed);
401 }
402
403 // Check the result for catastrophic energy non-conservation
404 // result = CheckResult(theMuPro, nucleus, result);
405
406 if(reentryCount>100) {
407 G4ExceptionDescription ed;
408 ed << "Call for " << cmptr->GetModelName() << G4endl;
409 ed << " Z= "
410 << nucleus.GetZ_asInt()
411 << " A= " << nucleus.GetA_asInt() << G4endl;
412 DumpState(aTrack,"ApplyYourself",ed);
413 ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
414 G4Exception("G4MuonicAtomDecay::DecayIt", "HAD_MAD_102",
415 FatalException, ed);
416 }
417 // Loop checking, 06-Aug-2015, Vladimir Ivanchenko
418 } while(result == nullptr);
419
420 // add delay time of capture (inter + intra)
421 G4int nsec = (G4int)result->GetNumberOfSecondaries();
422 for(G4int i=0; i<nsec; ++i) {
423 G4HadSecondary* sec = result->GetSecondary(i);
424 G4double ctime = sec->GetTime();
425 sec->SetTime(maDTime + ctime); // we add time0 in the next stage
426 #ifdef G4VERBOSE
427 if (GetVerboseLevel()>1) {
428 G4cout << "G4MuonicAtomDecay::DecayIt time set to: "
429 << (maDTime + ctime)/ns << "[ns]" << G4endl;
430 }
431 #endif
432 }
433
434 FillResult(result,aTrack);
435
436 ClearNumberOfInteractionLengthLeft();
437 return &theTotalResult;
438 }
439 }
440
441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
442
443 void G4MuonicAtomDecay::ProcessDescription(std::ostream& outFile) const
444 {
445 outFile << "MuonicAtom process where Muon decays in orbit or is captured by the nucleus." <<G4endl;
446 }
447
448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
449
450 void G4MuonicAtomDecay::FillResult(G4HadFinalState * aR, const G4Track & aT)
451 {
452 // based on G4HadronicProcess::FillResult
453
454 theTotalResult.ProposeLocalEnergyDeposit(aR->GetLocalEnergyDeposit());
455
456 G4double rotation = CLHEP::twopi*G4UniformRand();
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(fStopAndKill);
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()->GetProcessManager()
471 ->GetAtRestProcessVector()->size() > 0)
472 { theTotalResult.ProposeTrackStatus(fStopButAlive); }
473 else { theTotalResult.ProposeTrackStatus(fStopAndKill); } // check fixme
474
475 // primary is not killed apply rotation and Lorentz transformation
476 } else {
477 theTotalResult.ProposeTrackStatus(fAlive);
478 G4double mass = aT.GetParticleDefinition()->GetPDGMass();
479 G4double newE = efinal + mass;
480 G4double newP = std::sqrt(efinal*(efinal + 2*mass));
481 G4ThreeVector newPV = newP*aR->GetMomentumChange();
482 G4LorentzVector newP4(newE, newPV);
483 newP4.rotate(rotation, it);
484 newP4 *= aR->GetTrafoToLab();
485 theTotalResult.ProposeMomentumDirection(newP4.vect().unit());
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 after interaction",ed);
491 G4Exception("G4MuonicAtomDecay::FillResults", "HAD_MAD_103", JustWarning, ed);
492 }
493 #endif
494 if(newE < 0.0) { newE = 0.0; }
495 theTotalResult.ProposeEnergy( newE );
496 }
497 //G4cout << "FillResult: Efinal= " << efinal << " status= "
498 // << theTotalResult.GetTrackStatus()
499 // << " fKill= " << fStopAndKill << G4endl;
500
501 // check secondaries: apply rotation and Lorentz transformation
502 G4int nSec = (G4int)aR->GetNumberOfSecondaries();
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(i)->GetParticle()->Get4Momentum();
510 theM.rotate(rotation, it);
511 theM *= aR->GetTrafoToLab();
512 aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
513
514 // time of interaction starts from zero
515 G4double time = aR->GetSecondary(i)->GetTime();
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->GetSecondary(i)->GetParticle(),
522 time, aT.GetPosition());
523 track->SetCreatorModelID(aR->GetSecondary(i)->GetCreatorModelID());
524 G4double newWeight = weight*aR->GetSecondary(i)->GetWeight();
525 // G4cout << "#### ParticleDebug "
526 // <<GetProcessName()<<" "
527 //<<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" "
528 // <<aScaleFactor<<" "
529 // <<XBiasSurvivalProbability()<<" "
530 // <<XBiasSecondaryWeight()<<" "
531 // <<aT.GetWeight()<<" "
532 // <<aR->GetSecondary(i)->GetWeight()<<" "
533 // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
534 // <<G4endl;
535 track->SetWeight(newWeight);
536 track->SetTouchableHandle(aT.GetTouchableHandle());
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 energy",ed);
544 ed << "Secondary " << track->GetDefinition()->GetParticleName()
545 << G4endl;
546 G4Exception("G4MuonicAtomDecay::FillResults", "HAD_MAD_103",
547 JustWarning,ed);
548 }
549 }
550 #endif
551 }
552 }
553 aR->Clear();
554 }
555
556 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
557
558 void G4MuonicAtomDecay::DumpState(const G4Track& aTrack,
559 const G4String& method,
560 G4ExceptionDescription& ed)
561 {
562 ed << "Unrecoverable error in the method " << method << " of "
563 << GetProcessName() << G4endl;
564 ed << "TrackID= "<< aTrack.GetTrackID() << " ParentID= "
565 << aTrack.GetParentID()
566 << " " << aTrack.GetParticleDefinition()->GetParticleName()
567 << G4endl;
568 ed << "Ekin(GeV)= " << aTrack.GetKineticEnergy()/CLHEP::GeV
569 << "; direction= " << aTrack.GetMomentumDirection() << G4endl;
570 ed << "Position(mm)= " << aTrack.GetPosition()/CLHEP::mm << ";";
571
572 if (aTrack.GetMaterial()) {
573 ed << " material " << aTrack.GetMaterial()->GetName();
574 }
575 ed << G4endl;
576
577 if (aTrack.GetVolume()) {
578 ed << "PhysicalVolume <" << aTrack.GetVolume()->GetName()
579 << ">" << G4endl;
580 }
581 }
582
583 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
584
585 G4double G4MuonicAtomDecay::GetMeanFreePath(const G4Track& aTrack,G4double, G4ForceCondition*)
586 {
587 // based on G4Decay::GetMeanFreePath check; fixme
588
589 // get particle
590 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
591 const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
592 G4double aMass = aParticle->GetMass();
593 G4double aLife = aParticleDef->GetPDGLifeTime();
594
595 // returns the mean free path in GEANT4 internal units
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 life time ?
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 (= Ekin/mass)
610 G4double rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
611 const G4double HighestValue = 20.0; //
612 if ( rKineticEnergy > HighestValue) {
613 // gamma >> 1
614 pathlength = ( rKineticEnergy + 1.0)* aCtau;
615 } else if ( rKineticEnergy < DBL_MIN ) {
616 // too slow particle
617 #ifdef G4VERBOSE
618 if (GetVerboseLevel()>1) {
619 G4cout << "G4MuonicAtomDecay::GetMeanFreePath() !!particle stops!!";
620 G4cout << aParticleDef->GetParticleName() << G4endl;
621 G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
622 }
623 #endif
624 pathlength = DBL_MIN;
625 } else {
626 // beta <1
627 pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
628 }
629 }
630 return pathlength;
631 }
632