Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/stopping/src/G4MuonicAtomDecay.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/stopping/src/G4MuonicAtomDecay.cc (Version 11.3.0) and /processes/hadronic/stopping/src/G4MuonicAtomDecay.cc (Version 9.1.p1)


  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