ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuonicAtomDecay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MuonicAtomDecay.cc
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"
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"
61 #include "G4RandomDirection.hh"
62 #include "G4PhysicalConstants.hh"
63 #include "G4SystemOfUnits.hh"
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
68  const G4String& name)
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 
87 //{delete theTotalResult;}
88 {}
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
93 {
94  return ( a.GetParticleType() == "MuonicAtom" );
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
99 // void
100 // G4MuonicAtomDecay::PreparePhysicsTable(const G4ParticleDefinition& p)
101 // {
102 // G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this,&p);
103 // }
104 
105 // //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
106 
107 // void G4MuonicAtomDecay::BuildPhysicsTable(const G4ParticleDefinition& p)
108 // {
109 // G4HadronicProcessStore::Instance()->PrintInfo(&p);
110 // }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113 
115  const G4Track& aTrack, G4ForceCondition* condition)
116 {
117  *condition = NotForced;
118  // check if this is the beginning of tracking
121  }
122  return theNumberOfInteractionLengthLeft*GetMeanLifeTime(aTrack, condition);
123 }
124 
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
126 
129 {
130  *condition = NotForced;
131  return DBL_MAX; // this will need to be changed in future; fixme
132 }
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135 
138 {
139  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
140  G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
141  G4MuonicAtom* muatom = static_cast<G4MuonicAtom*>(aParticleDef);
142  G4double meanlife = muatom->GetPDGLifeTime();
143 #ifdef G4VERBOSE
144  if (GetVerboseLevel()>1) {
145  G4cout << "mean life time: "<< meanlife/ns << "[ns]" << G4endl;
146  }
147 #endif
148  return meanlife;
149 }
150 
151 
153  const G4Step&)
154 {
155 
156  // mainly based on G4HadronStoppingProcess & G4Decay
157  // if primary is not Alive then do nothing
158  theTotalResult.Clear(); // G4ParticleChange*
159  theTotalResult.Initialize(aTrack);
161  if(aTrack.GetTrackStatus() != fAlive &&
162  aTrack.GetTrackStatus() != fStopButAlive) {
163  return &theTotalResult;
164  }
165 
166  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
167  const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
168  G4MuonicAtom const* muatom = static_cast<G4MuonicAtom const*>(aParticleDef);
169  G4Ions const* baseion = muatom->GetBaseIon();
170  G4int Z = baseion->GetAtomicNumber();
171  G4double Zd = Z;
172  G4double KEnergy = G4MuonicAtomHelper::GetKShellEnergy(Zd); // fixme check
173  G4HadProjectile thePro(aTrack); // G4HadProjectile, here the muonic atom
174  thePro.SetBoundEnergy(KEnergy);
175 
176  G4ForceCondition* condition = nullptr; // it is unused in the following call anyway
177  G4double meanlife = GetMeanLifeTime(aTrack, condition);
178 
179  G4HadFinalState* result = nullptr; // result before converting to G4VParticleChange*
180  // G4int nSecondaries = 0;
181  // save track time and start from zero time
182  // G4double time0 = aTrack.GetGlobalTime(); FillResult does it
183  // see G4Decay DecayIt
184  // see time0 down below
185  thePro.SetGlobalTime(0.0);
186 
187  // do we need G4double fRemainderLifeTime; ???
188 
189  G4double maDTime = theNumberOfInteractionLengthLeft*meanlife; //fixme check
190 #ifdef G4VERBOSE
191  if (GetVerboseLevel()>1) {
192  G4cout << "G4MuonicAtomDecay::DecayIt time set to: "<< maDTime/ns << "[ns]" << G4endl;
193  }
194 #endif
195 
196  // decide on DIO or Capture
197 
198  G4double lambdac = 1./muatom->GetDIOLifeTime();
199  G4double lambdad = 1./muatom->GetNCLifeTime();
200  G4double lambda = lambdac + lambdad;
201 
202  if ( G4UniformRand()*lambda < lambdac) {
203  // if ( false ) { // force NC for testing
204 
205  // DIO
206  // result = dmptr->ApplyYourself(thePro, *nucleus); // not quite the reaction;
207  // using G4PhaseSpaceDecayChannel
208 
209 #ifdef G4VERBOSE
210  if (GetVerboseLevel()>0) {
211  G4cout << "G4MuonicAtomDecay::DecayIt: selected DIO mode" << G4endl;
212  }
213 #endif
214 
215  // decay table; we use it only for the DIO which is more of a decay
216  // code mostly copied from G4Decay
217 
218  G4DecayProducts* products = nullptr;
219  G4DecayTable *decaytable = aParticleDef->GetDecayTable();
220  G4VDecayChannel* decaychannel = nullptr;
221  G4double massParent = aParticle->GetMass();
222  decaychannel = decaytable->SelectADecayChannel(massParent);
223  if ( decaychannel ==0) {
224  // decay channel not found
226  ed << "Can not determine decay channel for "
227  << aParticleDef->GetParticleName() << G4endl
228  << " mass of dynamic particle: " << massParent/GeV << " (GEV)" << G4endl
229  << " dacay table has " << decaytable->entries() << " entries" << G4endl;
230  G4double checkedmass=massParent;
231  if (massParent < 0.) {
232  checkedmass=aParticleDef->GetPDGMass();
233  ed << "Using PDG mass ("<<checkedmass/GeV << "(GeV)) in IsOKWithParentMass" << G4endl;
234  }
235  for (G4int ic =0;ic <decaytable->entries();++ic) {
236  G4VDecayChannel * dc= decaytable->GetDecayChannel(ic);
237  ed << ic << ": BR " << dc->GetBR() << ", IsOK? "
238  << dc->IsOKWithParentMass(checkedmass)
239  << ", --> ";
240  G4int ndaughters=dc->GetNumberOfDaughters();
241  for (G4int id=0;id<ndaughters;++id) {
242  if (id>0) ed << " + "; // seperator, except for first
243  ed << dc->GetDaughterName(id);
244  }
245  ed << G4endl;
246  }
247  G4Exception("G4MuonicAtomDecay::DecayIt", "DECAY003", FatalException,ed);
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) {
261  ed << "No products are generated for "
262  << aParticleDef->GetParticleName();
263  G4Exception("G4MuonicAtomDecay::DecayIt","DECAY003",FatalException,ed);
264  }
265 #ifdef G4VERBOSE
266  if (GetVerboseLevel()>1) {
267  decaychannel->SetVerboseLevel(temp);
268  }
269 #endif
270 #ifdef G4VERBOSE
271  if (GetVerboseLevel()>2) {
272  if (! products->IsChecked() ) products->DumpInfo();
273  }
274 #endif
275  }
276 
277  // get parent particle information ...................................
278  G4double ParentEnergy = aParticle->GetTotalEnergy();
279  G4double ParentMass = aParticle->GetMass();
280  if (ParentEnergy < ParentMass) {
281  if (GetVerboseLevel()>0) {
282  G4cout << "G4MuonicAtomDecay::DecayIt : Total Energy is less than its mass" << G4endl;
283  G4cout << " Particle: " << aParticle->GetDefinition()->GetParticleName();
284  G4cout << " Energy:" << ParentEnergy/MeV << "[MeV]";
285  G4cout << " Mass:" << ParentMass/MeV << "[MeV]";
286  G4cout << G4endl;
287  }
288  G4Exception( "G4MuonicAtomDecay::DecayIt ",
289  "DECAY102",JustWarning,
290  "Total Energy is less than its mass");
291  ParentEnergy = ParentMass;
292  }
293 
294  G4ThreeVector ParentDirection(aParticle->GetMomentumDirection());
295 
296  //boost all decay products to laboratory frame
297  G4double energyDeposit = 0.0;
298  G4double finalGlobalTime = aTrack.GetGlobalTime();
299  G4double finalLocalTime = aTrack.GetLocalTime();
300  if (aTrack.GetTrackStatus() == fStopButAlive ){
301  // AtRest case
302  finalGlobalTime += maDTime;
303  finalLocalTime += maDTime;
304  energyDeposit += aParticle->GetKineticEnergy();
305  } else {
306  // PostStep case
307  products->Boost( ParentEnergy, ParentDirection);
308  }
309  // G4ParticleChangeForDecay fParticleChangeForDecay; // is it equivalent to G4ParticleChange* theTotalResult;
310 
311  //add products in theTotalResult
312  G4int numberOfSecondaries = products->entries();
313  theTotalResult.SetNumberOfSecondaries(numberOfSecondaries);
314 #ifdef G4VERBOSE
315  if (GetVerboseLevel()>1) {
316  G4cout << "G4MuonicAtomDecay::DecayIt : Decay vertex :";
317  G4cout << " Time: " << finalGlobalTime/ns << "[ns]";
318  G4cout << " X:" << (aTrack.GetPosition()).x() /cm << "[cm]";
319  G4cout << " Y:" << (aTrack.GetPosition()).y() /cm << "[cm]";
320  G4cout << " Z:" << (aTrack.GetPosition()).z() /cm << "[cm]";
321  G4cout << G4endl;
322  G4cout << "G4MuonicAtomDecay::DecayIt : decay products in Lab. Frame" << G4endl;
323  products->DumpInfo();
324  }
325 #endif
326  G4int index;
327  G4ThreeVector currentPosition;
328  const G4TouchableHandle thand = aTrack.GetTouchableHandle();
329  for (index=0; index < numberOfSecondaries; index++)
330  {
331  // get current position of the track
332  currentPosition = aTrack.GetPosition();
333  // create a new track object
334  G4Track* secondary = new G4Track( products->PopProducts(),
335  finalGlobalTime ,
336  currentPosition );
337  // switch on good for tracking flag
338  secondary->SetGoodForTrackingFlag();
339  secondary->SetTouchableHandle(thand);
340  // add the secondary track in the List
341  theTotalResult.AddSecondary(secondary);
342  }
343  delete products;
344 
345  // Kill the parent particle
348  theTotalResult.ProposeLocalTime( finalLocalTime );
349 
350  // Clear NumberOfInteractionLengthLeft
352 
353  return &theTotalResult ;
354 
355  } else { //either or
356 
357  // nuclearCapture
358 
359  // model
360  // need to be able to choose between preco or bertini; no good way to do it?
361  // hardcoded in the constructor for now
362 
363 #ifdef G4VERBOSE
364  if (GetVerboseLevel()>0) {
365  G4cout << "G4MuonicAtomDecay::DecayIt: selected NC mode" << G4endl;
366  }
367 #endif
368 
369  G4int A = baseion->GetAtomicMass();
370  // G4Nucleus* nucleus = GetTargetNucleusPointer(); // from G4HadronicProcess
371  G4Nucleus nucleus;
372  nucleus.SetParameters(A, Z);
373 
374  // we define a local projectile here which will be the orbiting muon
375  // we shall assume it is at rest; fixme
376 
377  // G4HadProjectile, here the muon
379  G4ThreeVector(0.,0.,0.)));
380  theMuPro.SetBoundEnergy(KEnergy);
381  theMuPro.SetGlobalTime(0.0);
382 
383  G4int reentryCount = 0; // this may be in the model already; check fixme <---
384  do {
385  // sample final state
386  // nuclear interaction should keep G4HadFinalState object
387  // model should define time of each secondary particle
388  try {
389  result = cmptr->ApplyYourself(theMuPro, nucleus); // muon and muonic atom nucleus
390  ++reentryCount;
391  }
392  catch(G4HadronicException & aR) {
394  ed << "Call for " << cmptr->GetModelName() << G4endl;
395  ed << " Z= "
396  << nucleus.GetZ_asInt()
397  << " A= " << nucleus.GetA_asInt() << G4endl;
398  DumpState(aTrack,"ApplyYourself",ed);
399  ed << " ApplyYourself failed" << G4endl;
400  G4Exception("G4MuonicAtomDecay::DecayIt", "HAD_MAD_101",
401  FatalException, ed);
402  }
403 
404  // Check the result for catastrophic energy non-conservation
405  // result = CheckResult(theMuPro, nucleus, result);
406 
407  if(reentryCount>100) {
409  ed << "Call for " << cmptr->GetModelName() << G4endl;
410  ed << " Z= "
411  << nucleus.GetZ_asInt()
412  << " A= " << nucleus.GetA_asInt() << G4endl;
413  DumpState(aTrack,"ApplyYourself",ed);
414  ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
415  G4Exception("G4MuonicAtomDecay::DecayIt", "HAD_MAD_102",
416  FatalException, ed);
417  }
418  // Loop checking, 06-Aug-2015, Vladimir Ivanchenko
419  } while(!result);
420 
421  // add delay time of capture (inter + intra)
422  G4int nsec = result->GetNumberOfSecondaries();
423  for(G4int i=0; i<nsec; ++i) {
424  G4HadSecondary* sec = result->GetSecondary(i);
425  G4double ctime = sec->GetTime();
426  sec->SetTime(maDTime + ctime); // we add time0 in the next stage
427 #ifdef G4VERBOSE
428  if (GetVerboseLevel()>1) {
429  G4cout << "G4MuonicAtomDecay::DecayIt time set to: "
430  << (maDTime + ctime)/ns << "[ns]" << G4endl;
431  }
432 #endif
433  }
434 
435  FillResult(result,aTrack);
436 
437  // delete result;// causes bad free check fixme; move to the class members?
438 
440  return &theTotalResult;
441 
442  }
443 
444 
445 }
446 
447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
448 
449 void G4MuonicAtomDecay::ProcessDescription(std::ostream& outFile) const
450 {
451  outFile << "MuonicAtom process where Muon decays in orbit or is captured by the nucleus." <<G4endl;
452 }
453 
454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
455 
457 {
458  // based on G4HadronicProcess::FillResult
459 
461 
462  G4double rotation = CLHEP::twopi*G4UniformRand();
463  G4ThreeVector it(0., 0., 1.);
464 
465  G4double efinal = aR->GetEnergyChange();
466  if(efinal < 0.0) { efinal = 0.0; }
467 
468  // check status of primary
469  if(aR->GetStatusChange() == stopAndKill) {
472 
473  // check its final energy
474  } else if(0.0 == efinal) {
477  ->GetAtRestProcessVector()->size() > 0)
479  else { theTotalResult.ProposeTrackStatus(fStopAndKill); } // check fixme
480 
481  // primary is not killed apply rotation and Lorentz transformation
482  } else {
485  G4double newE = efinal + mass;
486  G4double newP = std::sqrt(efinal*(efinal + 2*mass));
487  G4ThreeVector newPV = newP*aR->GetMomentumChange();
488  G4LorentzVector newP4(newE, newPV);
489  newP4.rotate(rotation, it);
490  newP4 *= aR->GetTrafoToLab();
492  newE = newP4.e() - mass;
493 #ifdef G4VERBOSE
494  if (GetVerboseLevel()>1 && newE <= 0.0) {
496  DumpState(aT,"Primary has zero energy after interaction",ed);
497  G4Exception("G4MuonicAtomDecay::FillResults", "HAD_MAD_103", JustWarning, ed);
498  }
499 #endif
500  if(newE < 0.0) { newE = 0.0; }
502  }
503  //G4cout << "FillResult: Efinal= " << efinal << " status= "
504  // << theTotalResult.GetTrackStatus()
505  // << " fKill= " << fStopAndKill << G4endl;
506 
507  // check secondaries: apply rotation and Lorentz transformation
508  G4int nSec = aR->GetNumberOfSecondaries();
510  G4double weight = aT.GetWeight();
511 
512  if (nSec > 0) {
513  G4double time0 = aT.GetGlobalTime();
514  for (G4int i = 0; i < nSec; ++i) {
516  theM.rotate(rotation, it);
517  theM *= aR->GetTrafoToLab();
518  aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
519 
520  // time of interaction starts from zero
521  G4double time = aR->GetSecondary(i)->GetTime();
522  if (time < 0.0) { time = 0.0; }
523 
524  // take into account global time
525  time += time0;
526 
527  G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
528  time, aT.GetPosition());
530  G4double newWeight = weight*aR->GetSecondary(i)->GetWeight();
531  // G4cout << "#### ParticleDebug "
532  // <<GetProcessName()<<" "
533  //<<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" "
534  // <<aScaleFactor<<" "
535  // <<XBiasSurvivalProbability()<<" "
536  // <<XBiasSecondaryWeight()<<" "
537  // <<aT.GetWeight()<<" "
538  // <<aR->GetSecondary(i)->GetWeight()<<" "
539  // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
540  // <<G4endl;
541  track->SetWeight(newWeight);
544 #ifdef G4VERBOSE
545  if (GetVerboseLevel()>1) {
546  G4double e = track->GetKineticEnergy();
547  if (e <= 0.0) {
549  DumpState(aT,"Secondary has zero energy",ed);
550  ed << "Secondary " << track->GetDefinition()->GetParticleName()
551  << G4endl;
552  G4Exception("G4MuonicAtomDecay::FillResults", "HAD_MAD_103",
553  JustWarning,ed);
554  }
555  }
556 #endif
557  }
558  }
559  aR->Clear();
560 }
561 
562 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
563 
565  const G4String& method,
567 {
568  ed << "Unrecoverable error in the method " << method << " of "
569  << GetProcessName() << G4endl;
570  ed << "TrackID= "<< aTrack.GetTrackID() << " ParentID= "
571  << aTrack.GetParentID()
572  << " " << aTrack.GetParticleDefinition()->GetParticleName()
573  << G4endl;
574  ed << "Ekin(GeV)= " << aTrack.GetKineticEnergy()/CLHEP::GeV
575  << "; direction= " << aTrack.GetMomentumDirection() << G4endl;
576  ed << "Position(mm)= " << aTrack.GetPosition()/CLHEP::mm << ";";
577 
578  if (aTrack.GetMaterial()) {
579  ed << " material " << aTrack.GetMaterial()->GetName();
580  }
581  ed << G4endl;
582 
583  if (aTrack.GetVolume()) {
584  ed << "PhysicalVolume <" << aTrack.GetVolume()->GetName()
585  << ">" << G4endl;
586  }
587 }
588 
589 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
590 
592 {
593  // based on G4Decay::GetMeanFreePath check; fixme
594 
595  // get particle
596  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
597  const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
598  G4double aMass = aParticle->GetMass();
599  G4double aLife = aParticleDef->GetPDGLifeTime();
600 
601  // returns the mean free path in GEANT4 internal units
602  G4double pathlength;
603  G4double aCtau = c_light * aLife;
604 
605  // check if the particle is stable?
606  if (aParticleDef->GetPDGStable()) {
607  pathlength = DBL_MAX;
608 
609  //check if the particle has very short life time ?
610  } else if (aCtau < DBL_MIN) {
611  pathlength = DBL_MIN;
612 
613  } else {
614  //calculate the mean free path
615  // by using normalized kinetic energy (= Ekin/mass)
616  G4double rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
617  const G4double HighestValue = 20.0; //
618  if ( rKineticEnergy > HighestValue) {
619  // gamma >> 1
620  pathlength = ( rKineticEnergy + 1.0)* aCtau;
621  } else if ( rKineticEnergy < DBL_MIN ) {
622  // too slow particle
623 #ifdef G4VERBOSE
624  if (GetVerboseLevel()>1) {
625  G4cout << "G4MuonicAtomDecay::GetMeanFreePath() !!particle stops!!";
626  G4cout << aParticleDef->GetParticleName() << G4endl;
627  G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
628  }
629 #endif
630  pathlength = DBL_MIN;
631  } else {
632  // beta <1
633  pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
634  }
635  }
636  return pathlength;
637 }