| Geant4 Cross Reference |
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 = 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 = 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