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 11.1)


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