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 // $Id: G4EnergyLossTables.cc 77290 2013-11-22 10:57:50Z gcosmo $ 27 // 28 // 28 // ------------------------------------------- 29 // ------------------------------------------------------------------- 29 // first version created by P.Urban , 06/04/19 30 // first version created by P.Urban , 06/04/1998 30 // modifications + "precise" functions added b 31 // modifications + "precise" functions added by L.Urban , 27/05/98 31 // modifications , TOF functions , 26/10/98, L 32 // modifications , TOF functions , 26/10/98, L.Urban 32 // cache mechanism in order to gain time, 11/0 33 // cache mechanism in order to gain time, 11/02/99, L.Urban 33 // bug fixed , 12/04/99 , L.Urban 34 // bug fixed , 12/04/99 , L.Urban 34 // 10.11.99: moved from RWT hash dictionary to 35 // 10.11.99: moved from RWT hash dictionary to STL map, G.Barrand, M.Maire 35 // 27.09.01 L.Urban , bug fixed (negative ener 36 // 27.09.01 L.Urban , bug fixed (negative energy deposit) 36 // 26.10.01 all static functions moved from .i 37 // 26.10.01 all static functions moved from .icc files (mma) 37 // 15.01.03 Add interfaces required for "cut p 38 // 15.01.03 Add interfaces required for "cut per region" (V.Ivanchenko) 38 // 12.03.03 Add warnings to obsolete interface 39 // 12.03.03 Add warnings to obsolete interfaces (V.Ivanchenko) 39 // 10.04.03 Add call to G4LossTableManager is 40 // 10.04.03 Add call to G4LossTableManager is particle is not registered (V.Ivanchenko) 40 // 41 // 41 // ------------------------------------------- 42 // ------------------------------------------------------------------- 42 43 43 #include "G4EnergyLossTables.hh" 44 #include "G4EnergyLossTables.hh" 44 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh" 45 #include "G4MaterialCutsCouple.hh" 46 #include "G4MaterialCutsCouple.hh" 46 #include "G4RegionStore.hh" 47 #include "G4RegionStore.hh" 47 #include "G4LossTableManager.hh" 48 #include "G4LossTableManager.hh" 48 49 49 50 50 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 51 << 52 /* 52 G4EnergyLossTablesHelper *G4EnergyLossTables:: << 53 G4ThreadLocal G4EnergyLossTablesHelper *G4EnergyLossTables::t = 0; 53 G4EnergyLossTablesHelper *G4EnergyLossTables:: << 54 G4ThreadLocal G4EnergyLossTablesHelper *G4EnergyLossTables::null_loss = 0; 54 G4ParticleDefinition* G4EnergyLossTables::last << 55 G4ThreadLocal G4ParticleDefinition* G4EnergyLossTables::lastParticle = 0; >> 56 G4ThreadLocal G4double G4EnergyLossTables::QQPositron = 1.0; // e_squared >> 57 G4ThreadLocal G4double G4EnergyLossTables::Chargesquare ; >> 58 G4ThreadLocal G4int G4EnergyLossTables::oldIndex = -1 ; >> 59 G4ThreadLocal G4double G4EnergyLossTables::rmin = 0. ; >> 60 G4ThreadLocal G4double G4EnergyLossTables::rmax = 0. ; >> 61 G4ThreadLocal G4double G4EnergyLossTables::Thigh = 0. ; >> 62 G4ThreadLocal G4int G4EnergyLossTables::let_counter = 0; >> 63 G4ThreadLocal G4int G4EnergyLossTables::let_max_num_warnings = 100; >> 64 G4ThreadLocal G4bool G4EnergyLossTables::first_loss = true; >> 65 >> 66 G4ThreadLocal G4EnergyLossTables::helper_map *G4EnergyLossTables::dict = 0; >> 67 */ >> 68 G4EnergyLossTablesHelper *G4EnergyLossTables::t = 0; >> 69 G4EnergyLossTablesHelper *G4EnergyLossTables::null_loss = 0; >> 70 G4ParticleDefinition* G4EnergyLossTables::lastParticle = 0; 55 G4double G4EnergyLossTables::QQPositron = 1.0; 71 G4double G4EnergyLossTables::QQPositron = 1.0; // e_squared 56 G4double G4EnergyLossTables::Chargesquare ; 72 G4double G4EnergyLossTables::Chargesquare ; 57 G4int G4EnergyLossTables::oldIndex = -1 ; 73 G4int G4EnergyLossTables::oldIndex = -1 ; 58 G4double G4EnergyLossTables::rmin = 0. ; 74 G4double G4EnergyLossTables::rmin = 0. ; 59 G4double G4EnergyLossTables::rmax = 0. ; 75 G4double G4EnergyLossTables::rmax = 0. ; 60 G4double G4EnergyLossTables::Thigh = 0. ; 76 G4double G4EnergyLossTables::Thigh = 0. ; 61 G4int G4EnergyLossTables::let_counter = 0; 77 G4int G4EnergyLossTables::let_counter = 0; 62 G4int G4EnergyLossTables::let_max_num_warni 78 G4int G4EnergyLossTables::let_max_num_warnings = 100; 63 G4bool G4EnergyLossTables::first_loss = true 79 G4bool G4EnergyLossTables::first_loss = true; 64 80 65 G4EnergyLossTables::helper_map *G4EnergyLossTa << 81 G4EnergyLossTables::helper_map *G4EnergyLossTables::dict = 0; 66 82 67 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 68 84 69 G4EnergyLossTablesHelper::G4EnergyLossTablesHe 85 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper( 70 const G4PhysicsTable* aDEDXTable, 86 const G4PhysicsTable* aDEDXTable, 71 const G4PhysicsTable* aRangeTable, 87 const G4PhysicsTable* aRangeTable, 72 const G4PhysicsTable* anInverseRangeTable, 88 const G4PhysicsTable* anInverseRangeTable, 73 const G4PhysicsTable* aLabTimeTable, 89 const G4PhysicsTable* aLabTimeTable, 74 const G4PhysicsTable* aProperTimeTable, 90 const G4PhysicsTable* aProperTimeTable, 75 G4double aLowestKineticEnergy, 91 G4double aLowestKineticEnergy, 76 G4double aHighestKineticEnergy, 92 G4double aHighestKineticEnergy, 77 G4double aMassRatio, 93 G4double aMassRatio, 78 G4int aNumberOfBins) 94 G4int aNumberOfBins) 79 : 95 : 80 theDEDXTable(aDEDXTable), theRangeTable(aRan 96 theDEDXTable(aDEDXTable), theRangeTable(aRangeTable), 81 theInverseRangeTable(anInverseRangeTable), 97 theInverseRangeTable(anInverseRangeTable), 82 theLabTimeTable(aLabTimeTable), 98 theLabTimeTable(aLabTimeTable), 83 theProperTimeTable(aProperTimeTable), 99 theProperTimeTable(aProperTimeTable), 84 theLowestKineticEnergy(aLowestKineticEnergy) 100 theLowestKineticEnergy(aLowestKineticEnergy), 85 theHighestKineticEnergy(aHighestKineticEnerg 101 theHighestKineticEnergy(aHighestKineticEnergy), 86 theMassRatio(aMassRatio), 102 theMassRatio(aMassRatio), 87 theNumberOfBins(aNumberOfBins) 103 theNumberOfBins(aNumberOfBins) 88 { } 104 { } 89 105 90 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 91 107 92 G4EnergyLossTablesHelper::G4EnergyLossTablesHe 108 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper() 93 { 109 { 94 theLowestKineticEnergy = 0.0; 110 theLowestKineticEnergy = 0.0; 95 theHighestKineticEnergy= 0.0; 111 theHighestKineticEnergy= 0.0; 96 theMassRatio = 0.0; 112 theMassRatio = 0.0; 97 theNumberOfBins = 0; 113 theNumberOfBins = 0; 98 theDEDXTable = theRangeTable = theInverseRan << 99 = theProperTimeTable = nullptr; << 100 } 114 } 101 115 102 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 103 117 104 void G4EnergyLossTables::Register( 118 void G4EnergyLossTables::Register( 105 const G4ParticleDefinition* p, 119 const G4ParticleDefinition* p, 106 const G4PhysicsTable* tDEDX, 120 const G4PhysicsTable* tDEDX, 107 const G4PhysicsTable* tRange, 121 const G4PhysicsTable* tRange, 108 const G4PhysicsTable* tInverseRange, 122 const G4PhysicsTable* tInverseRange, 109 const G4PhysicsTable* tLabTime, 123 const G4PhysicsTable* tLabTime, 110 const G4PhysicsTable* tProperTime, 124 const G4PhysicsTable* tProperTime, 111 G4double lowestKineticEnergy, 125 G4double lowestKineticEnergy, 112 G4double highestKineticEnergy, 126 G4double highestKineticEnergy, 113 G4double massRatio, 127 G4double massRatio, 114 G4int NumberOfBins) 128 G4int NumberOfBins) 115 { 129 { 116 if (!dict) dict = new G4EnergyLossTables::he 130 if (!dict) dict = new G4EnergyLossTables::helper_map; 117 if (!null_loss) null_loss = new G4EnergyLoss 131 if (!null_loss) null_loss = new G4EnergyLossTablesHelper; 118 if (!t) t = new G4EnergyLossTablesHelper; 132 if (!t) t = new G4EnergyLossTablesHelper; 119 133 120 (*dict)[p]= G4EnergyLossTablesHelper(tDEDX, 134 (*dict)[p]= G4EnergyLossTablesHelper(tDEDX, tRange,tInverseRange, 121 tLabTime,tProperTime,lowes 135 tLabTime,tProperTime,lowestKineticEnergy, 122 highestKineticEnergy, massRatio,Number 136 highestKineticEnergy, massRatio,NumberOfBins); 123 137 124 *t = GetTables(p) ; // important for cach 138 *t = GetTables(p) ; // important for cache !!!!! 125 lastParticle = (G4ParticleDefinition*) p ; 139 lastParticle = (G4ParticleDefinition*) p ; 126 Chargesquare = (p->GetPDGCharge())*(p->GetPD 140 Chargesquare = (p->GetPDGCharge())*(p->GetPDGCharge())/ 127 QQPositron ; 141 QQPositron ; 128 if (first_loss ) { 142 if (first_loss ) { 129 *null_loss = G4EnergyLossTablesHelper( << 143 *null_loss = G4EnergyLossTablesHelper(0, 0, 0, 0, 0, 0.0, 0.0, 0.0, 0); 130 nullptr, nullptr, nullptr, nu << 131 first_loss = false; 144 first_loss = false; 132 } 145 } 133 } 146 } 134 147 135 //....oooOO0OOooo........oooOO0OOooo........oo 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 136 149 137 const G4PhysicsTable* G4EnergyLossTables::GetD 150 const G4PhysicsTable* G4EnergyLossTables::GetDEDXTable( 138 const G4ParticleDefinition* p) 151 const G4ParticleDefinition* p) 139 { 152 { 140 if (!dict) dict = new G4EnergyLossTables::he 153 if (!dict) dict = new G4EnergyLossTables::helper_map; 141 helper_map::iterator it; 154 helper_map::iterator it; 142 if((it=dict->find(p))==dict->end()) return n << 155 if((it=dict->find(p))==dict->end()) return 0; 143 return (*it).second.theDEDXTable; 156 return (*it).second.theDEDXTable; 144 } 157 } 145 158 146 //....oooOO0OOooo........oooOO0OOooo........oo 159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 147 160 148 const G4PhysicsTable* G4EnergyLossTables::GetR 161 const G4PhysicsTable* G4EnergyLossTables::GetRangeTable( 149 const G4ParticleDefinition* p) 162 const G4ParticleDefinition* p) 150 { 163 { 151 if (!dict) dict = new G4EnergyLossTables::he 164 if (!dict) dict = new G4EnergyLossTables::helper_map; 152 helper_map::iterator it; 165 helper_map::iterator it; 153 if((it=dict->find(p))==dict->end()) return n << 166 if((it=dict->find(p))==dict->end()) return 0; 154 return (*it).second.theRangeTable; 167 return (*it).second.theRangeTable; 155 } 168 } 156 169 157 //....oooOO0OOooo........oooOO0OOooo........oo 170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 158 171 159 const G4PhysicsTable* G4EnergyLossTables::GetI 172 const G4PhysicsTable* G4EnergyLossTables::GetInverseRangeTable( 160 const G4ParticleDefinition* p) 173 const G4ParticleDefinition* p) 161 { 174 { 162 if (!dict) dict = new G4EnergyLossTables::he 175 if (!dict) dict = new G4EnergyLossTables::helper_map; 163 helper_map::iterator it; 176 helper_map::iterator it; 164 if((it=dict->find(p))==dict->end()) return n << 177 if((it=dict->find(p))==dict->end()) return 0; 165 return (*it).second.theInverseRangeTable; 178 return (*it).second.theInverseRangeTable; 166 } 179 } 167 180 168 //....oooOO0OOooo........oooOO0OOooo........oo 181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 169 182 170 const G4PhysicsTable* G4EnergyLossTables::GetL 183 const G4PhysicsTable* G4EnergyLossTables::GetLabTimeTable( 171 const G4ParticleDefinition* p) 184 const G4ParticleDefinition* p) 172 { 185 { 173 if (!dict) dict = new G4EnergyLossTables::he 186 if (!dict) dict = new G4EnergyLossTables::helper_map; 174 helper_map::iterator it; 187 helper_map::iterator it; 175 if((it=dict->find(p))==dict->end()) return n << 188 if((it=dict->find(p))==dict->end()) return 0; 176 return (*it).second.theLabTimeTable; 189 return (*it).second.theLabTimeTable; 177 } 190 } 178 191 179 //....oooOO0OOooo........oooOO0OOooo........oo 192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 180 193 181 const G4PhysicsTable* G4EnergyLossTables::GetP 194 const G4PhysicsTable* G4EnergyLossTables::GetProperTimeTable( 182 const G4ParticleDefinition* p) 195 const G4ParticleDefinition* p) 183 { 196 { 184 if (!dict) dict = new G4EnergyLossTables::he 197 if (!dict) dict = new G4EnergyLossTables::helper_map; 185 helper_map::iterator it; 198 helper_map::iterator it; 186 if((it=dict->find(p))==dict->end()) return n << 199 if((it=dict->find(p))==dict->end()) return 0; 187 return (*it).second.theProperTimeTable; 200 return (*it).second.theProperTimeTable; 188 } 201 } 189 202 190 //....oooOO0OOooo........oooOO0OOooo........oo 203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 191 204 192 G4EnergyLossTablesHelper G4EnergyLossTables::G 205 G4EnergyLossTablesHelper G4EnergyLossTables::GetTables( 193 const G4ParticleDefinition* p) 206 const G4ParticleDefinition* p) 194 { 207 { 195 if (!dict) dict = new G4EnergyLossTables::he 208 if (!dict) dict = new G4EnergyLossTables::helper_map; 196 if (!null_loss) null_loss = new G4EnergyLoss 209 if (!null_loss) null_loss = new G4EnergyLossTablesHelper; 197 210 198 helper_map::iterator it; 211 helper_map::iterator it; 199 if ((it=dict->find(p))==dict->end()) { 212 if ((it=dict->find(p))==dict->end()) { 200 return *null_loss; 213 return *null_loss; 201 } 214 } 202 return (*it).second; 215 return (*it).second; 203 } 216 } 204 217 205 //....oooOO0OOooo........oooOO0OOooo........oo 218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 206 219 207 G4double G4EnergyLossTables::GetDEDX( 220 G4double G4EnergyLossTables::GetDEDX( 208 const G4ParticleDefinition *aParticle, 221 const G4ParticleDefinition *aParticle, 209 G4double KineticEnergy, 222 G4double KineticEnergy, 210 const G4Material *aMaterial) 223 const G4Material *aMaterial) 211 { 224 { 212 if (!t) t = new G4EnergyLossTablesHelper; 225 if (!t) t = new G4EnergyLossTablesHelper; 213 226 214 CPRWarning(); 227 CPRWarning(); 215 if(aParticle != (const G4ParticleDefinition* 228 if(aParticle != (const G4ParticleDefinition*) lastParticle) 216 { 229 { 217 *t= GetTables(aParticle); 230 *t= GetTables(aParticle); 218 lastParticle = (G4ParticleDefinition*) aPa 231 lastParticle = (G4ParticleDefinition*) aParticle ; 219 Chargesquare = (aParticle->GetPDGCharge()) 232 Chargesquare = (aParticle->GetPDGCharge())* 220 (aParticle->GetPDGCharge()) 233 (aParticle->GetPDGCharge())/ 221 QQPositron ; 234 QQPositron ; 222 oldIndex = -1 ; 235 oldIndex = -1 ; 223 } 236 } 224 const G4PhysicsTable* dEdxTable= t->theDEDX 237 const G4PhysicsTable* dEdxTable= t->theDEDXTable; 225 if (!dEdxTable) { 238 if (!dEdxTable) { 226 ParticleHaveNoLoss(aParticle,"dEdx"); 239 ParticleHaveNoLoss(aParticle,"dEdx"); 227 return 0.0; 240 return 0.0; 228 } 241 } 229 242 230 G4int materialIndex = (G4int)aMaterial->GetI << 243 G4int materialIndex = aMaterial->GetIndex(); 231 G4double scaledKineticEnergy = KineticEnergy 244 G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio; 232 G4double dEdx; 245 G4double dEdx; 233 G4bool isOut; 246 G4bool isOut; 234 247 235 if (scaledKineticEnergy<t->theLowestKineticE 248 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 236 249 237 dEdx =(*dEdxTable)(materialIndex)->GetVal 250 dEdx =(*dEdxTable)(materialIndex)->GetValue( 238 t->theLowestKineticEnergy,isOut) 251 t->theLowestKineticEnergy,isOut) 239 *std::sqrt(scaledKineticEnergy/t->t 252 *std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy); 240 253 241 } else if (scaledKineticEnergy>t->theHighest 254 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 242 255 243 dEdx = (*dEdxTable)(materialIndex)->GetVa 256 dEdx = (*dEdxTable)(materialIndex)->GetValue( 244 t->theHighestKineticEnergy,isOut); 257 t->theHighestKineticEnergy,isOut); 245 258 246 } else { 259 } else { 247 260 248 dEdx = (*dEdxTable)(materialIndex)->GetVal 261 dEdx = (*dEdxTable)(materialIndex)->GetValue( 249 scaledKineticEnergy,isOut); 262 scaledKineticEnergy,isOut); 250 263 251 } 264 } 252 265 253 return dEdx*Chargesquare; 266 return dEdx*Chargesquare; 254 } 267 } 255 268 256 //....oooOO0OOooo........oooOO0OOooo........oo 269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 257 270 258 G4double G4EnergyLossTables::GetLabTime( 271 G4double G4EnergyLossTables::GetLabTime( 259 const G4ParticleDefinition *aParticle, 272 const G4ParticleDefinition *aParticle, 260 G4double KineticEnergy, 273 G4double KineticEnergy, 261 const G4Material *aMaterial) 274 const G4Material *aMaterial) 262 { 275 { 263 if (!t) t = new G4EnergyLossTablesHelper; 276 if (!t) t = new G4EnergyLossTablesHelper; 264 277 265 CPRWarning(); 278 CPRWarning(); 266 if(aParticle != (const G4ParticleDefinition* 279 if(aParticle != (const G4ParticleDefinition*) lastParticle) 267 { 280 { 268 *t= GetTables(aParticle); 281 *t= GetTables(aParticle); 269 lastParticle = (G4ParticleDefinition*) aPa 282 lastParticle = (G4ParticleDefinition*) aParticle ; 270 oldIndex = -1 ; 283 oldIndex = -1 ; 271 } 284 } 272 const G4PhysicsTable* labtimeTable= t->theLa 285 const G4PhysicsTable* labtimeTable= t->theLabTimeTable; 273 if (!labtimeTable) { 286 if (!labtimeTable) { 274 ParticleHaveNoLoss(aParticle,"LabTime"); 287 ParticleHaveNoLoss(aParticle,"LabTime"); 275 return 0.0; 288 return 0.0; 276 } 289 } 277 290 278 const G4double parlowen=0.4 , ppar=0.5-parlo 291 const G4double parlowen=0.4 , ppar=0.5-parlowen ; 279 G4int materialIndex = (G4int)aMaterial->GetI << 292 G4int materialIndex = aMaterial->GetIndex(); 280 G4double scaledKineticEnergy = KineticEnergy 293 G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio; 281 G4double time; 294 G4double time; 282 G4bool isOut; 295 G4bool isOut; 283 296 284 if (scaledKineticEnergy<t->theLowestKineticE 297 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 285 298 286 time = std::exp(ppar*std::log(scaledKinet 299 time = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))* 287 (*labtimeTable)(materialIndex)->Ge 300 (*labtimeTable)(materialIndex)->GetValue( 288 t->theLowestKineticEnergy,isOut) 301 t->theLowestKineticEnergy,isOut); 289 302 290 303 291 } else if (scaledKineticEnergy>t->theHighest 304 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 292 305 293 time = (*labtimeTable)(materialIndex)->Ge 306 time = (*labtimeTable)(materialIndex)->GetValue( 294 t->theHighestKineticEnergy,isOut 307 t->theHighestKineticEnergy,isOut); 295 308 296 } else { 309 } else { 297 310 298 time = (*labtimeTable)(materialIndex)->Get 311 time = (*labtimeTable)(materialIndex)->GetValue( 299 scaledKineticEnergy,isOut); 312 scaledKineticEnergy,isOut); 300 313 301 } 314 } 302 315 303 return time/t->theMassRatio ; 316 return time/t->theMassRatio ; 304 } 317 } 305 318 306 //....oooOO0OOooo........oooOO0OOooo........oo 319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 307 320 308 G4double G4EnergyLossTables::GetDeltaLabTime( 321 G4double G4EnergyLossTables::GetDeltaLabTime( 309 const G4ParticleDefinition *aParticle, 322 const G4ParticleDefinition *aParticle, 310 G4double KineticEnergyStart, 323 G4double KineticEnergyStart, 311 G4double KineticEnergyEnd, 324 G4double KineticEnergyEnd, 312 const G4Material *aMaterial) 325 const G4Material *aMaterial) 313 { 326 { 314 if (!t) t = new G4EnergyLossTablesHelper; 327 if (!t) t = new G4EnergyLossTablesHelper; 315 328 316 CPRWarning(); 329 CPRWarning(); 317 if(aParticle != (const G4ParticleDefinition* 330 if(aParticle != (const G4ParticleDefinition*) lastParticle) 318 { 331 { 319 *t= GetTables(aParticle); 332 *t= GetTables(aParticle); 320 lastParticle = (G4ParticleDefinition*) aPa 333 lastParticle = (G4ParticleDefinition*) aParticle ; 321 oldIndex = -1 ; 334 oldIndex = -1 ; 322 } 335 } 323 const G4PhysicsTable* labtimeTable= t->theLa 336 const G4PhysicsTable* labtimeTable= t->theLabTimeTable; 324 if (!labtimeTable) { 337 if (!labtimeTable) { 325 ParticleHaveNoLoss(aParticle,"LabTime"); 338 ParticleHaveNoLoss(aParticle,"LabTime"); 326 return 0.0; 339 return 0.0; 327 } 340 } 328 341 329 const G4double parlowen=0.4 , ppar=0.5-parlo 342 const G4double parlowen=0.4 , ppar=0.5-parlowen ; 330 const G4double dToverT = 0.05 , facT = 1. -d 343 const G4double dToverT = 0.05 , facT = 1. -dToverT ; 331 G4double timestart,timeend,deltatime,dTT; 344 G4double timestart,timeend,deltatime,dTT; 332 G4bool isOut; 345 G4bool isOut; 333 346 334 G4int materialIndex = (G4int)aMaterial->GetI << 347 G4int materialIndex = aMaterial->GetIndex(); 335 G4double scaledKineticEnergy = KineticEnergy 348 G4double scaledKineticEnergy = KineticEnergyStart*t->theMassRatio; 336 349 337 if (scaledKineticEnergy<t->theLowestKineticE 350 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 338 351 339 timestart = std::exp(ppar*std::log(scaled 352 timestart = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))* 340 (*labtimeTable)(materialIndex) 353 (*labtimeTable)(materialIndex)->GetValue( 341 t->theLowestKineticEnergy,isOu 354 t->theLowestKineticEnergy,isOut); 342 355 343 356 344 } else if (scaledKineticEnergy>t->theHighest 357 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 345 358 346 timestart = (*labtimeTable)(materialIndex 359 timestart = (*labtimeTable)(materialIndex)->GetValue( 347 t->theHighestKineticEnergy,isO 360 t->theHighestKineticEnergy,isOut); 348 361 349 } else { 362 } else { 350 363 351 timestart = (*labtimeTable)(materialIndex) 364 timestart = (*labtimeTable)(materialIndex)->GetValue( 352 scaledKineticEnergy,isOut); 365 scaledKineticEnergy,isOut); 353 366 354 } 367 } 355 368 356 dTT = (KineticEnergyStart - KineticEnergyEnd 369 dTT = (KineticEnergyStart - KineticEnergyEnd)/KineticEnergyStart ; 357 370 358 if( dTT < dToverT ) 371 if( dTT < dToverT ) 359 scaledKineticEnergy = facT*KineticEnergySt 372 scaledKineticEnergy = facT*KineticEnergyStart*t->theMassRatio; 360 else 373 else 361 scaledKineticEnergy = KineticEnergyEnd*t-> 374 scaledKineticEnergy = KineticEnergyEnd*t->theMassRatio; 362 375 363 if (scaledKineticEnergy<t->theLowestKineticE 376 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 364 377 365 timeend = std::exp(ppar*std::log(scaledKi 378 timeend = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))* 366 (*labtimeTable)(materialIndex) 379 (*labtimeTable)(materialIndex)->GetValue( 367 t->theLowestKineticEnergy,isOu 380 t->theLowestKineticEnergy,isOut); 368 381 369 382 370 } else if (scaledKineticEnergy>t->theHighest 383 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 371 384 372 timeend = (*labtimeTable)(materialIndex)- 385 timeend = (*labtimeTable)(materialIndex)->GetValue( 373 t->theHighestKineticEnergy,isO 386 t->theHighestKineticEnergy,isOut); 374 387 375 } else { 388 } else { 376 389 377 timeend = (*labtimeTable)(materialIndex)-> 390 timeend = (*labtimeTable)(materialIndex)->GetValue( 378 scaledKineticEnergy,isOut); 391 scaledKineticEnergy,isOut); 379 392 380 } 393 } 381 394 382 deltatime = timestart - timeend ; 395 deltatime = timestart - timeend ; 383 396 384 if( dTT < dToverT ) 397 if( dTT < dToverT ) 385 deltatime *= dTT/dToverT; 398 deltatime *= dTT/dToverT; 386 399 387 return deltatime/t->theMassRatio ; 400 return deltatime/t->theMassRatio ; 388 } 401 } 389 402 390 //....oooOO0OOooo........oooOO0OOooo........oo 403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 391 404 392 G4double G4EnergyLossTables::GetProperTime( 405 G4double G4EnergyLossTables::GetProperTime( 393 const G4ParticleDefinition *aParticle, 406 const G4ParticleDefinition *aParticle, 394 G4double KineticEnergy, 407 G4double KineticEnergy, 395 const G4Material *aMaterial) 408 const G4Material *aMaterial) 396 { 409 { 397 if (!t) t = new G4EnergyLossTablesHelper; 410 if (!t) t = new G4EnergyLossTablesHelper; 398 411 399 CPRWarning(); 412 CPRWarning(); 400 if(aParticle != (const G4ParticleDefinition* 413 if(aParticle != (const G4ParticleDefinition*) lastParticle) 401 { 414 { 402 *t= GetTables(aParticle); 415 *t= GetTables(aParticle); 403 lastParticle = (G4ParticleDefinition*) aPa 416 lastParticle = (G4ParticleDefinition*) aParticle ; 404 oldIndex = -1 ; 417 oldIndex = -1 ; 405 } 418 } 406 const G4PhysicsTable* propertimeTable= t->th 419 const G4PhysicsTable* propertimeTable= t->theProperTimeTable; 407 if (!propertimeTable) { 420 if (!propertimeTable) { 408 ParticleHaveNoLoss(aParticle,"ProperTime") 421 ParticleHaveNoLoss(aParticle,"ProperTime"); 409 return 0.0; 422 return 0.0; 410 } 423 } 411 424 412 const G4double parlowen=0.4 , ppar=0.5-parlo 425 const G4double parlowen=0.4 , ppar=0.5-parlowen ; 413 G4int materialIndex = (G4int)aMaterial->GetI << 426 G4int materialIndex = aMaterial->GetIndex(); 414 G4double scaledKineticEnergy = KineticEnergy 427 G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio; 415 G4double time; 428 G4double time; 416 G4bool isOut; 429 G4bool isOut; 417 430 418 if (scaledKineticEnergy<t->theLowestKineticE 431 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 419 432 420 time = std::exp(ppar*std::log(scaledKinet 433 time = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))* 421 (*propertimeTable)(materialIndex)- 434 (*propertimeTable)(materialIndex)->GetValue( 422 t->theLowestKineticEnergy,isOut) 435 t->theLowestKineticEnergy,isOut); 423 436 424 437 425 } else if (scaledKineticEnergy>t->theHighest 438 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 426 439 427 time = (*propertimeTable)(materialIndex)- 440 time = (*propertimeTable)(materialIndex)->GetValue( 428 t->theHighestKineticEnergy,isOut 441 t->theHighestKineticEnergy,isOut); 429 442 430 } else { 443 } else { 431 444 432 time = (*propertimeTable)(materialIndex)-> 445 time = (*propertimeTable)(materialIndex)->GetValue( 433 scaledKineticEnergy,isOut); 446 scaledKineticEnergy,isOut); 434 447 435 } 448 } 436 449 437 return time/t->theMassRatio ; 450 return time/t->theMassRatio ; 438 } 451 } 439 452 440 //....oooOO0OOooo........oooOO0OOooo........oo 453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 441 454 442 G4double G4EnergyLossTables::GetDeltaProperTim 455 G4double G4EnergyLossTables::GetDeltaProperTime( 443 const G4ParticleDefinition *aParticle, 456 const G4ParticleDefinition *aParticle, 444 G4double KineticEnergyStart, 457 G4double KineticEnergyStart, 445 G4double KineticEnergyEnd, 458 G4double KineticEnergyEnd, 446 const G4Material *aMaterial) 459 const G4Material *aMaterial) 447 { 460 { 448 if (!t) t = new G4EnergyLossTablesHelper; 461 if (!t) t = new G4EnergyLossTablesHelper; 449 462 450 CPRWarning(); 463 CPRWarning(); 451 if(aParticle != (const G4ParticleDefinition* 464 if(aParticle != (const G4ParticleDefinition*) lastParticle) 452 { 465 { 453 *t= GetTables(aParticle); 466 *t= GetTables(aParticle); 454 lastParticle = (G4ParticleDefinition*) aPa 467 lastParticle = (G4ParticleDefinition*) aParticle ; 455 oldIndex = -1 ; 468 oldIndex = -1 ; 456 } 469 } 457 const G4PhysicsTable* propertimeTable= t->th 470 const G4PhysicsTable* propertimeTable= t->theProperTimeTable; 458 if (!propertimeTable) { 471 if (!propertimeTable) { 459 ParticleHaveNoLoss(aParticle,"ProperTime") 472 ParticleHaveNoLoss(aParticle,"ProperTime"); 460 return 0.0; 473 return 0.0; 461 } 474 } 462 475 463 const G4double parlowen=0.4 , ppar=0.5-parlo 476 const G4double parlowen=0.4 , ppar=0.5-parlowen ; 464 const G4double dToverT = 0.05 , facT = 1. -d 477 const G4double dToverT = 0.05 , facT = 1. -dToverT ; 465 G4double timestart,timeend,deltatime,dTT; 478 G4double timestart,timeend,deltatime,dTT; 466 G4bool isOut; 479 G4bool isOut; 467 480 468 G4int materialIndex = (G4int)aMaterial->GetI << 481 G4int materialIndex = aMaterial->GetIndex(); 469 G4double scaledKineticEnergy = KineticEnergy 482 G4double scaledKineticEnergy = KineticEnergyStart*t->theMassRatio; 470 483 471 if (scaledKineticEnergy<t->theLowestKineticE 484 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 472 485 473 timestart = std::exp(ppar*std::log(scaled 486 timestart = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))* 474 (*propertimeTable)(materialInd 487 (*propertimeTable)(materialIndex)->GetValue( 475 t->theLowestKineticEnergy,isOu 488 t->theLowestKineticEnergy,isOut); 476 489 477 490 478 } else if (scaledKineticEnergy>t->theHighest 491 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 479 492 480 timestart = (*propertimeTable)(materialIn 493 timestart = (*propertimeTable)(materialIndex)->GetValue( 481 t->theHighestKineticEnergy,isO 494 t->theHighestKineticEnergy,isOut); 482 495 483 } else { 496 } else { 484 497 485 timestart = (*propertimeTable)(materialInd 498 timestart = (*propertimeTable)(materialIndex)->GetValue( 486 scaledKineticEnergy,isOut); 499 scaledKineticEnergy,isOut); 487 500 488 } 501 } 489 502 490 dTT = (KineticEnergyStart - KineticEnergyEnd 503 dTT = (KineticEnergyStart - KineticEnergyEnd)/KineticEnergyStart ; 491 504 492 if( dTT < dToverT ) 505 if( dTT < dToverT ) 493 scaledKineticEnergy = facT*KineticEnergySt 506 scaledKineticEnergy = facT*KineticEnergyStart*t->theMassRatio; 494 else 507 else 495 scaledKineticEnergy = KineticEnergyEnd*t-> 508 scaledKineticEnergy = KineticEnergyEnd*t->theMassRatio; 496 509 497 if (scaledKineticEnergy<t->theLowestKineticE 510 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 498 511 499 timeend = std::exp(ppar*std::log(scaledKi 512 timeend = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))* 500 (*propertimeTable)(materialInd 513 (*propertimeTable)(materialIndex)->GetValue( 501 t->theLowestKineticEnergy,isOu 514 t->theLowestKineticEnergy,isOut); 502 515 503 516 504 } else if (scaledKineticEnergy>t->theHighest 517 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 505 518 506 timeend = (*propertimeTable)(materialInde 519 timeend = (*propertimeTable)(materialIndex)->GetValue( 507 t->theHighestKineticEnergy,isO 520 t->theHighestKineticEnergy,isOut); 508 521 509 } else { 522 } else { 510 523 511 timeend = (*propertimeTable)(materialIndex 524 timeend = (*propertimeTable)(materialIndex)->GetValue( 512 scaledKineticEnergy,isOut); 525 scaledKineticEnergy,isOut); 513 526 514 } 527 } 515 528 516 deltatime = timestart - timeend ; 529 deltatime = timestart - timeend ; 517 530 518 if( dTT < dToverT ) 531 if( dTT < dToverT ) 519 deltatime *= dTT/dToverT ; 532 deltatime *= dTT/dToverT ; 520 533 521 return deltatime/t->theMassRatio ; 534 return deltatime/t->theMassRatio ; 522 } 535 } 523 536 524 //....oooOO0OOooo........oooOO0OOooo........oo 537 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 525 538 526 G4double G4EnergyLossTables::GetRange( 539 G4double G4EnergyLossTables::GetRange( 527 const G4ParticleDefinition *aParticle, 540 const G4ParticleDefinition *aParticle, 528 G4double KineticEnergy, 541 G4double KineticEnergy, 529 const G4Material *aMaterial) 542 const G4Material *aMaterial) 530 { 543 { 531 if (!t) t = new G4EnergyLossTablesHelper; 544 if (!t) t = new G4EnergyLossTablesHelper; 532 545 533 CPRWarning(); 546 CPRWarning(); 534 if(aParticle != (const G4ParticleDefinition* 547 if(aParticle != (const G4ParticleDefinition*) lastParticle) 535 { 548 { 536 *t= GetTables(aParticle); 549 *t= GetTables(aParticle); 537 lastParticle = (G4ParticleDefinition*) aPa 550 lastParticle = (G4ParticleDefinition*) aParticle ; 538 Chargesquare = (aParticle->GetPDGCharge()) 551 Chargesquare = (aParticle->GetPDGCharge())* 539 (aParticle->GetPDGCharge()) 552 (aParticle->GetPDGCharge())/ 540 QQPositron ; 553 QQPositron ; 541 oldIndex = -1 ; 554 oldIndex = -1 ; 542 } 555 } 543 const G4PhysicsTable* rangeTable= t->theRang 556 const G4PhysicsTable* rangeTable= t->theRangeTable; 544 const G4PhysicsTable* dEdxTable= t->theDEDX 557 const G4PhysicsTable* dEdxTable= t->theDEDXTable; 545 if (!rangeTable) { 558 if (!rangeTable) { 546 ParticleHaveNoLoss(aParticle,"Range"); 559 ParticleHaveNoLoss(aParticle,"Range"); 547 return 0.0; 560 return 0.0; 548 } 561 } 549 562 550 G4int materialIndex = (G4int)aMaterial->GetI << 563 G4int materialIndex = aMaterial->GetIndex(); 551 G4double scaledKineticEnergy = KineticEnergy 564 G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio; 552 G4double Range; 565 G4double Range; 553 G4bool isOut; 566 G4bool isOut; 554 567 555 if (scaledKineticEnergy<t->theLowestKineticE 568 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 556 569 557 Range = std::sqrt(scaledKineticEnergy/t->t 570 Range = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy)* 558 (*rangeTable)(materialIndex)->GetV 571 (*rangeTable)(materialIndex)->GetValue( 559 t->theLowestKineticEnergy,isOut) 572 t->theLowestKineticEnergy,isOut); 560 573 561 } else if (scaledKineticEnergy>t->theHighest 574 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 562 575 563 Range = (*rangeTable)(materialIndex)->GetV 576 Range = (*rangeTable)(materialIndex)->GetValue( 564 t->theHighestKineticEnergy,isOut)+ 577 t->theHighestKineticEnergy,isOut)+ 565 (scaledKineticEnergy-t->theHighest 578 (scaledKineticEnergy-t->theHighestKineticEnergy)/ 566 (*dEdxTable)(materialIndex)->GetVa 579 (*dEdxTable)(materialIndex)->GetValue( 567 t->theHighestKineticEnergy,isOut 580 t->theHighestKineticEnergy,isOut); 568 581 569 } else { 582 } else { 570 583 571 Range = (*rangeTable)(materialIndex)->GetV 584 Range = (*rangeTable)(materialIndex)->GetValue( 572 scaledKineticEnergy,isOut); 585 scaledKineticEnergy,isOut); 573 586 574 } 587 } 575 588 576 return Range/(Chargesquare*t->theMassRatio); 589 return Range/(Chargesquare*t->theMassRatio); 577 } 590 } 578 591 579 //....oooOO0OOooo........oooOO0OOooo........oo 592 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 580 593 581 G4double G4EnergyLossTables::GetPreciseEnergyF 594 G4double G4EnergyLossTables::GetPreciseEnergyFromRange( 582 const G4P 595 const G4ParticleDefinition *aParticle, 583 G4d 596 G4double range, 584 const G4M 597 const G4Material *aMaterial) 585 // it returns the value of the kinetic energy 598 // it returns the value of the kinetic energy for a given range 586 { 599 { 587 if (!t) t = new G4EnergyLossTablesHelper; 600 if (!t) t = new G4EnergyLossTablesHelper; 588 601 589 CPRWarning(); 602 CPRWarning(); 590 if( aParticle != (const G4ParticleDefinition 603 if( aParticle != (const G4ParticleDefinition*) lastParticle) 591 { 604 { 592 *t= GetTables(aParticle); 605 *t= GetTables(aParticle); 593 lastParticle = (G4ParticleDefinition*) aPa 606 lastParticle = (G4ParticleDefinition*) aParticle; 594 Chargesquare = (aParticle->GetPDGCharge()) 607 Chargesquare = (aParticle->GetPDGCharge())* 595 (aParticle->GetPDGCharge()) 608 (aParticle->GetPDGCharge())/ 596 QQPositron ; 609 QQPositron ; 597 oldIndex = -1 ; 610 oldIndex = -1 ; 598 } 611 } 599 const G4PhysicsTable* dEdxTable= t->theDEDX 612 const G4PhysicsTable* dEdxTable= t->theDEDXTable; 600 const G4PhysicsTable* inverseRangeTable= t- 613 const G4PhysicsTable* inverseRangeTable= t->theInverseRangeTable; 601 if (!inverseRangeTable) { 614 if (!inverseRangeTable) { 602 ParticleHaveNoLoss(aParticle,"InverseRange 615 ParticleHaveNoLoss(aParticle,"InverseRange"); 603 return 0.0; 616 return 0.0; 604 } 617 } 605 618 606 G4double scaledrange,scaledKineticEnergy ; 619 G4double scaledrange,scaledKineticEnergy ; 607 G4bool isOut ; 620 G4bool isOut ; 608 621 609 G4int materialIndex = (G4int)aMaterial->GetI << 622 G4int materialIndex = aMaterial->GetIndex() ; 610 623 611 if(materialIndex != oldIndex) 624 if(materialIndex != oldIndex) 612 { 625 { 613 oldIndex = materialIndex ; 626 oldIndex = materialIndex ; 614 rmin = (*inverseRangeTable)(materialIndex) 627 rmin = (*inverseRangeTable)(materialIndex)-> 615 GetLowEdgeEnergy 628 GetLowEdgeEnergy(0) ; 616 rmax = (*inverseRangeTable)(materialIndex) 629 rmax = (*inverseRangeTable)(materialIndex)-> 617 GetLowEdgeEnergy(t->theNumb 630 GetLowEdgeEnergy(t->theNumberOfBins-2) ; 618 Thigh = (*inverseRangeTable)(materialIndex 631 Thigh = (*inverseRangeTable)(materialIndex)-> 619 GetValue(rmax,is 632 GetValue(rmax,isOut) ; 620 } 633 } 621 634 622 scaledrange = range*Chargesquare*t->theMassR 635 scaledrange = range*Chargesquare*t->theMassRatio ; 623 636 624 if(scaledrange < rmin) 637 if(scaledrange < rmin) 625 { 638 { 626 scaledKineticEnergy = t->theLowestKineticE 639 scaledKineticEnergy = t->theLowestKineticEnergy* 627 scaledrange*scaledrange/(rm 640 scaledrange*scaledrange/(rmin*rmin) ; 628 } 641 } 629 else 642 else 630 { 643 { 631 if(scaledrange < rmax) 644 if(scaledrange < rmax) 632 { 645 { 633 scaledKineticEnergy = (*inverseRangeTabl 646 scaledKineticEnergy = (*inverseRangeTable)(materialIndex)-> 634 GetValue( scaled 647 GetValue( scaledrange,isOut) ; 635 } 648 } 636 else 649 else 637 { 650 { 638 scaledKineticEnergy = Thigh + 651 scaledKineticEnergy = Thigh + 639 (scaledrange-rmax)* 652 (scaledrange-rmax)* 640 (*dEdxTable)(materialInd 653 (*dEdxTable)(materialIndex)-> 641 GetValue(Thig 654 GetValue(Thigh,isOut) ; 642 } 655 } 643 } 656 } 644 657 645 return scaledKineticEnergy/t->theMassRatio ; 658 return scaledKineticEnergy/t->theMassRatio ; 646 } 659 } 647 660 648 //....oooOO0OOooo........oooOO0OOooo........oo 661 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 649 662 650 G4double G4EnergyLossTables::GetPreciseDEDX( 663 G4double G4EnergyLossTables::GetPreciseDEDX( 651 const G4ParticleDefinition *aParticle, 664 const G4ParticleDefinition *aParticle, 652 G4double KineticEnergy, 665 G4double KineticEnergy, 653 const G4Material *aMaterial) 666 const G4Material *aMaterial) 654 { 667 { 655 if (!t) t = new G4EnergyLossTablesHelper; 668 if (!t) t = new G4EnergyLossTablesHelper; 656 669 657 CPRWarning(); 670 CPRWarning(); 658 if( aParticle != (const G4ParticleDefinition 671 if( aParticle != (const G4ParticleDefinition*) lastParticle) 659 { 672 { 660 *t= GetTables(aParticle); 673 *t= GetTables(aParticle); 661 lastParticle = (G4ParticleDefinition*) aPa 674 lastParticle = (G4ParticleDefinition*) aParticle; 662 Chargesquare = (aParticle->GetPDGCharge()) 675 Chargesquare = (aParticle->GetPDGCharge())* 663 (aParticle->GetPDGCharge()) 676 (aParticle->GetPDGCharge())/ 664 QQPositron ; 677 QQPositron ; 665 oldIndex = -1 ; 678 oldIndex = -1 ; 666 } 679 } 667 const G4PhysicsTable* dEdxTable= t->theDEDX 680 const G4PhysicsTable* dEdxTable= t->theDEDXTable; 668 if (!dEdxTable) { 681 if (!dEdxTable) { 669 ParticleHaveNoLoss(aParticle,"dEdx"); 682 ParticleHaveNoLoss(aParticle,"dEdx"); 670 return 0.0; 683 return 0.0; 671 } 684 } 672 685 673 G4int materialIndex = (G4int)aMaterial->GetI << 686 G4int materialIndex = aMaterial->GetIndex(); 674 G4double scaledKineticEnergy = KineticEnergy 687 G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio; 675 G4double dEdx; 688 G4double dEdx; 676 G4bool isOut; 689 G4bool isOut; 677 690 678 if (scaledKineticEnergy<t->theLowestKineticE 691 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 679 692 680 dEdx = std::sqrt(scaledKineticEnergy/t->t 693 dEdx = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy) 681 *(*dEdxTable)(materialIndex)->GetV 694 *(*dEdxTable)(materialIndex)->GetValue( 682 t->theLowestKineticEnergy,isOut) 695 t->theLowestKineticEnergy,isOut); 683 696 684 } else if (scaledKineticEnergy>t->theHighest 697 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 685 698 686 dEdx = (*dEdxTable)(materialIndex)->GetVa 699 dEdx = (*dEdxTable)(materialIndex)->GetValue( 687 t->theHighestKineticEnergy,isOut); 700 t->theHighestKineticEnergy,isOut); 688 701 689 } else { 702 } else { 690 703 691 dEdx = (*dEdxTable)(materialIndex)->GetV 704 dEdx = (*dEdxTable)(materialIndex)->GetValue( 692 scaledKineticEnergy, 705 scaledKineticEnergy,isOut) ; 693 706 694 } 707 } 695 708 696 return dEdx*Chargesquare; 709 return dEdx*Chargesquare; 697 } 710 } 698 711 699 //....oooOO0OOooo........oooOO0OOooo........oo 712 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 700 713 701 G4double G4EnergyLossTables::GetPreciseRangeF 714 G4double G4EnergyLossTables::GetPreciseRangeFromEnergy( 702 const G4ParticleDefinition *aParticle, 715 const G4ParticleDefinition *aParticle, 703 G4double KineticEnergy, 716 G4double KineticEnergy, 704 const G4Material *aMaterial) 717 const G4Material *aMaterial) 705 { 718 { 706 if (!t) t = new G4EnergyLossTablesHelper; 719 if (!t) t = new G4EnergyLossTablesHelper; 707 720 708 CPRWarning(); 721 CPRWarning(); 709 if( aParticle != (const G4ParticleDefinition 722 if( aParticle != (const G4ParticleDefinition*) lastParticle) 710 { 723 { 711 *t= GetTables(aParticle); 724 *t= GetTables(aParticle); 712 lastParticle = (G4ParticleDefinition*) aPa 725 lastParticle = (G4ParticleDefinition*) aParticle; 713 Chargesquare = (aParticle->GetPDGCharge()) 726 Chargesquare = (aParticle->GetPDGCharge())* 714 (aParticle->GetPDGCharge()) 727 (aParticle->GetPDGCharge())/ 715 QQPositron ; 728 QQPositron ; 716 oldIndex = -1 ; 729 oldIndex = -1 ; 717 } 730 } 718 const G4PhysicsTable* rangeTable= t->theRang 731 const G4PhysicsTable* rangeTable= t->theRangeTable; 719 const G4PhysicsTable* dEdxTable= t->theDEDX 732 const G4PhysicsTable* dEdxTable= t->theDEDXTable; 720 if (!rangeTable) { 733 if (!rangeTable) { 721 ParticleHaveNoLoss(aParticle,"Range"); 734 ParticleHaveNoLoss(aParticle,"Range"); 722 return 0.0; 735 return 0.0; 723 } 736 } 724 G4int materialIndex = (G4int)aMaterial->GetI << 737 G4int materialIndex = aMaterial->GetIndex(); 725 738 726 G4double Thighr = t->theHighestKineticEnergy 739 G4double Thighr = t->theHighestKineticEnergy*t->theLowestKineticEnergy/ 727 (*rangeTable)(materialIndex 740 (*rangeTable)(materialIndex)-> 728 GetLowEdgeEnergy(1) ; 741 GetLowEdgeEnergy(1) ; 729 742 730 G4double scaledKineticEnergy = KineticEnergy 743 G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio; 731 G4double Range; 744 G4double Range; 732 G4bool isOut; 745 G4bool isOut; 733 746 734 if (scaledKineticEnergy<t->theLowestKineticE 747 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 735 748 736 Range = std::sqrt(scaledKineticEnergy/t->t 749 Range = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy)* 737 (*rangeTable)(materialIndex)->GetV 750 (*rangeTable)(materialIndex)->GetValue( 738 t->theLowestKineticEnergy,isOut) 751 t->theLowestKineticEnergy,isOut); 739 752 740 } else if (scaledKineticEnergy>Thighr) { 753 } else if (scaledKineticEnergy>Thighr) { 741 754 742 Range = (*rangeTable)(materialIndex)->GetV 755 Range = (*rangeTable)(materialIndex)->GetValue( 743 Thighr,isOut)+ 756 Thighr,isOut)+ 744 (scaledKineticEnergy-Thighr)/ 757 (scaledKineticEnergy-Thighr)/ 745 (*dEdxTable)(materialIndex)->GetVa 758 (*dEdxTable)(materialIndex)->GetValue( 746 Thighr,isOut); 759 Thighr,isOut); 747 760 748 } else { 761 } else { 749 762 750 Range = (*rangeTable)(materialIndex)->Get 763 Range = (*rangeTable)(materialIndex)->GetValue( 751 scaledKineticEnergy,isO 764 scaledKineticEnergy,isOut) ; 752 765 753 } 766 } 754 767 755 return Range/(Chargesquare*t->theMassRatio); 768 return Range/(Chargesquare*t->theMassRatio); 756 } 769 } 757 770 758 //....oooOO0OOooo........oooOO0OOooo........oo 771 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 759 772 760 G4double G4EnergyLossTables::GetDEDX( 773 G4double G4EnergyLossTables::GetDEDX( 761 const G4ParticleDefinition *aParticle, 774 const G4ParticleDefinition *aParticle, 762 G4double KineticEnergy, 775 G4double KineticEnergy, 763 const G4MaterialCutsCouple *couple, 776 const G4MaterialCutsCouple *couple, 764 G4bool check) 777 G4bool check) 765 { 778 { 766 if (!t) t = new G4EnergyLossTablesHelper; 779 if (!t) t = new G4EnergyLossTablesHelper; 767 780 768 if(aParticle != (const G4ParticleDefinition* 781 if(aParticle != (const G4ParticleDefinition*) lastParticle) 769 { 782 { 770 *t= GetTables(aParticle); 783 *t= GetTables(aParticle); 771 lastParticle = (G4ParticleDefinition*) aPa 784 lastParticle = (G4ParticleDefinition*) aParticle ; 772 Chargesquare = (aParticle->GetPDGCharge()) 785 Chargesquare = (aParticle->GetPDGCharge())* 773 (aParticle->GetPDGCharge()) 786 (aParticle->GetPDGCharge())/ 774 QQPositron ; 787 QQPositron ; 775 oldIndex = -1 ; 788 oldIndex = -1 ; 776 } 789 } 777 const G4PhysicsTable* dEdxTable= t->theDEDX 790 const G4PhysicsTable* dEdxTable= t->theDEDXTable; 778 791 779 if (!dEdxTable ) { 792 if (!dEdxTable ) { 780 if (check) return G4LossTableManager::Inst 793 if (check) return G4LossTableManager::Instance()->GetDEDX(aParticle,KineticEnergy,couple); 781 else ParticleHaveNoLoss(aParticle, " 794 else ParticleHaveNoLoss(aParticle, "dEdx"); 782 return 0.0; 795 return 0.0; 783 } 796 } 784 797 785 G4int materialIndex = couple->GetIndex(); 798 G4int materialIndex = couple->GetIndex(); 786 G4double scaledKineticEnergy = KineticEnergy 799 G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio; 787 G4double dEdx; 800 G4double dEdx; 788 G4bool isOut; 801 G4bool isOut; 789 802 790 if (scaledKineticEnergy<t->theLowestKineticE 803 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 791 804 792 dEdx =(*dEdxTable)(materialIndex)->GetVal 805 dEdx =(*dEdxTable)(materialIndex)->GetValue( 793 t->theLowestKineticEnergy,isOut) 806 t->theLowestKineticEnergy,isOut) 794 *std::sqrt(scaledKineticEnergy/t->t 807 *std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy); 795 808 796 } else if (scaledKineticEnergy>t->theHighest 809 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 797 810 798 dEdx = (*dEdxTable)(materialIndex)->GetVa 811 dEdx = (*dEdxTable)(materialIndex)->GetValue( 799 t->theHighestKineticEnergy,isOut); 812 t->theHighestKineticEnergy,isOut); 800 813 801 } else { 814 } else { 802 815 803 dEdx = (*dEdxTable)(materialIndex)->GetVal 816 dEdx = (*dEdxTable)(materialIndex)->GetValue( 804 scaledKineticEnergy,isOut); 817 scaledKineticEnergy,isOut); 805 818 806 } 819 } 807 820 808 return dEdx*Chargesquare; 821 return dEdx*Chargesquare; 809 } 822 } 810 823 811 //....oooOO0OOooo........oooOO0OOooo........oo 824 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 812 825 813 G4double G4EnergyLossTables::GetRange( 826 G4double G4EnergyLossTables::GetRange( 814 const G4ParticleDefinition *aParticle, 827 const G4ParticleDefinition *aParticle, 815 G4double KineticEnergy, 828 G4double KineticEnergy, 816 const G4MaterialCutsCouple *couple, 829 const G4MaterialCutsCouple *couple, 817 G4bool check) 830 G4bool check) 818 { 831 { 819 if (!t) t = new G4EnergyLossTablesHelper; 832 if (!t) t = new G4EnergyLossTablesHelper; 820 833 821 if(aParticle != (const G4ParticleDefinition* 834 if(aParticle != (const G4ParticleDefinition*) lastParticle) 822 { 835 { 823 *t= GetTables(aParticle); 836 *t= GetTables(aParticle); 824 lastParticle = (G4ParticleDefinition*) aPa 837 lastParticle = (G4ParticleDefinition*) aParticle ; 825 Chargesquare = (aParticle->GetPDGCharge()) 838 Chargesquare = (aParticle->GetPDGCharge())* 826 (aParticle->GetPDGCharge()) 839 (aParticle->GetPDGCharge())/ 827 QQPositron ; 840 QQPositron ; 828 oldIndex = -1 ; 841 oldIndex = -1 ; 829 } 842 } 830 const G4PhysicsTable* rangeTable= t->theRang 843 const G4PhysicsTable* rangeTable= t->theRangeTable; 831 const G4PhysicsTable* dEdxTable= t->theDEDX 844 const G4PhysicsTable* dEdxTable= t->theDEDXTable; 832 if (!rangeTable) { 845 if (!rangeTable) { 833 if(check) return G4LossTableManager::Insta 846 if(check) return G4LossTableManager::Instance()->GetRange(aParticle,KineticEnergy,couple); 834 else return DBL_MAX; 847 else return DBL_MAX; 835 //ParticleHaveNoLoss(aParticle,"Range"); 848 //ParticleHaveNoLoss(aParticle,"Range"); 836 } 849 } 837 850 838 G4int materialIndex = couple->GetIndex(); 851 G4int materialIndex = couple->GetIndex(); 839 G4double scaledKineticEnergy = KineticEnergy 852 G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio; 840 G4double Range; 853 G4double Range; 841 G4bool isOut; 854 G4bool isOut; 842 855 843 if (scaledKineticEnergy<t->theLowestKineticE 856 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 844 857 845 Range = std::sqrt(scaledKineticEnergy/t->t 858 Range = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy)* 846 (*rangeTable)(materialIndex)->GetV 859 (*rangeTable)(materialIndex)->GetValue( 847 t->theLowestKineticEnergy,isOut) 860 t->theLowestKineticEnergy,isOut); 848 861 849 } else if (scaledKineticEnergy>t->theHighest 862 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 850 863 851 Range = (*rangeTable)(materialIndex)->GetV 864 Range = (*rangeTable)(materialIndex)->GetValue( 852 t->theHighestKineticEnergy,isOut)+ 865 t->theHighestKineticEnergy,isOut)+ 853 (scaledKineticEnergy-t->theHighest 866 (scaledKineticEnergy-t->theHighestKineticEnergy)/ 854 (*dEdxTable)(materialIndex)->GetVa 867 (*dEdxTable)(materialIndex)->GetValue( 855 t->theHighestKineticEnergy,isOut 868 t->theHighestKineticEnergy,isOut); 856 869 857 } else { 870 } else { 858 871 859 Range = (*rangeTable)(materialIndex)->GetV 872 Range = (*rangeTable)(materialIndex)->GetValue( 860 scaledKineticEnergy,isOut); 873 scaledKineticEnergy,isOut); 861 874 862 } 875 } 863 876 864 return Range/(Chargesquare*t->theMassRatio); 877 return Range/(Chargesquare*t->theMassRatio); 865 } 878 } 866 879 867 //....oooOO0OOooo........oooOO0OOooo........oo 880 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 868 881 869 G4double G4EnergyLossTables::GetPreciseEnergyF 882 G4double G4EnergyLossTables::GetPreciseEnergyFromRange( 870 const G4P 883 const G4ParticleDefinition *aParticle, 871 G4d 884 G4double range, 872 const G4M 885 const G4MaterialCutsCouple *couple, 873 G4bool check) 886 G4bool check) 874 // it returns the value of the kinetic energy 887 // it returns the value of the kinetic energy for a given range 875 { 888 { 876 if (!t) t = new G4EnergyLossTablesHelper; 889 if (!t) t = new G4EnergyLossTablesHelper; 877 890 878 if( aParticle != (const G4ParticleDefinition 891 if( aParticle != (const G4ParticleDefinition*) lastParticle) 879 { 892 { 880 *t= GetTables(aParticle); 893 *t= GetTables(aParticle); 881 lastParticle = (G4ParticleDefinition*) aPa 894 lastParticle = (G4ParticleDefinition*) aParticle; 882 Chargesquare = (aParticle->GetPDGCharge()) 895 Chargesquare = (aParticle->GetPDGCharge())* 883 (aParticle->GetPDGCharge()) 896 (aParticle->GetPDGCharge())/ 884 QQPositron ; 897 QQPositron ; 885 oldIndex = -1 ; 898 oldIndex = -1 ; 886 } 899 } 887 const G4PhysicsTable* dEdxTable= t->theDEDX 900 const G4PhysicsTable* dEdxTable= t->theDEDXTable; 888 const G4PhysicsTable* inverseRangeTable= t- 901 const G4PhysicsTable* inverseRangeTable= t->theInverseRangeTable; 889 902 890 if (!inverseRangeTable) { 903 if (!inverseRangeTable) { 891 if(check) return G4LossTableManager::Insta 904 if(check) return G4LossTableManager::Instance()->GetEnergy(aParticle,range,couple); 892 else return DBL_MAX; 905 else return DBL_MAX; 893 // else ParticleHaveNoLoss(aPartic 906 // else ParticleHaveNoLoss(aParticle,"InverseRange"); 894 } 907 } 895 908 896 G4double scaledrange,scaledKineticEnergy ; 909 G4double scaledrange,scaledKineticEnergy ; 897 G4bool isOut ; 910 G4bool isOut ; 898 911 899 G4int materialIndex = couple->GetIndex() ; 912 G4int materialIndex = couple->GetIndex() ; 900 913 901 if(materialIndex != oldIndex) 914 if(materialIndex != oldIndex) 902 { 915 { 903 oldIndex = materialIndex ; 916 oldIndex = materialIndex ; 904 rmin = (*inverseRangeTable)(materialIndex) 917 rmin = (*inverseRangeTable)(materialIndex)-> 905 GetLowEdgeEnergy 918 GetLowEdgeEnergy(0) ; 906 rmax = (*inverseRangeTable)(materialIndex) 919 rmax = (*inverseRangeTable)(materialIndex)-> 907 GetLowEdgeEnergy(t->theNumb 920 GetLowEdgeEnergy(t->theNumberOfBins-2) ; 908 Thigh = (*inverseRangeTable)(materialIndex 921 Thigh = (*inverseRangeTable)(materialIndex)-> 909 GetValue(rmax,is 922 GetValue(rmax,isOut) ; 910 } 923 } 911 924 912 scaledrange = range*Chargesquare*t->theMassR 925 scaledrange = range*Chargesquare*t->theMassRatio ; 913 926 914 if(scaledrange < rmin) 927 if(scaledrange < rmin) 915 { 928 { 916 scaledKineticEnergy = t->theLowestKineticE 929 scaledKineticEnergy = t->theLowestKineticEnergy* 917 scaledrange*scaledrange/(rm 930 scaledrange*scaledrange/(rmin*rmin) ; 918 } 931 } 919 else 932 else 920 { 933 { 921 if(scaledrange < rmax) 934 if(scaledrange < rmax) 922 { 935 { 923 scaledKineticEnergy = (*inverseRangeTabl 936 scaledKineticEnergy = (*inverseRangeTable)(materialIndex)-> 924 GetValue( scaled 937 GetValue( scaledrange,isOut) ; 925 } 938 } 926 else 939 else 927 { 940 { 928 scaledKineticEnergy = Thigh + 941 scaledKineticEnergy = Thigh + 929 (scaledrange-rmax)* 942 (scaledrange-rmax)* 930 (*dEdxTable)(materialInd 943 (*dEdxTable)(materialIndex)-> 931 GetValue(Thig 944 GetValue(Thigh,isOut) ; 932 } 945 } 933 } 946 } 934 947 935 return scaledKineticEnergy/t->theMassRatio ; 948 return scaledKineticEnergy/t->theMassRatio ; 936 } 949 } 937 950 938 //....oooOO0OOooo........oooOO0OOooo........oo 951 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 939 952 940 G4double G4EnergyLossTables::GetPreciseDEDX( 953 G4double G4EnergyLossTables::GetPreciseDEDX( 941 const G4ParticleDefinition *aParticle, 954 const G4ParticleDefinition *aParticle, 942 G4double KineticEnergy, 955 G4double KineticEnergy, 943 const G4MaterialCutsCouple *couple) 956 const G4MaterialCutsCouple *couple) 944 { 957 { 945 if (!t) t = new G4EnergyLossTablesHelper; 958 if (!t) t = new G4EnergyLossTablesHelper; 946 959 947 if( aParticle != (const G4ParticleDefinition 960 if( aParticle != (const G4ParticleDefinition*) lastParticle) 948 { 961 { 949 *t= GetTables(aParticle); 962 *t= GetTables(aParticle); 950 lastParticle = (G4ParticleDefinition*) aPa 963 lastParticle = (G4ParticleDefinition*) aParticle; 951 Chargesquare = (aParticle->GetPDGCharge()) 964 Chargesquare = (aParticle->GetPDGCharge())* 952 (aParticle->GetPDGCharge()) 965 (aParticle->GetPDGCharge())/ 953 QQPositron ; 966 QQPositron ; 954 oldIndex = -1 ; 967 oldIndex = -1 ; 955 } 968 } 956 const G4PhysicsTable* dEdxTable= t->theDEDX 969 const G4PhysicsTable* dEdxTable= t->theDEDXTable; 957 if ( !dEdxTable ) 970 if ( !dEdxTable ) 958 return G4LossTableManager::Instance()->Get 971 return G4LossTableManager::Instance()->GetDEDX(aParticle,KineticEnergy,couple); 959 972 960 G4int materialIndex = couple->GetIndex(); 973 G4int materialIndex = couple->GetIndex(); 961 G4double scaledKineticEnergy = KineticEnergy 974 G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio; 962 G4double dEdx; 975 G4double dEdx; 963 G4bool isOut; 976 G4bool isOut; 964 977 965 if (scaledKineticEnergy<t->theLowestKineticE 978 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 966 979 967 dEdx = std::sqrt(scaledKineticEnergy/t->t 980 dEdx = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy) 968 *(*dEdxTable)(materialIndex)->GetV 981 *(*dEdxTable)(materialIndex)->GetValue( 969 t->theLowestKineticEnergy,isOut) 982 t->theLowestKineticEnergy,isOut); 970 983 971 } else if (scaledKineticEnergy>t->theHighest 984 } else if (scaledKineticEnergy>t->theHighestKineticEnergy) { 972 985 973 dEdx = (*dEdxTable)(materialIndex)->GetVa 986 dEdx = (*dEdxTable)(materialIndex)->GetValue( 974 t->theHighestKineticEnergy,isOut); 987 t->theHighestKineticEnergy,isOut); 975 988 976 } else { 989 } else { 977 990 978 dEdx = (*dEdxTable)(materialIndex)->GetV 991 dEdx = (*dEdxTable)(materialIndex)->GetValue( 979 scaledKineticEnergy, 992 scaledKineticEnergy,isOut) ; 980 993 981 } 994 } 982 995 983 return dEdx*Chargesquare; 996 return dEdx*Chargesquare; 984 } 997 } 985 998 986 //....oooOO0OOooo........oooOO0OOooo........oo 999 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 987 1000 988 G4double G4EnergyLossTables::GetPreciseRangeFr 1001 G4double G4EnergyLossTables::GetPreciseRangeFromEnergy( 989 const G4ParticleDefinition *aParticle, 1002 const G4ParticleDefinition *aParticle, 990 G4double KineticEnergy, 1003 G4double KineticEnergy, 991 const G4MaterialCutsCouple *couple) 1004 const G4MaterialCutsCouple *couple) 992 { 1005 { 993 if (!t) t = new G4EnergyLossTablesHelper; 1006 if (!t) t = new G4EnergyLossTablesHelper; 994 1007 995 if( aParticle != (const G4ParticleDefinition 1008 if( aParticle != (const G4ParticleDefinition*) lastParticle) 996 { 1009 { 997 *t= GetTables(aParticle); 1010 *t= GetTables(aParticle); 998 lastParticle = (G4ParticleDefinition*) aPa 1011 lastParticle = (G4ParticleDefinition*) aParticle; 999 Chargesquare = (aParticle->GetPDGCharge()) 1012 Chargesquare = (aParticle->GetPDGCharge())* 1000 (aParticle->GetPDGCharge() 1013 (aParticle->GetPDGCharge())/ 1001 QQPositron ; 1014 QQPositron ; 1002 oldIndex = -1 ; 1015 oldIndex = -1 ; 1003 } 1016 } 1004 const G4PhysicsTable* rangeTable= t->theRan 1017 const G4PhysicsTable* rangeTable= t->theRangeTable; 1005 const G4PhysicsTable* dEdxTable= t->theDED 1018 const G4PhysicsTable* dEdxTable= t->theDEDXTable; 1006 if ( !dEdxTable || !rangeTable) 1019 if ( !dEdxTable || !rangeTable) 1007 return G4LossTableManager::Instance()->Ge 1020 return G4LossTableManager::Instance()->GetDEDX(aParticle,KineticEnergy,couple); 1008 1021 1009 G4int materialIndex = couple->GetIndex(); 1022 G4int materialIndex = couple->GetIndex(); 1010 1023 1011 G4double Thighr = t->theHighestKineticEnerg 1024 G4double Thighr = t->theHighestKineticEnergy*t->theLowestKineticEnergy/ 1012 (*rangeTable)(materialInde 1025 (*rangeTable)(materialIndex)-> 1013 GetLowEdgeEnergy(1) ; 1026 GetLowEdgeEnergy(1) ; 1014 1027 1015 G4double scaledKineticEnergy = KineticEnerg 1028 G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio; 1016 G4double Range; 1029 G4double Range; 1017 G4bool isOut; 1030 G4bool isOut; 1018 1031 1019 if (scaledKineticEnergy<t->theLowestKinetic 1032 if (scaledKineticEnergy<t->theLowestKineticEnergy) { 1020 1033 1021 Range = std::sqrt(scaledKineticEnergy/t-> 1034 Range = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy)* 1022 (*rangeTable)(materialIndex)->Get 1035 (*rangeTable)(materialIndex)->GetValue( 1023 t->theLowestKineticEnergy,isOut 1036 t->theLowestKineticEnergy,isOut); 1024 1037 1025 } else if (scaledKineticEnergy>Thighr) { 1038 } else if (scaledKineticEnergy>Thighr) { 1026 1039 1027 Range = (*rangeTable)(materialIndex)->Get 1040 Range = (*rangeTable)(materialIndex)->GetValue( 1028 Thighr,isOut)+ 1041 Thighr,isOut)+ 1029 (scaledKineticEnergy-Thighr)/ 1042 (scaledKineticEnergy-Thighr)/ 1030 (*dEdxTable)(materialIndex)->GetV 1043 (*dEdxTable)(materialIndex)->GetValue( 1031 Thighr,isOut); 1044 Thighr,isOut); 1032 1045 1033 } else { 1046 } else { 1034 1047 1035 Range = (*rangeTable)(materialIndex)->Ge 1048 Range = (*rangeTable)(materialIndex)->GetValue( 1036 scaledKineticEnergy,is 1049 scaledKineticEnergy,isOut) ; 1037 1050 1038 } 1051 } 1039 1052 1040 return Range/(Chargesquare*t->theMassRatio) 1053 return Range/(Chargesquare*t->theMassRatio); 1041 } 1054 } 1042 1055 1043 //....oooOO0OOooo........oooOO0OOooo........o 1056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1044 1057 1045 void G4EnergyLossTables::CPRWarning() 1058 void G4EnergyLossTables::CPRWarning() 1046 { 1059 { 1047 if (let_counter < let_max_num_warnings) { 1060 if (let_counter < let_max_num_warnings) { 1048 1061 1049 G4cout << G4endl; 1062 G4cout << G4endl; 1050 G4cout << "##### G4EnergyLossTable WARNIN 1063 G4cout << "##### G4EnergyLossTable WARNING: The obsolete interface is used!" << G4endl; 1051 G4cout << "##### RESULTS ARE NOT GARANTEE 1064 G4cout << "##### RESULTS ARE NOT GARANTEED!" << G4endl; 1052 G4cout << "##### Please, substitute G4Mat 1065 G4cout << "##### Please, substitute G4Material by G4MaterialCutsCouple" << G4endl; 1053 G4cout << "##### Obsolete interface will 1066 G4cout << "##### Obsolete interface will be removed soon" << G4endl; 1054 G4cout << G4endl; 1067 G4cout << G4endl; 1055 let_counter++; 1068 let_counter++; 1056 1069 1057 } else if (let_counter == let_max_num_warni 1070 } else if (let_counter == let_max_num_warnings) { 1058 1071 1059 G4cout << "##### G4EnergyLossTable WARNIN 1072 G4cout << "##### G4EnergyLossTable WARNING closed" << G4endl; 1060 let_counter++; 1073 let_counter++; 1061 } 1074 } 1062 } 1075 } 1063 1076 1064 //....oooOO0OOooo........oooOO0OOooo........o 1077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1065 1078 1066 void 1079 void 1067 G4EnergyLossTables::ParticleHaveNoLoss(const 1080 G4EnergyLossTables::ParticleHaveNoLoss(const G4ParticleDefinition*, 1068 const G4String& /*q*/) 1081 const G4String& /*q*/) 1069 { 1082 { 1070 } 1083 } 1071 1084 1072 //....oooOO0OOooo........oooOO0OOooo........o 1085 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1073 1086