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 = (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