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 source file 29 // GEANT4 class source file 30 // 30 // 31 // Class: G4IonDEDXHandler 31 // Class: G4IonDEDXHandler 32 // 32 // 33 // Author: Anton Lechner (Anton. 33 // Author: Anton Lechner (Anton.Lechner@cern.ch) 34 // 34 // 35 // First implementation: 11. 03. 2009 35 // First implementation: 11. 03. 2009 36 // 36 // 37 // Modifications: 12. 11 .2009 - Function Buil 37 // Modifications: 12. 11 .2009 - Function BuildDEDXTable: Using adapted build 38 // methods of st 38 // methods of stopping power classes according 39 // to interface 39 // to interface change in G4VIonDEDXTable. 40 // Function Upda 40 // Function UpdateCacheValue: Using adapted 41 // ScalingFactor 41 // ScalingFactorEnergy function according to 42 // interface cha 42 // interface change in G4VIonDEDXScaling- 43 // Algorithm (AL 43 // Algorithm (AL) 44 // 44 // 45 // Class description: 45 // Class description: 46 // Ion dE/dx table handler. 46 // Ion dE/dx table handler. 47 // 47 // 48 // Comments: 48 // Comments: 49 // 49 // 50 // =========================================== 50 // =========================================================================== 51 51 52 #include <iomanip> 52 #include <iomanip> 53 53 54 #include "G4IonDEDXHandler.hh" 54 #include "G4IonDEDXHandler.hh" 55 #include "G4SystemOfUnits.hh" 55 #include "G4SystemOfUnits.hh" 56 #include "G4VIonDEDXTable.hh" 56 #include "G4VIonDEDXTable.hh" 57 #include "G4VIonDEDXScalingAlgorithm.hh" 57 #include "G4VIonDEDXScalingAlgorithm.hh" 58 #include "G4ParticleDefinition.hh" 58 #include "G4ParticleDefinition.hh" 59 #include "G4Material.hh" 59 #include "G4Material.hh" 60 #include "G4PhysicsFreeVector.hh" << 60 #include "G4LPhysicsFreeVector.hh" 61 #include "G4Exp.hh" 61 #include "G4Exp.hh" 62 62 >> 63 //#define PRINT_DEBUG >> 64 >> 65 63 // ########################################### 66 // ######################################################################### 64 67 65 G4IonDEDXHandler::G4IonDEDXHandler( 68 G4IonDEDXHandler::G4IonDEDXHandler( 66 G4VIonDEDXTable* ionTable, << 69 G4VIonDEDXTable* ionTable, 67 G4VIonDEDXScalingAlgorithm* ionAlgo << 70 G4VIonDEDXScalingAlgorithm* ionAlgorithm, 68 const G4String& name, << 71 const G4String& name, 69 G4int maxCacheSize, << 72 G4int maxCacheSize, 70 G4bool splines) : << 73 G4bool splines) : 71 table(ionTable), 74 table(ionTable), 72 algorithm(ionAlgorithm), 75 algorithm(ionAlgorithm), 73 tableName(name), 76 tableName(name), 74 useSplines(splines), 77 useSplines(splines), 75 maxCacheEntries(maxCacheSize) { 78 maxCacheEntries(maxCacheSize) { 76 79 77 if(table == nullptr) { << 80 if(table == 0) { 78 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan 81 G4cerr << "G4IonDEDXHandler::G4IonDEDXHandler() " 79 << " Pointer to G4VIonDEDXTable ob 82 << " Pointer to G4VIonDEDXTable object is null-pointer." 80 << G4endl; 83 << G4endl; 81 } 84 } 82 85 83 if(algorithm == nullptr) { << 86 if(algorithm == 0) { 84 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan 87 G4cerr << "G4IonDEDXHandler::G4IonDEDXHandler() " 85 << " Pointer to G4VIonDEDXScalingA 88 << " Pointer to G4VIonDEDXScalingAlgorithm object is null-pointer." 86 << G4endl; 89 << G4endl; 87 } 90 } 88 91 89 if(maxCacheEntries <= 0) { 92 if(maxCacheEntries <= 0) { 90 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan 93 G4cerr << "G4IonDEDXHandler::G4IonDEDXHandler() " 91 << " Cache size <=0. Resetting to 94 << " Cache size <=0. Resetting to 5." 92 << G4endl; 95 << G4endl; 93 maxCacheEntries = 5; 96 maxCacheEntries = 5; 94 } 97 } 95 } 98 } 96 99 97 // ########################################### 100 // ######################################################################### 98 101 99 G4IonDEDXHandler::~G4IonDEDXHandler() { 102 G4IonDEDXHandler::~G4IonDEDXHandler() { 100 103 101 ClearCache(); 104 ClearCache(); 102 105 103 // All stopping power vectors built accordin 106 // All stopping power vectors built according to Bragg's addivitiy rule 104 // are deleted. All other stopping power vec 107 // are deleted. All other stopping power vectors are expected to be 105 // deleted by their creator class (sub-class 108 // deleted by their creator class (sub-class of G4VIonDEDXTable). 106 // DEDXTableBraggRule::iterator iter = stopp 109 // DEDXTableBraggRule::iterator iter = stoppingPowerTableBragg.begin(); 107 // DEDXTableBraggRule::iterator iter_end = s 110 // DEDXTableBraggRule::iterator iter_end = stoppingPowerTableBragg.end(); 108 111 109 // for(;iter != iter_end; iter++) delete it 112 // for(;iter != iter_end; iter++) delete iter -> second; 110 stoppingPowerTableBragg.clear(); 113 stoppingPowerTableBragg.clear(); 111 114 112 stoppingPowerTable.clear(); 115 stoppingPowerTable.clear(); 113 116 114 if(table != nullptr) << 117 if(table != 0) delete table; 115 delete table; << 118 if(algorithm != 0) delete algorithm; 116 if(algorithm != nullptr) << 117 delete algorithm; << 118 } 119 } 119 120 120 // ########################################### 121 // ######################################################################### 121 122 122 G4bool G4IonDEDXHandler::IsApplicable( 123 G4bool G4IonDEDXHandler::IsApplicable( 123 const G4ParticleDefinition* p 124 const G4ParticleDefinition* particle, // Projectile (ion) 124 const G4Material* material) { 125 const G4Material* material) { // Target material 125 126 126 G4bool isApplicable = true; 127 G4bool isApplicable = true; 127 128 128 if(table == nullptr || algorithm == nullptr) << 129 if(table == 0 || algorithm == 0) { 129 isApplicable = false; 130 isApplicable = false; 130 } 131 } 131 else { 132 else { 132 133 133 G4int atomicNumberIon = particle -> GetAt 134 G4int atomicNumberIon = particle -> GetAtomicNumber(); 134 G4int atomicNumberBase = 135 G4int atomicNumberBase = 135 algorithm -> AtomicNumberBaseI 136 algorithm -> AtomicNumberBaseIon(atomicNumberIon, material); 136 137 137 G4IonKey key = std::make_pair(atomicNumbe 138 G4IonKey key = std::make_pair(atomicNumberBase, material); 138 139 139 DEDXTable::iterator iter = stoppingPowerT 140 DEDXTable::iterator iter = stoppingPowerTable.find(key); 140 if(iter == stoppingPowerTable.end()) isAp 141 if(iter == stoppingPowerTable.end()) isApplicable = false; 141 } 142 } 142 143 143 return isApplicable; 144 return isApplicable; 144 } 145 } 145 146 146 // ########################################### 147 // ######################################################################### 147 148 148 G4double G4IonDEDXHandler::GetDEDX( 149 G4double G4IonDEDXHandler::GetDEDX( 149 const G4ParticleDefinition* p 150 const G4ParticleDefinition* particle, // Projectile (ion) 150 const G4Material* material, 151 const G4Material* material, // Target material 151 G4double kineticEnergy) { 152 G4double kineticEnergy) { // Kinetic energy of projectile 152 153 153 G4double dedx = 0.0; 154 G4double dedx = 0.0; 154 155 155 G4CacheValue value = GetCacheValue(particle, 156 G4CacheValue value = GetCacheValue(particle, material); 156 157 157 if(kineticEnergy <= 0.0) dedx = 0.0; 158 if(kineticEnergy <= 0.0) dedx = 0.0; 158 else if(value.dedxVector != 0) { 159 else if(value.dedxVector != 0) { 159 160 160 G4bool b; 161 G4bool b; 161 G4double factor = value.density; 162 G4double factor = value.density; 162 163 163 factor *= algorithm -> ScalingFactorDEDX( 164 factor *= algorithm -> ScalingFactorDEDX(particle, 164 m 165 material, 165 k 166 kineticEnergy); 166 G4double scaledKineticEnergy = kineticEne 167 G4double scaledKineticEnergy = kineticEnergy * value.energyScaling; 167 168 168 if(scaledKineticEnergy < value.lowerEnerg 169 if(scaledKineticEnergy < value.lowerEnergyEdge) { 169 170 170 factor *= std::sqrt(scaledKineticEnerg 171 factor *= std::sqrt(scaledKineticEnergy / value.lowerEnergyEdge); 171 scaledKineticEnergy = value.lowerEnerg 172 scaledKineticEnergy = value.lowerEnergyEdge; 172 } 173 } 173 174 174 dedx = factor * value.dedxVector -> GetVa 175 dedx = factor * value.dedxVector -> GetValue(scaledKineticEnergy, b); 175 176 176 if(dedx < 0.0) dedx = 0.0; 177 if(dedx < 0.0) dedx = 0.0; 177 } 178 } 178 else dedx = 0.0; 179 else dedx = 0.0; 179 180 180 #ifdef PRINT_DEBUG 181 #ifdef PRINT_DEBUG 181 G4cout << "G4IonDEDXHandler::GetDEDX() E 182 G4cout << "G4IonDEDXHandler::GetDEDX() E = " 182 << kineticEnergy / MeV << " MeV * 183 << kineticEnergy / MeV << " MeV * " 183 << value.energyScaling << " = " 184 << value.energyScaling << " = " 184 << kineticEnergy * value.energySca 185 << kineticEnergy * value.energyScaling / MeV 185 << " MeV, dE/dx = " << dedx / MeV 186 << " MeV, dE/dx = " << dedx / MeV * cm << " MeV/cm" 186 << ", material = " << material -> 187 << ", material = " << material -> GetName() 187 << G4endl; 188 << G4endl; 188 #endif 189 #endif 189 190 190 return dedx; 191 return dedx; 191 } 192 } 192 193 193 // ########################################### 194 // ######################################################################### 194 195 195 G4bool G4IonDEDXHandler::BuildDEDXTable( 196 G4bool G4IonDEDXHandler::BuildDEDXTable( 196 const G4ParticleDefinition* p 197 const G4ParticleDefinition* particle, // Projectile (ion) 197 const G4Material* material) { 198 const G4Material* material) { // Target material 198 199 199 G4int atomicNumberIon = particle -> GetAtomi 200 G4int atomicNumberIon = particle -> GetAtomicNumber(); 200 201 201 G4bool isApplicable = BuildDEDXTable(atomicN 202 G4bool isApplicable = BuildDEDXTable(atomicNumberIon, material); 202 203 203 return isApplicable; 204 return isApplicable; 204 } 205 } 205 206 206 207 207 // ########################################### 208 // ######################################################################### 208 209 209 G4bool G4IonDEDXHandler::BuildDEDXTable( 210 G4bool G4IonDEDXHandler::BuildDEDXTable( 210 G4int atomicNumberIon, 211 G4int atomicNumberIon, // Projectile (ion) 211 const G4Material* material) { 212 const G4Material* material) { // Target material 212 213 213 G4bool isApplicable = true; 214 G4bool isApplicable = true; 214 215 215 if(table == 0 || algorithm == 0) { 216 if(table == 0 || algorithm == 0) { 216 isApplicable = false; 217 isApplicable = false; 217 return isApplicable; 218 return isApplicable; 218 } 219 } 219 220 220 G4int atomicNumberBase = 221 G4int atomicNumberBase = 221 algorithm -> AtomicNumberBaseI 222 algorithm -> AtomicNumberBaseIon(atomicNumberIon, material); 222 223 223 // Checking if vector is already built, and 224 // Checking if vector is already built, and returns if this is indeed 224 // the case 225 // the case 225 G4IonKey key = std::make_pair(atomicNumberBa 226 G4IonKey key = std::make_pair(atomicNumberBase, material); 226 227 227 auto iter = stoppingPowerTable.find(key); << 228 DEDXTable::iterator iter = stoppingPowerTable.find(key); 228 if(iter != stoppingPowerTable.end()) return 229 if(iter != stoppingPowerTable.end()) return isApplicable; 229 230 230 // Checking if table contains stopping power 231 // Checking if table contains stopping power vector for given material name 231 // or chemical formula 232 // or chemical formula 232 const G4String& chemFormula = material -> Ge 233 const G4String& chemFormula = material -> GetChemicalFormula(); 233 const G4String& materialName = material -> G 234 const G4String& materialName = material -> GetName(); 234 235 235 isApplicable = table -> BuildPhysicsVector(a 236 isApplicable = table -> BuildPhysicsVector(atomicNumberBase, chemFormula); 236 237 237 if(isApplicable) { 238 if(isApplicable) { 238 stoppingPowerTable[key] = 239 stoppingPowerTable[key] = 239 table -> GetPhysicsVector(atomic 240 table -> GetPhysicsVector(atomicNumberBase, chemFormula); 240 return isApplicable; 241 return isApplicable; 241 } 242 } 242 243 243 isApplicable = table -> BuildPhysicsVector(a 244 isApplicable = table -> BuildPhysicsVector(atomicNumberBase, materialName); 244 if(isApplicable) { 245 if(isApplicable) { 245 stoppingPowerTable[key] = 246 stoppingPowerTable[key] = 246 table -> GetPhysicsVector(atomic 247 table -> GetPhysicsVector(atomicNumberBase, materialName); 247 return isApplicable; 248 return isApplicable; 248 } 249 } 249 250 250 // Building the stopping power vector based 251 // Building the stopping power vector based on Bragg's additivity rule 251 const G4ElementVector* elementVector = mater 252 const G4ElementVector* elementVector = material -> GetElementVector() ; 252 253 253 std::vector<G4PhysicsVector*> dEdxTable; 254 std::vector<G4PhysicsVector*> dEdxTable; 254 255 255 size_t nmbElements = material -> GetNumberOf 256 size_t nmbElements = material -> GetNumberOfElements(); 256 257 257 for(size_t i = 0; i < nmbElements; i++) { 258 for(size_t i = 0; i < nmbElements; i++) { 258 259 259 G4int atomicNumberMat = G4int((*elementV 260 G4int atomicNumberMat = G4int((*elementVector)[i] -> GetZ()); 260 261 261 isApplicable = table -> BuildPhysicsVect 262 isApplicable = table -> BuildPhysicsVector(atomicNumberBase, atomicNumberMat); 262 263 263 if(isApplicable) { 264 if(isApplicable) { 264 265 265 G4PhysicsVector* dEdx = 266 G4PhysicsVector* dEdx = 266 table -> GetPhysicsVector(at 267 table -> GetPhysicsVector(atomicNumberBase, atomicNumberMat); 267 dEdxTable.push_back(dEdx); 268 dEdxTable.push_back(dEdx); 268 } 269 } 269 else { 270 else { 270 271 271 dEdxTable.clear(); 272 dEdxTable.clear(); 272 break; 273 break; 273 } 274 } 274 } 275 } 275 276 276 if(isApplicable) { 277 if(isApplicable) { 277 278 278 if(dEdxTable.size() > 0) { 279 if(dEdxTable.size() > 0) { 279 280 280 size_t nmbdEdxBins = dEdxTable[0] -> G 281 size_t nmbdEdxBins = dEdxTable[0] -> GetVectorLength(); 281 G4double lowerEdge = dEdxTable[0] -> G 282 G4double lowerEdge = dEdxTable[0] -> GetLowEdgeEnergy(0); 282 G4double upperEdge = dEdxTable[0] -> G 283 G4double upperEdge = dEdxTable[0] -> GetLowEdgeEnergy(nmbdEdxBins-1); 283 284 284 G4PhysicsFreeVector* dEdxBragg = << 285 G4LPhysicsFreeVector* dEdxBragg = 285 new G4PhysicsFreeVector(nmbdEdxBins, << 286 new G4LPhysicsFreeVector(nmbdEdxBins, 286 lowerEdge, << 287 lowerEdge, 287 upperEdge, << 288 upperEdge); 288 useSplines); << 289 289 290 const G4double* massFractionVector = m 290 const G4double* massFractionVector = material -> GetFractionVector(); 291 291 292 G4bool b; 292 G4bool b; 293 for(size_t j = 0; j < nmbdEdxBins; j++ 293 for(size_t j = 0; j < nmbdEdxBins; j++) { 294 294 295 G4double edge = dEdxTable[0] -> Ge 295 G4double edge = dEdxTable[0] -> GetLowEdgeEnergy(j); 296 296 297 G4double value = 0.0; 297 G4double value = 0.0; 298 for(size_t i = 0; i < nmbElements; i++ 298 for(size_t i = 0; i < nmbElements; i++) { 299 299 300 value += (dEdxTable[i] -> GetV 300 value += (dEdxTable[i] -> GetValue(edge ,b)) * 301 301 massFractionVector[i]; 302 } 302 } 303 303 304 dEdxBragg -> PutValues(j, edge, va 304 dEdxBragg -> PutValues(j, edge, value); 305 } 305 } 306 if (useSplines) << 306 dEdxBragg -> SetSpline(useSplines); 307 dEdxBragg -> FillSecondDerivatives(); << 308 307 309 #ifdef PRINT_DEBUG 308 #ifdef PRINT_DEBUG 310 G4cout << "G4IonDEDXHandler::BuildPhys 309 G4cout << "G4IonDEDXHandler::BuildPhysicsVector() for ion with Z=" 311 << atomicNumberBase << " in " 310 << atomicNumberBase << " in " 312 << material -> GetName() 311 << material -> GetName() 313 << G4endl; 312 << G4endl; 314 313 315 G4cout << *dEdxBragg; 314 G4cout << *dEdxBragg; 316 #endif 315 #endif 317 316 318 stoppingPowerTable[key] = dEdxBragg; 317 stoppingPowerTable[key] = dEdxBragg; 319 stoppingPowerTableBragg[key] = dEdxBragg; 318 stoppingPowerTableBragg[key] = dEdxBragg; 320 } 319 } 321 } 320 } 322 321 323 ClearCache(); 322 ClearCache(); 324 323 325 return isApplicable; 324 return isApplicable; 326 } 325 } 327 326 328 // ########################################### 327 // ######################################################################### 329 328 330 G4CacheValue G4IonDEDXHandler::UpdateCacheValu 329 G4CacheValue G4IonDEDXHandler::UpdateCacheValue( 331 const G4ParticleDefinition* part 330 const G4ParticleDefinition* particle, // Projectile (ion) 332 const G4Material* material) { 331 const G4Material* material) { // Target material 333 332 334 G4CacheValue value; 333 G4CacheValue value; 335 334 336 G4int atomicNumberIon = particle -> GetAtomi 335 G4int atomicNumberIon = particle -> GetAtomicNumber(); 337 G4int atomicNumberBase = 336 G4int atomicNumberBase = 338 algorithm -> AtomicNumberBaseI 337 algorithm -> AtomicNumberBaseIon(atomicNumberIon, material); 339 338 340 G4IonKey key = std::make_pair(atomicNumberBa 339 G4IonKey key = std::make_pair(atomicNumberBase, material); 341 340 342 DEDXTable::iterator iter = stoppingPowerTabl 341 DEDXTable::iterator iter = stoppingPowerTable.find(key); 343 342 344 if(iter != stoppingPowerTable.end()) { 343 if(iter != stoppingPowerTable.end()) { 345 value.dedxVector = iter -> second; 344 value.dedxVector = iter -> second; 346 345 347 G4double nmbNucleons = G4double(particle 346 G4double nmbNucleons = G4double(particle -> GetAtomicMass()); 348 value.energyScaling = 347 value.energyScaling = 349 algorithm -> ScalingFactorEnergy(pa 348 algorithm -> ScalingFactorEnergy(particle, material) / nmbNucleons; 350 349 351 size_t nmbdEdxBins = value.dedxVector -> 350 size_t nmbdEdxBins = value.dedxVector -> GetVectorLength(); 352 value.lowerEnergyEdge = value.dedxVector 351 value.lowerEnergyEdge = value.dedxVector -> GetLowEdgeEnergy(0); 353 352 354 value.upperEnergyEdge = 353 value.upperEnergyEdge = 355 value.dedxVector -> Get 354 value.dedxVector -> GetLowEdgeEnergy(nmbdEdxBins-1); 356 value.density = material -> GetDensity(); 355 value.density = material -> GetDensity(); 357 } 356 } 358 else { 357 else { 359 value.dedxVector = 0; 358 value.dedxVector = 0; 360 value.energyScaling = 0.0; 359 value.energyScaling = 0.0; 361 value.lowerEnergyEdge = 0.0; 360 value.lowerEnergyEdge = 0.0; 362 value.upperEnergyEdge = 0.0; 361 value.upperEnergyEdge = 0.0; 363 value.density = 0.0; 362 value.density = 0.0; 364 } 363 } 365 364 366 #ifdef PRINT_DEBUG 365 #ifdef PRINT_DEBUG 367 G4cout << "G4IonDEDXHandler::UpdateCacheValu 366 G4cout << "G4IonDEDXHandler::UpdateCacheValue() for " 368 << particle -> GetParticleName() << " 367 << particle -> GetParticleName() << " in " 369 << material -> GetName() 368 << material -> GetName() 370 << G4endl; 369 << G4endl; 371 #endif 370 #endif 372 371 373 return value; 372 return value; 374 } 373 } 375 374 376 // ########################################### 375 // ######################################################################### 377 376 378 G4CacheValue G4IonDEDXHandler::GetCacheValue( 377 G4CacheValue G4IonDEDXHandler::GetCacheValue( 379 const G4ParticleDefinition* part 378 const G4ParticleDefinition* particle, // Projectile (ion) 380 const G4Material* material) { 379 const G4Material* material) { // Target material 381 380 382 G4CacheKey key = std::make_pair(particle, ma 381 G4CacheKey key = std::make_pair(particle, material); 383 382 384 G4CacheEntry entry; 383 G4CacheEntry entry; 385 CacheEntryList::iterator* pointerIter = 384 CacheEntryList::iterator* pointerIter = 386 (CacheEntryList::iterator*) 385 (CacheEntryList::iterator*) cacheKeyPointers[key]; 387 386 388 if(!pointerIter) { 387 if(!pointerIter) { 389 entry.value = UpdateCacheValue(particle, 388 entry.value = UpdateCacheValue(particle, material); 390 389 391 entry.key = key; 390 entry.key = key; 392 cacheEntries.push_front(entry); 391 cacheEntries.push_front(entry); 393 392 394 CacheEntryList::iterator* pointerIter1 = 393 CacheEntryList::iterator* pointerIter1 = 395 new CacheEn 394 new CacheEntryList::iterator(); 396 *pointerIter1 = cacheEntries.begin(); 395 *pointerIter1 = cacheEntries.begin(); 397 cacheKeyPointers[key] = pointerIter1; 396 cacheKeyPointers[key] = pointerIter1; 398 397 399 if(G4int(cacheEntries.size()) > maxCache 398 if(G4int(cacheEntries.size()) > maxCacheEntries) { 400 399 401 G4CacheEntry lastEntry = cacheEntries.bac 400 G4CacheEntry lastEntry = cacheEntries.back(); 402 401 403 void* pointerIter2 = cacheKeyPointers 402 void* pointerIter2 = cacheKeyPointers[lastEntry.key]; 404 CacheEntryList::iterator* listPointer 403 CacheEntryList::iterator* listPointerIter = 405 (CacheEntryList::ite 404 (CacheEntryList::iterator*) pointerIter2; 406 405 407 cacheEntries.erase(*listPointerIter); 406 cacheEntries.erase(*listPointerIter); 408 407 409 delete listPointerIter; 408 delete listPointerIter; 410 cacheKeyPointers.erase(lastEntry.key) 409 cacheKeyPointers.erase(lastEntry.key); 411 } 410 } 412 } 411 } 413 else { 412 else { 414 entry = *(*pointerIter); 413 entry = *(*pointerIter); 415 // Cache entries are currently not re-or 414 // Cache entries are currently not re-ordered. 416 // Uncomment for activating re-ordering: 415 // Uncomment for activating re-ordering: 417 // cacheEntries.erase(*pointerIter) 416 // cacheEntries.erase(*pointerIter); 418 // cacheEntries.push_front(entry); 417 // cacheEntries.push_front(entry); 419 // *pointerIter = cacheEntries.begi 418 // *pointerIter = cacheEntries.begin(); 420 } 419 } 421 return entry.value; 420 return entry.value; 422 } 421 } 423 422 424 // ########################################### 423 // ######################################################################### 425 424 426 void G4IonDEDXHandler::ClearCache() { 425 void G4IonDEDXHandler::ClearCache() { 427 426 428 CacheIterPointerMap::iterator iter = cacheKe 427 CacheIterPointerMap::iterator iter = cacheKeyPointers.begin(); 429 CacheIterPointerMap::iterator iter_end = cac 428 CacheIterPointerMap::iterator iter_end = cacheKeyPointers.end(); 430 429 431 for(;iter != iter_end; iter++) { 430 for(;iter != iter_end; iter++) { 432 void* pointerIter = iter -> second; 431 void* pointerIter = iter -> second; 433 CacheEntryList::iterator* listPointerIte 432 CacheEntryList::iterator* listPointerIter = 434 (CacheEntryList::ite 433 (CacheEntryList::iterator*) pointerIter; 435 434 436 delete listPointerIter; 435 delete listPointerIter; 437 } 436 } 438 437 439 cacheEntries.clear(); 438 cacheEntries.clear(); 440 cacheKeyPointers.clear(); 439 cacheKeyPointers.clear(); 441 } 440 } 442 441 443 // ########################################### 442 // ######################################################################### 444 443 445 void G4IonDEDXHandler::PrintDEDXTable( 444 void G4IonDEDXHandler::PrintDEDXTable( 446 const G4ParticleDefinition* 445 const G4ParticleDefinition* particle, // Projectile (ion) 447 const G4Material* material, 446 const G4Material* material, // Target material 448 G4double lowerBoundary, 447 G4double lowerBoundary, // Minimum energy per nucleon 449 G4double upperBoundary, 448 G4double upperBoundary, // Maximum energy per nucleon 450 G4int nmbBins, 449 G4int nmbBins, // Number of bins 451 G4bool logScaleEnergy) { 450 G4bool logScaleEnergy) { // Logarithmic scaling of energy 452 451 453 G4double atomicMassNumber = particle -> GetA 452 G4double atomicMassNumber = particle -> GetAtomicMass(); 454 G4double materialDensity = material -> GetDe 453 G4double materialDensity = material -> GetDensity(); 455 454 456 G4cout << "# dE/dx table for " << particle - 455 G4cout << "# dE/dx table for " << particle -> GetParticleName() 457 << " in material " << material -> Get 456 << " in material " << material -> GetName() 458 << " of density " << materialDensity 457 << " of density " << materialDensity / g * cm3 459 << " g/cm3" 458 << " g/cm3" 460 << G4endl 459 << G4endl 461 << "# Projectile mass number A1 = " < 460 << "# Projectile mass number A1 = " << atomicMassNumber 462 << G4endl 461 << G4endl 463 << "# Energy range (per nucleon) of t 462 << "# Energy range (per nucleon) of tabulation: " 464 << GetLowerEnergyEdge(particle, mater 463 << GetLowerEnergyEdge(particle, material) / atomicMassNumber / MeV 465 << " - " 464 << " - " 466 << GetUpperEnergyEdge(particle, mater 465 << GetUpperEnergyEdge(particle, material) / atomicMassNumber / MeV 467 << " MeV" 466 << " MeV" 468 << G4endl 467 << G4endl 469 << "# ------------------------------- 468 << "# ------------------------------------------------------" 470 << G4endl; 469 << G4endl; 471 G4cout << "#" 470 G4cout << "#" 472 << std::setw(13) << std::right << "E" 471 << std::setw(13) << std::right << "E" 473 << std::setw(14) << "E/A1" 472 << std::setw(14) << "E/A1" 474 << std::setw(14) << "dE/dx" 473 << std::setw(14) << "dE/dx" 475 << std::setw(14) << "1/rho*dE/dx" 474 << std::setw(14) << "1/rho*dE/dx" 476 << G4endl; 475 << G4endl; 477 G4cout << "#" 476 G4cout << "#" 478 << std::setw(13) << std::right << "(M 477 << std::setw(13) << std::right << "(MeV)" 479 << std::setw(14) << "(MeV)" 478 << std::setw(14) << "(MeV)" 480 << std::setw(14) << "(MeV/cm)" 479 << std::setw(14) << "(MeV/cm)" 481 << std::setw(14) << "(MeV*cm2/mg)" 480 << std::setw(14) << "(MeV*cm2/mg)" 482 << G4endl 481 << G4endl 483 << "# ------------------------------- 482 << "# ------------------------------------------------------" 484 << G4endl; 483 << G4endl; 485 484 486 //G4CacheValue value = GetCacheValue(particl 485 //G4CacheValue value = GetCacheValue(particle, material); 487 486 488 G4double energyLowerBoundary = lowerBoundary 487 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber; 489 G4double energyUpperBoundary = upperBoundary 488 G4double energyUpperBoundary = upperBoundary * atomicMassNumber; 490 489 491 if(logScaleEnergy) { 490 if(logScaleEnergy) { 492 491 493 energyLowerBoundary = std::log(energyLowe 492 energyLowerBoundary = std::log(energyLowerBoundary); 494 energyUpperBoundary = std::log(energyUppe 493 energyUpperBoundary = std::log(energyUpperBoundary); 495 } 494 } 496 495 497 G4double deltaEnergy = (energyUpperBoundary 496 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) / 498 497 G4double(nmbBins); 499 498 500 G4cout.precision(6); 499 G4cout.precision(6); 501 for(int i = 0; i < nmbBins + 1; i++) { 500 for(int i = 0; i < nmbBins + 1; i++) { 502 501 503 G4double energy = energyLowerBoundary + 502 G4double energy = energyLowerBoundary + i * deltaEnergy; 504 if(logScaleEnergy) energy = G4Exp(energy 503 if(logScaleEnergy) energy = G4Exp(energy); 505 504 506 G4double loss = GetDEDX(particle, materi 505 G4double loss = GetDEDX(particle, material, energy); 507 506 508 G4cout << std::setw(14) << std::right << 507 G4cout << std::setw(14) << std::right << energy / MeV 509 << std::setw(14) << energy / atom 508 << std::setw(14) << energy / atomicMassNumber / MeV 510 << std::setw(14) << loss / MeV * 509 << std::setw(14) << loss / MeV * cm 511 << std::setw(14) << loss / materi 510 << std::setw(14) << loss / materialDensity / (MeV*cm2/(0.001*g)) 512 << G4endl; 511 << G4endl; 513 } 512 } 514 } 513 } 515 514 516 // ########################################### 515 // ######################################################################### 517 516 518 G4double G4IonDEDXHandler::GetLowerEnergyEdge( 517 G4double G4IonDEDXHandler::GetLowerEnergyEdge( 519 const G4ParticleDefinition* p 518 const G4ParticleDefinition* particle, // Projectile (ion) 520 const G4Material* material) { 519 const G4Material* material) { // Target material 521 520 522 G4double edge = 0.0; 521 G4double edge = 0.0; 523 522 524 G4CacheValue value = GetCacheValue(particle, 523 G4CacheValue value = GetCacheValue(particle, material); 525 524 526 if(value.energyScaling > 0) 525 if(value.energyScaling > 0) 527 edge = value.lowerEnergyEdge / value 526 edge = value.lowerEnergyEdge / value.energyScaling; 528 527 529 return edge; 528 return edge; 530 } 529 } 531 530 532 // ########################################### 531 // ######################################################################### 533 532 534 G4double G4IonDEDXHandler::GetUpperEnergyEdge( 533 G4double G4IonDEDXHandler::GetUpperEnergyEdge( 535 const G4ParticleDefinition* p 534 const G4ParticleDefinition* particle, // Projectile (ion) 536 const G4Material* material) { 535 const G4Material* material) { // Target material 537 536 538 G4double edge = 0.0; 537 G4double edge = 0.0; 539 538 540 G4CacheValue value = GetCacheValue(particle, 539 G4CacheValue value = GetCacheValue(particle, material); 541 540 542 if(value.energyScaling > 0) 541 if(value.energyScaling > 0) 543 edge = value.upperEnergyEdge / value 542 edge = value.upperEnergyEdge / value.energyScaling; 544 543 545 return edge; 544 return edge; 546 } 545 } 547 546 548 // ########################################### 547 // ######################################################################### 549 548 550 G4String G4IonDEDXHandler::GetName() { 549 G4String G4IonDEDXHandler::GetName() { 551 550 552 return tableName; 551 return tableName; 553 } 552 } 554 553 555 // ########################################### 554 // ######################################################################### 556 555