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 // $Id: G4IonParametrisedLossModel.cc,v 1.10 2010-11-04 12:21:48 vnivanch Exp $ >> 27 // GEANT4 tag $Name: geant4-09-04-patch-01 $ 26 // 28 // 27 // =========================================== 29 // =========================================================================== 28 // GEANT4 class source file 30 // GEANT4 class source file 29 // 31 // 30 // Class: G4IonParametrisedLoss 32 // Class: G4IonParametrisedLossModel 31 // 33 // 32 // Base class: G4VEmModel (utils) 34 // Base class: G4VEmModel (utils) 33 // << 35 // 34 // Author: Anton Lechner (Anton. 36 // Author: Anton Lechner (Anton.Lechner@cern.ch) 35 // 37 // 36 // First implementation: 10. 11. 2008 38 // First implementation: 10. 11. 2008 37 // 39 // 38 // Modifications: 03. 02. 2009 - Bug fix itera 40 // Modifications: 03. 02. 2009 - Bug fix iterators (AL) 39 // 11. 03. 2009 - Introduced ne 41 // 11. 03. 2009 - Introduced new table handler(G4IonDEDXHandler) 40 // and modified 42 // and modified method to add/remove tables 41 // (tables are n 43 // (tables are now built in init. phase), 42 // Minor bug fix 44 // Minor bug fix in ComputeDEDXPerVolume (AL) 43 // 11. 05. 2009 - Introduced sc 45 // 11. 05. 2009 - Introduced scaling algorithm for heavier ions: 44 // G4IonDEDXScal 46 // G4IonDEDXScalingICRU73 (AL) 45 // 12. 11. 2009 - Moved from or 47 // 12. 11. 2009 - Moved from original ICRU 73 classes to new 46 // class (G4IonS 48 // class (G4IonStoppingData), which is capable 47 // of reading st 49 // of reading stopping power data files stored 48 // in G4LEDATA ( 50 // in G4LEDATA (requires G4EMLOW6.8 or higher). 49 // Simultanesoul << 51 // Simultanesouly, the upper energy limit of 50 // ICRU 73 is in 52 // ICRU 73 is increased to 1 GeV/nucleon. 51 // - Removed nucle 53 // - Removed nuclear stopping from Corrections- 52 // AlongStep sin 54 // AlongStep since dedicated process was created. 53 // - Added functio << 55 // - Added function for switching off scaling 54 // of heavy ions << 56 // of heavy ions from ICRU 73 data 55 // - Minor fix in << 57 // - Minor fix in ComputeLossForStep function 56 // - Minor fix in 58 // - Minor fix in ComputeDEDXPerVolume (AL) 57 // 23. 11. 2009 - Changed energ 59 // 23. 11. 2009 - Changed energy loss limit from 0.15 to 0.01 58 // to improve ac 60 // to improve accuracy for large steps (AL) 59 // 24. 11. 2009 - Bug fix: Rang << 61 // 24. 11. 2009 - Bug fix: Range calculation corrected if same 60 // materials app 62 // materials appears with different cuts in diff. 61 // regions (adde 63 // regions (added UpdateRangeCache function and 62 // modified Buil 64 // modified BuildRangeVector, ComputeLossForStep 63 // functions acc 65 // functions accordingly, added new cache param.) 64 // - Removed GetRa << 66 // - Removed GetRange function (AL) 65 // 04. 11. 2010 - Moved virtual 67 // 04. 11. 2010 - Moved virtual methods to the source (VI) 66 // 68 // 67 // 69 // 68 // Class description: 70 // Class description: 69 // Model for computing the energy loss of i << 71 // Model for computing the energy loss of ions by employing a 70 // parameterisation of dE/dx tables (by def << 72 // parameterisation of dE/dx tables (by default ICRU 73 tables). For 71 // ion-material combinations and/or project << 73 // ion-material combinations and/or projectile energies not covered 72 // by this model, the G4BraggIonModel and G 74 // by this model, the G4BraggIonModel and G4BetheBloch models are 73 // employed. 75 // employed. 74 // 76 // 75 // Comments: 77 // Comments: 76 // 78 // 77 // =========================================== << 79 // =========================================================================== >> 80 >> 81 78 #include "G4IonParametrisedLossModel.hh" 82 #include "G4IonParametrisedLossModel.hh" 79 #include "G4PhysicalConstants.hh" << 83 #include "G4LPhysicsFreeVector.hh" 80 #include "G4SystemOfUnits.hh" << 81 #include "G4PhysicsFreeVector.hh" << 82 #include "G4IonStoppingData.hh" 84 #include "G4IonStoppingData.hh" 83 #include "G4VIonDEDXTable.hh" 85 #include "G4VIonDEDXTable.hh" 84 #include "G4VIonDEDXScalingAlgorithm.hh" 86 #include "G4VIonDEDXScalingAlgorithm.hh" 85 #include "G4IonDEDXScalingICRU73.hh" 87 #include "G4IonDEDXScalingICRU73.hh" 86 #include "G4BraggIonModel.hh" 88 #include "G4BraggIonModel.hh" 87 #include "G4BetheBlochModel.hh" 89 #include "G4BetheBlochModel.hh" 88 #include "G4ProductionCutsTable.hh" 90 #include "G4ProductionCutsTable.hh" 89 #include "G4ParticleChangeForLoss.hh" 91 #include "G4ParticleChangeForLoss.hh" 90 #include "G4LossTableManager.hh" 92 #include "G4LossTableManager.hh" 91 #include "G4EmParameters.hh" << 92 #include "G4GenericIon.hh" 93 #include "G4GenericIon.hh" 93 #include "G4Electron.hh" 94 #include "G4Electron.hh" 94 #include "G4DeltaAngle.hh" << 95 #include "Randomize.hh" 95 #include "Randomize.hh" 96 #include "G4Exp.hh" << 97 96 98 //#define PRINT_TABLE_BUILT 97 //#define PRINT_TABLE_BUILT >> 98 >> 99 99 // ########################################### 100 // ######################################################################### 100 101 101 G4IonParametrisedLossModel::G4IonParametrisedL 102 G4IonParametrisedLossModel::G4IonParametrisedLossModel( 102 const G4ParticleDefinition*, << 103 const G4ParticleDefinition*, 103 const G4String& nam) << 104 const G4String& name) 104 : G4VEmModel(nam), << 105 : G4VEmModel(name), 105 braggIonModel(0), 106 braggIonModel(0), 106 betheBlochModel(0), 107 betheBlochModel(0), 107 nmbBins(90), 108 nmbBins(90), 108 nmbSubBins(100), 109 nmbSubBins(100), 109 particleChangeLoss(0), 110 particleChangeLoss(0), 110 corrFactor(1.0), 111 corrFactor(1.0), 111 energyLossLimit(0.01), 112 energyLossLimit(0.01), 112 cutEnergies(0), << 113 cutEnergies(0) 113 isInitialised(false) << 114 { 114 { 115 genericIon = G4GenericIon::Definition(); 115 genericIon = G4GenericIon::Definition(); 116 genericIonPDGMass = genericIon->GetPDGMass() << 116 genericIonPDGMass = genericIon -> GetPDGMass(); 117 corrections = G4LossTableManager::Instance() << 117 corrections = G4LossTableManager::Instance() -> EmCorrections(); >> 118 >> 119 // The upper limit of the current model is set to 100 TeV >> 120 SetHighEnergyLimit(100.0 * TeV); 118 121 119 // The Bragg ion and Bethe Bloch models are 122 // The Bragg ion and Bethe Bloch models are instantiated 120 braggIonModel = new G4BraggIonModel(); 123 braggIonModel = new G4BraggIonModel(); 121 betheBlochModel = new G4BetheBlochModel(); 124 betheBlochModel = new G4BetheBlochModel(); 122 125 >> 126 // By default ICRU 73 stopping power tables are loaded: >> 127 AddDEDXTable("ICRU73", >> 128 new G4IonStoppingData("ion_stopping_data/icru73"), >> 129 new G4IonDEDXScalingICRU73()); >> 130 123 // The boundaries for the range tables are s 131 // The boundaries for the range tables are set 124 lowerEnergyEdgeIntegr = 0.025 * MeV; 132 lowerEnergyEdgeIntegr = 0.025 * MeV; 125 upperEnergyEdgeIntegr = betheBlochModel -> H 133 upperEnergyEdgeIntegr = betheBlochModel -> HighEnergyLimit(); 126 134 127 // Cache parameters are set 135 // Cache parameters are set 128 cacheParticle = 0; 136 cacheParticle = 0; 129 cacheMass = 0; 137 cacheMass = 0; 130 cacheElecMassRatio = 0; 138 cacheElecMassRatio = 0; 131 cacheChargeSquare = 0; 139 cacheChargeSquare = 0; 132 140 133 // Cache parameters are set 141 // Cache parameters are set 134 rangeCacheParticle = 0; << 142 rangeCacheParticle = 0; 135 rangeCacheMatCutsCouple = 0; << 143 rangeCacheMatCutsCouple = 0; 136 rangeCacheEnergyRange = 0; << 144 rangeCacheEnergyRange = 0; 137 rangeCacheRangeEnergy = 0; 145 rangeCacheRangeEnergy = 0; 138 146 139 // Cache parameters are set 147 // Cache parameters are set 140 dedxCacheParticle = 0; 148 dedxCacheParticle = 0; 141 dedxCacheMaterial = 0; 149 dedxCacheMaterial = 0; 142 dedxCacheEnergyCut = 0; 150 dedxCacheEnergyCut = 0; 143 dedxCacheIter = lossTableList.end(); 151 dedxCacheIter = lossTableList.end(); 144 dedxCacheTransitionEnergy = 0.0; << 152 dedxCacheTransitionEnergy = 0.0; 145 dedxCacheTransitionFactor = 0.0; 153 dedxCacheTransitionFactor = 0.0; 146 dedxCacheGenIonMassRatio = 0.0; 154 dedxCacheGenIonMassRatio = 0.0; 147 << 148 // default generator << 149 SetAngularDistribution(new G4DeltaAngle()); << 150 } 155 } 151 156 152 // ########################################### 157 // ######################################################################### 153 158 154 G4IonParametrisedLossModel::~G4IonParametrised 159 G4IonParametrisedLossModel::~G4IonParametrisedLossModel() { 155 160 >> 161 // Range vs energy table objects are deleted and the container is cleared >> 162 RangeEnergyTable::iterator iterRange = r.begin(); >> 163 RangeEnergyTable::iterator iterRange_end = r.end(); >> 164 >> 165 for(;iterRange != iterRange_end; iterRange++) delete iterRange -> second; >> 166 r.clear(); >> 167 >> 168 // Energy vs range table objects are deleted and the container is cleared >> 169 EnergyRangeTable::iterator iterEnergy = E.begin(); >> 170 EnergyRangeTable::iterator iterEnergy_end = E.end(); >> 171 >> 172 for(;iterEnergy != iterEnergy_end; iterEnergy++) delete iterEnergy -> second; >> 173 E.clear(); >> 174 156 // dE/dx table objects are deleted and the c 175 // dE/dx table objects are deleted and the container is cleared 157 LossTableList::iterator iterTables = lossTab 176 LossTableList::iterator iterTables = lossTableList.begin(); 158 LossTableList::iterator iterTables_end = los 177 LossTableList::iterator iterTables_end = lossTableList.end(); 159 178 160 for(;iterTables != iterTables_end; ++iterTab << 179 for(;iterTables != iterTables_end; iterTables++) delete *iterTables; 161 lossTableList.clear(); 180 lossTableList.clear(); 162 181 163 // range table << 182 // The Bragg ion and Bethe Bloch objects are deleted 164 RangeEnergyTable::iterator itr = r.begin(); << 183 delete betheBlochModel; 165 RangeEnergyTable::iterator itr_end = r.end() << 184 delete braggIonModel; 166 for(;itr != itr_end; ++itr) { delete itr->se << 167 r.clear(); << 168 << 169 // inverse range << 170 EnergyRangeTable::iterator ite = E.begin(); << 171 EnergyRangeTable::iterator ite_end = E.end() << 172 for(;ite != ite_end; ++ite) { delete ite->se << 173 E.clear(); << 174 << 175 } 185 } 176 186 177 // ########################################### 187 // ######################################################################### 178 188 179 G4double G4IonParametrisedLossModel::MinEnergy 189 G4double G4IonParametrisedLossModel::MinEnergyCut( 180 const G 190 const G4ParticleDefinition*, 181 const G 191 const G4MaterialCutsCouple* couple) { 182 192 183 return couple->GetMaterial()->GetIonisation( << 193 return couple -> GetMaterial() -> GetIonisation() -> >> 194 GetMeanExcitationEnergy(); 184 } 195 } 185 196 186 // ########################################### 197 // ######################################################################### 187 198 188 G4double G4IonParametrisedLossModel::MaxSecond 199 G4double G4IonParametrisedLossModel::MaxSecondaryEnergy( 189 const G4ParticleD 200 const G4ParticleDefinition* particle, 190 G4double kineticE 201 G4double kineticEnergy) { 191 202 192 // ############## Maximum energy of secondar 203 // ############## Maximum energy of secondaries ########################## 193 // Function computes maximum energy of secon << 204 // Function computes maximum energy of secondary electrons which are 194 // released by an ion 205 // released by an ion 195 // 206 // 196 // See Geant4 physics reference manual (vers 207 // See Geant4 physics reference manual (version 9.1), section 9.1.1 197 // << 208 // 198 // Ref.: W.M. Yao et al, Jour. of Phys. G 33 209 // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1. 199 // C.Caso et al. (Part. Data Group), E 210 // C.Caso et al. (Part. Data Group), Europ. Phys. Jour. C 3 1 (1998). 200 // B. Rossi, High energy particles, Ne 211 // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952). 201 // 212 // 202 // (Implementation adapted from G4BraggIonMo 213 // (Implementation adapted from G4BraggIonModel) 203 214 204 if(particle != cacheParticle) UpdateCache(pa 215 if(particle != cacheParticle) UpdateCache(particle); 205 216 206 G4double tau = kineticEnergy/cacheMass; 217 G4double tau = kineticEnergy/cacheMass; 207 G4double tmax = 2.0 * electron_mass_c2 * tau 218 G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) / 208 (1. + 2.0 * (tau + 1.) * cac << 219 (1. + 2.0 * (tau + 1.) * cacheElecMassRatio + 209 cacheElecMassRatio * cacheEl 220 cacheElecMassRatio * cacheElecMassRatio); 210 221 211 return tmax; 222 return tmax; 212 } 223 } 213 224 214 // ########################################### 225 // ######################################################################### 215 226 216 G4double G4IonParametrisedLossModel::GetCharge 227 G4double G4IonParametrisedLossModel::GetChargeSquareRatio( 217 const G4ParticleD 228 const G4ParticleDefinition* particle, 218 const G4Material* material, 229 const G4Material* material, 219 G4double kineticE 230 G4double kineticEnergy) { // Kinetic energy 220 231 221 G4double chargeSquareRatio = corrections -> 232 G4double chargeSquareRatio = corrections -> 222 Effective 233 EffectiveChargeSquareRatio(particle, 223 234 material, 224 235 kineticEnergy); 225 corrFactor = chargeSquareRatio * << 236 corrFactor = chargeSquareRatio * 226 corrections -> Effectiv 237 corrections -> EffectiveChargeCorrection(particle, 227 238 material, 228 239 kineticEnergy); 229 return corrFactor; 240 return corrFactor; 230 } 241 } 231 242 232 // ########################################### 243 // ######################################################################### 233 244 234 G4double G4IonParametrisedLossModel::GetPartic 245 G4double G4IonParametrisedLossModel::GetParticleCharge( 235 const G4ParticleD 246 const G4ParticleDefinition* particle, 236 const G4Material* material, 247 const G4Material* material, 237 G4double kineticE << 248 G4double kineticEnergy) { // Kinetic energy 238 249 239 return corrections -> GetParticleCharge(part 250 return corrections -> GetParticleCharge(particle, material, kineticEnergy); 240 } 251 } 241 252 242 // ########################################### 253 // ######################################################################### 243 254 244 void G4IonParametrisedLossModel::Initialise( 255 void G4IonParametrisedLossModel::Initialise( 245 const G 256 const G4ParticleDefinition* particle, 246 const G 257 const G4DataVector& cuts) { 247 258 248 // Cached parameters are reset 259 // Cached parameters are reset 249 cacheParticle = 0; 260 cacheParticle = 0; 250 cacheMass = 0; 261 cacheMass = 0; 251 cacheElecMassRatio = 0; 262 cacheElecMassRatio = 0; 252 cacheChargeSquare = 0; 263 cacheChargeSquare = 0; 253 << 264 254 // Cached parameters are reset 265 // Cached parameters are reset 255 rangeCacheParticle = 0; << 266 rangeCacheParticle = 0; 256 rangeCacheMatCutsCouple = 0; << 267 rangeCacheMatCutsCouple = 0; 257 rangeCacheEnergyRange = 0; << 268 rangeCacheEnergyRange = 0; 258 rangeCacheRangeEnergy = 0; 269 rangeCacheRangeEnergy = 0; 259 270 260 // Cached parameters are reset 271 // Cached parameters are reset 261 dedxCacheParticle = 0; 272 dedxCacheParticle = 0; 262 dedxCacheMaterial = 0; 273 dedxCacheMaterial = 0; 263 dedxCacheEnergyCut = 0; 274 dedxCacheEnergyCut = 0; 264 dedxCacheIter = lossTableList.end(); 275 dedxCacheIter = lossTableList.end(); 265 dedxCacheTransitionEnergy = 0.0; << 276 dedxCacheTransitionEnergy = 0.0; 266 dedxCacheTransitionFactor = 0.0; 277 dedxCacheTransitionFactor = 0.0; 267 dedxCacheGenIonMassRatio = 0.0; 278 dedxCacheGenIonMassRatio = 0.0; 268 279 269 // By default ICRU 73 stopping power tables << 270 if(!isInitialised) { << 271 G4bool icru90 = G4EmParameters::Instance() << 272 isInitialised = true; << 273 AddDEDXTable("ICRU73", << 274 new G4IonStoppingData("ion_stopping_data/ << 275 new G4IonDEDXScalingICRU73()); << 276 } << 277 // The cache of loss tables is cleared 280 // The cache of loss tables is cleared 278 LossTableList::iterator iterTables = lossTab 281 LossTableList::iterator iterTables = lossTableList.begin(); 279 LossTableList::iterator iterTables_end = los 282 LossTableList::iterator iterTables_end = lossTableList.end(); 280 283 281 for(;iterTables != iterTables_end; ++iterTab << 284 for(;iterTables != iterTables_end; iterTables++) 282 (*iterTables) -> ClearCache(); << 285 (*iterTables) -> ClearCache(); 283 } << 284 286 285 // Range vs energy and energy vs range vecto 287 // Range vs energy and energy vs range vectors from previous runs are 286 // cleared 288 // cleared 287 RangeEnergyTable::iterator iterRange = r.beg 289 RangeEnergyTable::iterator iterRange = r.begin(); 288 RangeEnergyTable::iterator iterRange_end = r 290 RangeEnergyTable::iterator iterRange_end = r.end(); 289 291 290 for(;iterRange != iterRange_end; ++iterRange << 292 for(;iterRange != iterRange_end; iterRange++) delete iterRange -> second; 291 delete iterRange->second; << 292 } << 293 r.clear(); 293 r.clear(); 294 294 295 EnergyRangeTable::iterator iterEnergy = E.be 295 EnergyRangeTable::iterator iterEnergy = E.begin(); 296 EnergyRangeTable::iterator iterEnergy_end = 296 EnergyRangeTable::iterator iterEnergy_end = E.end(); 297 297 298 for(;iterEnergy != iterEnergy_end; ++iterEne << 298 for(;iterEnergy != iterEnergy_end; iterEnergy++) delete iterEnergy -> second; 299 delete iterEnergy->second; << 300 } << 301 E.clear(); 299 E.clear(); 302 300 303 // The cut energies << 301 // The cut energies are (re)loaded 304 cutEnergies = cuts; << 302 size_t size = cuts.size(); >> 303 cutEnergies.clear(); >> 304 for(size_t i = 0; i < size; i++) cutEnergies.push_back(cuts[i]); 305 305 306 // All dE/dx vectors are built 306 // All dE/dx vectors are built 307 const G4ProductionCutsTable* coupleTable= 307 const G4ProductionCutsTable* coupleTable= 308 G4ProductionCutsTable::Ge 308 G4ProductionCutsTable::GetProductionCutsTable(); 309 G4int nmbCouples = (G4int)coupleTable->GetTa << 309 size_t nmbCouples = coupleTable -> GetTableSize(); 310 310 311 #ifdef PRINT_TABLE_BUILT 311 #ifdef PRINT_TABLE_BUILT 312 G4cout << "G4IonParametrisedLossModel::Ini 312 G4cout << "G4IonParametrisedLossModel::Initialise():" 313 << " Building dE/dx vectors:" 313 << " Building dE/dx vectors:" 314 << G4endl; << 314 << G4endl; 315 #endif 315 #endif 316 316 317 for (G4int i = 0; i < nmbCouples; ++i) { << 317 for (size_t i = 0; i < nmbCouples; i++) { 318 318 319 const G4MaterialCutsCouple* couple = coupl << 319 const G4MaterialCutsCouple* couple = 320 const G4Material* material = couple->GetMa << 320 coupleTable -> GetMaterialCutsCouple(i); 321 321 322 for(G4int atomicNumberIon = 3; atomicNumbe << 322 const G4Material* material = couple -> GetMaterial(); >> 323 // G4ProductionCuts* productionCuts = couple -> GetProductionCuts(); >> 324 >> 325 for(G4int atomicNumberIon = 3; atomicNumberIon < 102; atomicNumberIon++) { 323 326 324 LossTableList::iterator iter = lossTabl 327 LossTableList::iterator iter = lossTableList.begin(); 325 LossTableList::iterator iter_end = loss 328 LossTableList::iterator iter_end = lossTableList.end(); 326 329 327 for(;iter != iter_end; ++iter) { << 330 for(;iter != iter_end; iter++) { 328 331 329 if(*iter == 0) { 332 if(*iter == 0) { 330 G4cout << "G4IonParametrisedLoss 333 G4cout << "G4IonParametrisedLossModel::Initialise():" 331 << " Skipping illegal tab 334 << " Skipping illegal table." 332 << G4endl; << 335 << G4endl; 333 } 336 } 334 << 337 335 if((*iter)->BuildDEDXTable(atomicNumberIon << 338 G4bool isApplicable = >> 339 (*iter) -> BuildDEDXTable(atomicNumberIon, material); >> 340 if(isApplicable) { 336 341 337 #ifdef PRINT_TABLE_BUILT 342 #ifdef PRINT_TABLE_BUILT 338 G4cout << " Atomic Number Ion = 343 G4cout << " Atomic Number Ion = " << atomicNumberIon 339 << ", Material = " << mate 344 << ", Material = " << material -> GetName() 340 << ", Table = " << (*iter) 345 << ", Table = " << (*iter) -> GetName() 341 << G4endl; << 346 << G4endl; 342 #endif 347 #endif 343 break; << 348 break; 344 } 349 } 345 } 350 } 346 } 351 } 347 } 352 } 348 353 349 // The particle change object << 354 // The particle change object 350 if(! particleChangeLoss) { 355 if(! particleChangeLoss) { 351 particleChangeLoss = GetParticleChangeForL 356 particleChangeLoss = GetParticleChangeForLoss(); 352 braggIonModel->SetParticleChange(particleC << 353 betheBlochModel->SetParticleChange(particl << 354 } 357 } 355 << 358 356 // The G4BraggIonModel and G4BetheBlochModel 359 // The G4BraggIonModel and G4BetheBlochModel instances are initialised with 357 // the same settings as the current model: 360 // the same settings as the current model: 358 braggIonModel -> Initialise(particle, cuts); 361 braggIonModel -> Initialise(particle, cuts); 359 betheBlochModel -> Initialise(particle, cuts 362 betheBlochModel -> Initialise(particle, cuts); 360 } 363 } 361 364 362 // ########################################### 365 // ######################################################################### 363 366 364 G4double G4IonParametrisedLossModel::ComputeCr 367 G4double G4IonParametrisedLossModel::ComputeCrossSectionPerAtom( 365 const G4Par 368 const G4ParticleDefinition* particle, 366 G4double ki 369 G4double kineticEnergy, 367 G4double atomicNumber, << 370 G4double atomicNumber, 368 G4double, 371 G4double, 369 G4double cu 372 G4double cutEnergy, 370 G4double ma 373 G4double maxKinEnergy) { 371 374 372 // ############## Cross section per atom ## 375 // ############## Cross section per atom ################################ 373 // Function computes ionization cross sectio 376 // Function computes ionization cross section per atom 374 // 377 // 375 // See Geant4 physics reference manual (vers 378 // See Geant4 physics reference manual (version 9.1), section 9.1.3 376 // << 379 // 377 // Ref.: W.M. Yao et al, Jour. of Phys. G 33 380 // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1. 378 // B. Rossi, High energy particles, Ne 381 // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952). 379 // 382 // 380 // (Implementation adapted from G4BraggIonMo 383 // (Implementation adapted from G4BraggIonModel) 381 384 382 G4double crosssection = 0.0; 385 G4double crosssection = 0.0; 383 G4double tmax = MaxSecondaryEnergy(part 386 G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy); 384 G4double maxEnergy = std::min(tmax, maxKinEn 387 G4double maxEnergy = std::min(tmax, maxKinEnergy); 385 388 386 if(cutEnergy < tmax) { 389 if(cutEnergy < tmax) { 387 390 388 G4double energy = kineticEnergy + cacheM 391 G4double energy = kineticEnergy + cacheMass; 389 G4double betaSquared = kineticEnergy * << 392 G4double betaSquared = kineticEnergy * 390 (energy + 393 (energy + cacheMass) / (energy * energy); 391 394 392 crosssection = 1.0 / cutEnergy - 1.0 / ma << 395 crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy - 393 betaSquared * std::lo 396 betaSquared * std::log(maxEnergy / cutEnergy) / tmax; 394 397 395 crosssection *= twopi_mc2_rcl2 * cacheCha 398 crosssection *= twopi_mc2_rcl2 * cacheChargeSquare / betaSquared; 396 } 399 } 397 400 398 #ifdef PRINT_DEBUG_CS 401 #ifdef PRINT_DEBUG_CS 399 G4cout << "################################# 402 G4cout << "########################################################" 400 << G4endl 403 << G4endl 401 << "# G4IonParametrisedLossModel::Com 404 << "# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom" 402 << G4endl 405 << G4endl 403 << "# particle =" << particle -> GetP << 406 << "# particle =" << particle -> GetParticleName() 404 << G4endl 407 << G4endl 405 << "# cut(MeV) = " << cutEnergy/MeV << 408 << "# cut(MeV) = " << cutEnergy/MeV 406 << G4endl; 409 << G4endl; 407 410 408 G4cout << "#" 411 G4cout << "#" 409 << std::setw(13) << std::right << "E( 412 << std::setw(13) << std::right << "E(MeV)" 410 << std::setw(14) << "CS(um)" 413 << std::setw(14) << "CS(um)" 411 << std::setw(14) << "E_max_sec(MeV)" 414 << std::setw(14) << "E_max_sec(MeV)" 412 << G4endl 415 << G4endl 413 << "# ------------------------------- 416 << "# ------------------------------------------------------" 414 << G4endl; 417 << G4endl; 415 418 416 G4cout << std::setw(14) << std::right << kin 419 G4cout << std::setw(14) << std::right << kineticEnergy / MeV 417 << std::setw(14) << crosssection / (u 420 << std::setw(14) << crosssection / (um * um) 418 << std::setw(14) << tmax / MeV 421 << std::setw(14) << tmax / MeV 419 << G4endl; 422 << G4endl; 420 #endif 423 #endif 421 424 422 crosssection *= atomicNumber; 425 crosssection *= atomicNumber; 423 426 424 return crosssection; 427 return crosssection; 425 } 428 } 426 429 427 // ########################################### 430 // ######################################################################### 428 431 429 G4double G4IonParametrisedLossModel::CrossSect 432 G4double G4IonParametrisedLossModel::CrossSectionPerVolume( 430 const G4Material* material, 433 const G4Material* material, 431 const G4Par 434 const G4ParticleDefinition* particle, 432 G4double ki 435 G4double kineticEnergy, 433 G4double cu 436 G4double cutEnergy, 434 G4double ma 437 G4double maxEnergy) { 435 438 436 G4double nbElecPerVolume = material -> GetTo 439 G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume(); 437 G4double cross = ComputeCrossSectionPerAtom( << 440 G4double cross = ComputeCrossSectionPerAtom(particle, 438 << 441 kineticEnergy, 439 442 nbElecPerVolume, 0, 440 443 cutEnergy, 441 444 maxEnergy); 442 445 443 return cross; 446 return cross; 444 } 447 } 445 448 446 // ########################################### 449 // ######################################################################### 447 450 448 G4double G4IonParametrisedLossModel::ComputeDE 451 G4double G4IonParametrisedLossModel::ComputeDEDXPerVolume( 449 const G4 452 const G4Material* material, 450 const G4ParticleDefinition 453 const G4ParticleDefinition* particle, 451 G4double kineticEnergy, 454 G4double kineticEnergy, 452 G4double cutEnergy) { 455 G4double cutEnergy) { 453 456 454 // ############## dE/dx #################### 457 // ############## dE/dx ################################################## 455 // Function computes dE/dx values, where fol 458 // Function computes dE/dx values, where following rules are adopted: 456 // A. If the ion-material pair is covered 459 // A. If the ion-material pair is covered by any native ion data 457 // parameterisation, then: << 460 // parameterisation, then: 458 // * This parameterization is used for << 461 // * This parameterization is used for energies below a given energy 459 // limit, 462 // limit, 460 // * whereas above the limit the Bethe- << 463 // * whereas above the limit the Bethe-Bloch model is applied, in 461 // combination with an effective char << 464 // combination with an effective charge estimate and high order 462 // correction terms. 465 // correction terms. 463 // A smoothing procedure is applied to 466 // A smoothing procedure is applied to dE/dx values computed with 464 // the second approach. The smoothing f 467 // the second approach. The smoothing factor is based on the dE/dx 465 // values of both approaches at the tra 468 // values of both approaches at the transition energy (high order 466 // correction terms are included in the 469 // correction terms are included in the calculation of the transition 467 // factor). 470 // factor). 468 // B. If the particle is a generic ion, th 471 // B. If the particle is a generic ion, the BraggIon and Bethe-Bloch 469 // models are used and a smoothing proc 472 // models are used and a smoothing procedure is applied to values 470 // obtained with the second approach. 473 // obtained with the second approach. 471 // C. If the ion-material is not covered b 474 // C. If the ion-material is not covered by any ion data parameterization 472 // then: 475 // then: 473 // * The BraggIon model is used for ene << 476 // * The BraggIon model is used for energies below a given energy 474 // limit, 477 // limit, 475 // * whereas above the limit the Bethe- << 478 // * whereas above the limit the Bethe-Bloch model is applied, in 476 // combination with an effective char << 479 // combination with an effective charge estimate and high order 477 // correction terms. 480 // correction terms. 478 // Also in this case, a smoothing proce << 481 // Also in this case, a smoothing procedure is applied to dE/dx values 479 // computed with the second model. << 482 // computed with the second model. 480 483 481 G4double dEdx = 0.0; 484 G4double dEdx = 0.0; 482 UpdateDEDXCache(particle, material, cutEnerg 485 UpdateDEDXCache(particle, material, cutEnergy); 483 486 484 LossTableList::iterator iter = dedxCacheIter 487 LossTableList::iterator iter = dedxCacheIter; 485 488 486 if(iter != lossTableList.end()) { 489 if(iter != lossTableList.end()) { 487 490 488 G4double transitionEnergy = dedxCacheTran 491 G4double transitionEnergy = dedxCacheTransitionEnergy; 489 492 490 if(transitionEnergy > kineticEnergy) { 493 if(transitionEnergy > kineticEnergy) { 491 494 492 dEdx = (*iter) -> GetDEDX(particle, ma 495 dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy); 493 496 494 G4double dEdxDeltaRays = DeltaRayMeanE << 497 G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material, 495 par << 498 particle, 496 kin << 499 kineticEnergy, 497 cut 500 cutEnergy); 498 dEdx -= dEdxDeltaRays; << 501 dEdx -= dEdxDeltaRays; 499 } 502 } 500 else { 503 else { 501 G4double massRatio = dedxCacheGenIonMa 504 G4double massRatio = dedxCacheGenIonMassRatio; 502 << 505 503 G4double chargeSquare = << 506 G4double chargeSquare = 504 GetChargeSquareRatio(pa 507 GetChargeSquareRatio(particle, material, kineticEnergy); 505 508 506 G4double scaledKineticEnergy = kinetic 509 G4double scaledKineticEnergy = kineticEnergy * massRatio; 507 G4double scaledTransitionEnergy = tran 510 G4double scaledTransitionEnergy = transitionEnergy * massRatio; 508 511 509 G4double lowEnergyLimit = betheBlochMo 512 G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit(); 510 513 511 if(scaledTransitionEnergy >= lowEnergy 514 if(scaledTransitionEnergy >= lowEnergyLimit) { 512 515 513 dEdx = betheBlochModel -> ComputeDE 516 dEdx = betheBlochModel -> ComputeDEDXPerVolume( 514 material 517 material, genericIon, 515 scaledKineticEnergy, cutEner 518 scaledKineticEnergy, cutEnergy); 516 << 519 517 dEdx *= chargeSquare; 520 dEdx *= chargeSquare; 518 521 519 dEdx += corrections -> ComputeIonCo << 522 dEdx += corrections -> ComputeIonCorrections(particle, 520 523 material, kineticEnergy); 521 524 522 G4double factor = 1.0 + dedxCacheTr << 525 G4double factor = 1.0 + dedxCacheTransitionFactor / 523 526 kineticEnergy; 524 527 525 dEdx *= factor; 528 dEdx *= factor; 526 } 529 } 527 } 530 } 528 } 531 } 529 else { 532 else { 530 G4double massRatio = 1.0; 533 G4double massRatio = 1.0; 531 G4double chargeSquare = 1.0; 534 G4double chargeSquare = 1.0; 532 535 533 if(particle != genericIon) { 536 if(particle != genericIon) { 534 537 535 chargeSquare = GetChargeSquareRatio(pa 538 chargeSquare = GetChargeSquareRatio(particle, material, kineticEnergy); 536 massRatio = genericIonPDGMass / partic << 539 massRatio = genericIonPDGMass / particle -> GetPDGMass(); 537 } 540 } 538 541 539 G4double scaledKineticEnergy = kineticEne 542 G4double scaledKineticEnergy = kineticEnergy * massRatio; 540 543 541 G4double lowEnergyLimit = betheBlochModel 544 G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit(); 542 if(scaledKineticEnergy < lowEnergyLimit) 545 if(scaledKineticEnergy < lowEnergyLimit) { 543 dEdx = braggIonModel -> ComputeDEDXPer 546 dEdx = braggIonModel -> ComputeDEDXPerVolume( 544 material 547 material, genericIon, 545 scaledKineticEnergy, cutEner 548 scaledKineticEnergy, cutEnergy); 546 549 547 dEdx *= chargeSquare; 550 dEdx *= chargeSquare; 548 } 551 } 549 else { 552 else { 550 G4double dEdxLimitParam = braggIonMode 553 G4double dEdxLimitParam = braggIonModel -> ComputeDEDXPerVolume( 551 material 554 material, genericIon, 552 lowEnergyLimit, cutEnergy); 555 lowEnergyLimit, cutEnergy); 553 556 554 G4double dEdxLimitBetheBloch = betheBl 557 G4double dEdxLimitBetheBloch = betheBlochModel -> ComputeDEDXPerVolume( 555 material 558 material, genericIon, 556 lowEnergyLimit, cutEnergy); 559 lowEnergyLimit, cutEnergy); 557 560 558 if(particle != genericIon) { 561 if(particle != genericIon) { 559 G4double chargeSquareLowEnergyLimit << 562 G4double chargeSquareLowEnergyLimit = 560 GetChargeSquareRatio(particle, << 563 GetChargeSquareRatio(particle, material, 561 lowEnergyL 564 lowEnergyLimit / massRatio); 562 565 563 dEdxLimitParam *= chargeSquareLowEn 566 dEdxLimitParam *= chargeSquareLowEnergyLimit; 564 dEdxLimitBetheBloch *= chargeSquare 567 dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit; 565 568 566 dEdxLimitBetheBloch += << 569 dEdxLimitBetheBloch += 567 corrections -> ComputeIonC << 570 corrections -> ComputeIonCorrections(particle, 568 material 571 material, lowEnergyLimit / massRatio); 569 } 572 } 570 573 571 G4double factor = (1.0 + (dEdxLimitPar 574 G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0) 572 * lowEnergyLimi 575 * lowEnergyLimit / scaledKineticEnergy); 573 576 574 dEdx = betheBlochModel -> ComputeDEDXP 577 dEdx = betheBlochModel -> ComputeDEDXPerVolume( 575 material 578 material, genericIon, 576 scaledKineticEnergy, cutEner 579 scaledKineticEnergy, cutEnergy); 577 580 578 dEdx *= chargeSquare; 581 dEdx *= chargeSquare; 579 582 580 if(particle != genericIon) { 583 if(particle != genericIon) { 581 dEdx += corrections -> ComputeIonCo << 584 dEdx += corrections -> ComputeIonCorrections(particle, 582 material 585 material, kineticEnergy); 583 } 586 } 584 587 585 dEdx *= factor; 588 dEdx *= factor; 586 } 589 } 587 590 588 } 591 } 589 592 590 if (dEdx < 0.0) dEdx = 0.0; 593 if (dEdx < 0.0) dEdx = 0.0; 591 594 592 return dEdx; 595 return dEdx; 593 } 596 } 594 597 595 // ########################################### 598 // ######################################################################### 596 599 597 void G4IonParametrisedLossModel::PrintDEDXTabl 600 void G4IonParametrisedLossModel::PrintDEDXTable( 598 const G4ParticleDefinition* << 601 const G4ParticleDefinition* particle, // Projectile (ion) 599 const G4Material* material, 602 const G4Material* material, // Absorber material 600 G4double lowerBoundary, 603 G4double lowerBoundary, // Minimum energy per nucleon 601 G4double upperBoundary, 604 G4double upperBoundary, // Maximum energy per nucleon 602 G4int numBins, << 605 G4int nmbBins, // Number of bins 603 G4bool logScaleEnergy) { 606 G4bool logScaleEnergy) { // Logarithmic scaling of energy 604 607 605 G4double atomicMassNumber = particle -> GetA 608 G4double atomicMassNumber = particle -> GetAtomicMass(); 606 G4double materialDensity = material -> GetDe 609 G4double materialDensity = material -> GetDensity(); 607 610 608 G4cout << "# dE/dx table for " << particle - << 611 G4cout << "# dE/dx table for " << particle -> GetParticleName() 609 << " in material " << material -> Get 612 << " in material " << material -> GetName() 610 << " of density " << materialDensity 613 << " of density " << materialDensity / g * cm3 611 << " g/cm3" 614 << " g/cm3" 612 << G4endl 615 << G4endl 613 << "# Projectile mass number A1 = " < 616 << "# Projectile mass number A1 = " << atomicMassNumber 614 << G4endl 617 << G4endl 615 << "# ------------------------------- 618 << "# ------------------------------------------------------" 616 << G4endl; 619 << G4endl; 617 G4cout << "#" 620 G4cout << "#" 618 << std::setw(13) << std::right << "E" 621 << std::setw(13) << std::right << "E" 619 << std::setw(14) << "E/A1" 622 << std::setw(14) << "E/A1" 620 << std::setw(14) << "dE/dx" 623 << std::setw(14) << "dE/dx" 621 << std::setw(14) << "1/rho*dE/dx" 624 << std::setw(14) << "1/rho*dE/dx" 622 << G4endl; 625 << G4endl; 623 G4cout << "#" 626 G4cout << "#" 624 << std::setw(13) << std::right << "(M 627 << std::setw(13) << std::right << "(MeV)" 625 << std::setw(14) << "(MeV)" 628 << std::setw(14) << "(MeV)" 626 << std::setw(14) << "(MeV/cm)" 629 << std::setw(14) << "(MeV/cm)" 627 << std::setw(14) << "(MeV*cm2/mg)" 630 << std::setw(14) << "(MeV*cm2/mg)" 628 << G4endl 631 << G4endl 629 << "# ------------------------------- 632 << "# ------------------------------------------------------" 630 << G4endl; 633 << G4endl; 631 634 632 G4double energyLowerBoundary = lowerBoundary 635 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber; 633 G4double energyUpperBoundary = upperBoundary << 636 G4double energyUpperBoundary = upperBoundary * atomicMassNumber; 634 637 635 if(logScaleEnergy) { 638 if(logScaleEnergy) { 636 639 637 energyLowerBoundary = std::log(energyLowe 640 energyLowerBoundary = std::log(energyLowerBoundary); 638 energyUpperBoundary = std::log(energyUppe << 641 energyUpperBoundary = std::log(energyUpperBoundary); 639 } 642 } 640 643 641 G4double deltaEnergy = (energyUpperBoundary << 644 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) / 642 645 G4double(nmbBins); 643 646 644 for(int i = 0; i < numBins + 1; i++) { << 647 for(int i = 0; i < nmbBins + 1; i++) { 645 648 646 G4double energy = energyLowerBoundary + 649 G4double energy = energyLowerBoundary + i * deltaEnergy; 647 if(logScaleEnergy) energy = G4Exp(energy << 650 if(logScaleEnergy) energy = std::exp(energy); 648 651 649 G4double dedx = ComputeDEDXPerVolume(mat 652 G4double dedx = ComputeDEDXPerVolume(material, particle, energy, DBL_MAX); 650 G4cout.precision(6); 653 G4cout.precision(6); 651 G4cout << std::setw(14) << std::right << 654 G4cout << std::setw(14) << std::right << energy / MeV 652 << std::setw(14) << energy / atom 655 << std::setw(14) << energy / atomicMassNumber / MeV 653 << std::setw(14) << dedx / MeV * 656 << std::setw(14) << dedx / MeV * cm 654 << std::setw(14) << dedx / materi << 657 << std::setw(14) << dedx / materialDensity / (MeV*cm2/(0.001*g)) 655 << G4endl; 658 << G4endl; 656 } 659 } 657 } 660 } 658 661 659 // ########################################### 662 // ######################################################################### 660 663 661 void G4IonParametrisedLossModel::PrintDEDXTabl 664 void G4IonParametrisedLossModel::PrintDEDXTableHandlers( 662 const G4ParticleDefinition* << 665 const G4ParticleDefinition* particle, // Projectile (ion) 663 const G4Material* material, 666 const G4Material* material, // Absorber material 664 G4double lowerBoundary, 667 G4double lowerBoundary, // Minimum energy per nucleon 665 G4double upperBoundary, 668 G4double upperBoundary, // Maximum energy per nucleon 666 G4int numBins, << 669 G4int nmbBins, // Number of bins 667 G4bool logScaleEnergy) { 670 G4bool logScaleEnergy) { // Logarithmic scaling of energy 668 671 669 LossTableList::iterator iter = lossTableList 672 LossTableList::iterator iter = lossTableList.begin(); 670 LossTableList::iterator iter_end = lossTable 673 LossTableList::iterator iter_end = lossTableList.end(); 671 674 672 for(;iter != iter_end; iter++) { << 675 for(;iter != iter_end; iter++) { 673 G4bool isApplicable = (*iter) -> IsAppl 676 G4bool isApplicable = (*iter) -> IsApplicable(particle, material); 674 if(isApplicable) { << 677 if(isApplicable) { 675 (*iter) -> PrintDEDXTable(particle, material 678 (*iter) -> PrintDEDXTable(particle, material, 676 lowerBoundar << 679 lowerBoundary, upperBoundary, 677 numBins,logS << 680 nmbBins,logScaleEnergy); 678 break; 681 break; 679 } << 682 } 680 } 683 } 681 } 684 } 682 685 683 // ########################################### 686 // ######################################################################### 684 687 685 void G4IonParametrisedLossModel::SampleSeconda 688 void G4IonParametrisedLossModel::SampleSecondaries( 686 std::vector<G4Dyn 689 std::vector<G4DynamicParticle*>* secondaries, 687 const G4MaterialCutsCouple* couple, << 690 const G4MaterialCutsCouple*, 688 const G4DynamicParticle* particle, 691 const G4DynamicParticle* particle, 689 G4double cutKinEnergySec, 692 G4double cutKinEnergySec, 690 G4double userMaxKinEnergySec) { 693 G4double userMaxKinEnergySec) { >> 694 >> 695 691 // ############## Sampling of secondaries ## 696 // ############## Sampling of secondaries ################################# 692 // The probability density function (pdf) of << 697 // The probability density function (pdf) of the kinetic energy T of a 693 // secondary electron may be written as: 698 // secondary electron may be written as: 694 // pdf(T) = f(T) * g(T) 699 // pdf(T) = f(T) * g(T) 695 // where << 700 // where 696 // f(T) = (Tmax - Tcut) / (Tmax * Tcut) * 701 // f(T) = (Tmax - Tcut) / (Tmax * Tcut) * (1 / T^2) 697 // g(T) = 1 - beta^2 * T / Tmax 702 // g(T) = 1 - beta^2 * T / Tmax 698 // where Tmax is the maximum kinetic energy 703 // where Tmax is the maximum kinetic energy of the secondary, Tcut 699 // is the lower energy cut and beta is the k << 704 // is the lower energy cut and beta is the kinetic energy of the 700 // projectile. 705 // projectile. 701 // 706 // 702 // Sampling of the kinetic energy of a secon 707 // Sampling of the kinetic energy of a secondary electron: 703 // 1) T0 is sampled from f(T) using the cum 708 // 1) T0 is sampled from f(T) using the cumulated distribution function 704 // F(T) = int_Tcut^T f(T')dT' 709 // F(T) = int_Tcut^T f(T')dT' 705 // 2) T is accepted or rejected by evaluati 710 // 2) T is accepted or rejected by evaluating the rejection function g(T) 706 // at the sampled energy T0 against a ra 711 // at the sampled energy T0 against a randomly sampled value 707 // 712 // 708 // 713 // 709 // See Geant4 physics reference manual (vers 714 // See Geant4 physics reference manual (version 9.1), section 9.1.4 710 // 715 // 711 // 716 // 712 // Reference pdf: W.M. Yao et al, Jour. of P 717 // Reference pdf: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1. 713 // 718 // 714 // (Implementation adapted from G4BraggIonMo 719 // (Implementation adapted from G4BraggIonModel) 715 720 716 G4double rossiMaxKinEnergySec = MaxSecondary 721 G4double rossiMaxKinEnergySec = MaxSecondaryKinEnergy(particle); 717 G4double maxKinEnergySec = << 722 G4double maxKinEnergySec = 718 std::min(rossiMaxKinEn 723 std::min(rossiMaxKinEnergySec, userMaxKinEnergySec); 719 724 720 if(cutKinEnergySec >= maxKinEnergySec) retur 725 if(cutKinEnergySec >= maxKinEnergySec) return; 721 726 722 G4double kineticEnergy = particle -> GetKine 727 G4double kineticEnergy = particle -> GetKineticEnergy(); >> 728 G4ThreeVector direction = particle ->GetMomentumDirection(); 723 729 724 G4double energy = kineticEnergy + cacheMass 730 G4double energy = kineticEnergy + cacheMass; 725 G4double betaSquared = kineticEnergy * (ener << 731 G4double betaSquared = kineticEnergy * 726 / (energy * energy); << 732 (energy + cacheMass) / (energy * energy); 727 733 728 G4double kinEnergySec; 734 G4double kinEnergySec; 729 G4double grej; << 735 G4double g; 730 736 731 do { 737 do { 732 738 733 // Sampling kinetic energy from f(T) (usin 739 // Sampling kinetic energy from f(T) (using F(T)): 734 G4double xi = G4UniformRand(); 740 G4double xi = G4UniformRand(); 735 kinEnergySec = cutKinEnergySec * maxKinEne << 741 kinEnergySec = cutKinEnergySec * maxKinEnergySec / 736 (maxKinEnergySec * (1. 742 (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi); 737 743 738 // Deriving the value of the rejection fun << 744 // Deriving the value of the rejection function at the obtained kinetic 739 // energy: 745 // energy: 740 grej = 1.0 - betaSquared * kinEnergySec / << 746 g = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec; 741 747 742 if(grej > 1.0) { << 748 if(g > 1.0) { 743 G4cout << "G4IonParametrisedLossModel: 749 G4cout << "G4IonParametrisedLossModel::SampleSecondary Warning: " 744 << "Majorant 1.0 < " 750 << "Majorant 1.0 < " 745 << grej << " for e= " << kinEne << 751 << g << " for e= " << kinEnergySec 746 << G4endl; 752 << G4endl; 747 } 753 } 748 754 749 } while( G4UniformRand() >= grej ); << 755 } while( G4UniformRand() >= g ); 750 756 751 const G4Material* mat = couple->GetMaterial << 757 G4double momentumSec = 752 G4int Z = SelectRandomAtomNumber(mat); << 758 std::sqrt(kinEnergySec * (kinEnergySec + 2.0 * electron_mass_c2)); 753 759 754 const G4ParticleDefinition* electron = G4Ele << 760 G4double totMomentum = energy*std::sqrt(betaSquared); >> 761 G4double cost = kinEnergySec * (energy + electron_mass_c2) / >> 762 (momentumSec * totMomentum); >> 763 if(cost > 1.0) cost = 1.0; >> 764 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost)); 755 765 756 G4DynamicParticle* delta = new G4DynamicPart << 766 G4double phi = twopi * G4UniformRand() ; 757 GetAngularDistribution()->SampleDirection( << 758 << 759 << 760 767 761 secondaries->push_back(delta); << 768 G4ThreeVector directionSec(sint*std::cos(phi),sint*std::sin(phi), cost) ; >> 769 directionSec.rotateUz(direction); 762 770 763 // Change kinematics of primary particle << 771 // create G4DynamicParticle object for delta ray 764 G4ThreeVector direction = particle ->GetMome << 772 G4DynamicParticle* delta = new G4DynamicParticle(G4Electron::Definition(), 765 G4double totalMomentum = std::sqrt(kineticEn << 773 directionSec, >> 774 kinEnergySec); 766 775 767 G4ThreeVector finalP = totalMomentum*directi << 776 secondaries -> push_back(delta); 768 finalP = finalP.unit(); << 769 777 >> 778 // Change kinematics of primary particle 770 kineticEnergy -= kinEnergySec; 779 kineticEnergy -= kinEnergySec; >> 780 G4ThreeVector finalP = direction*totMomentum - directionSec*momentumSec; >> 781 finalP = finalP.unit(); 771 782 772 particleChangeLoss->SetProposedKineticEnergy << 783 particleChangeLoss -> SetProposedKineticEnergy(kineticEnergy); 773 particleChangeLoss->SetProposedMomentumDirec << 784 particleChangeLoss -> SetProposedMomentumDirection(finalP); 774 } 785 } 775 786 776 // ########################################### 787 // ######################################################################### 777 788 778 void G4IonParametrisedLossModel::UpdateRangeCa 789 void G4IonParametrisedLossModel::UpdateRangeCache( 779 const G4ParticleDefinitio 790 const G4ParticleDefinition* particle, 780 const G4MaterialCutsCoupl 791 const G4MaterialCutsCouple* matCutsCouple) { 781 792 782 // ############## Caching ################## 793 // ############## Caching ################################################## 783 // If the ion-material-cut combination is co 794 // If the ion-material-cut combination is covered by any native ion data 784 // parameterisation (for low energies), rang << 795 // parameterisation (for low energies), range vectors are computed 785 796 786 if(particle == rangeCacheParticle && << 797 if(particle == rangeCacheParticle && 787 matCutsCouple == rangeCacheMatCutsCouple) 798 matCutsCouple == rangeCacheMatCutsCouple) { 788 } 799 } 789 else{ 800 else{ 790 rangeCacheParticle = particle; 801 rangeCacheParticle = particle; 791 rangeCacheMatCutsCouple = matCutsCouple; 802 rangeCacheMatCutsCouple = matCutsCouple; 792 803 793 const G4Material* material = matCutsCoupl 804 const G4Material* material = matCutsCouple -> GetMaterial(); 794 LossTableList::iterator iter = IsApplicab 805 LossTableList::iterator iter = IsApplicable(particle, material); 795 806 796 // If any table is applicable, the transi 807 // If any table is applicable, the transition factor is computed: 797 if(iter != lossTableList.end()) { 808 if(iter != lossTableList.end()) { 798 809 799 // Build range-energy and energy-range 810 // Build range-energy and energy-range vectors if they don't exist 800 IonMatCouple ionMatCouple = std::make_ 811 IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple); 801 RangeEnergyTable::iterator iterRange = 812 RangeEnergyTable::iterator iterRange = r.find(ionMatCouple); 802 813 803 if(iterRange == r.end()) BuildRangeVec 814 if(iterRange == r.end()) BuildRangeVector(particle, matCutsCouple); 804 815 805 rangeCacheEnergyRange = E[ionMatCouple << 816 rangeCacheEnergyRange = E[ionMatCouple]; 806 rangeCacheRangeEnergy = r[ionMatCouple 817 rangeCacheRangeEnergy = r[ionMatCouple]; 807 } 818 } 808 else { 819 else { 809 rangeCacheEnergyRange = 0; << 820 rangeCacheEnergyRange = 0; 810 rangeCacheRangeEnergy = 0; 821 rangeCacheRangeEnergy = 0; 811 } 822 } 812 } 823 } 813 } 824 } 814 825 815 // ########################################### 826 // ######################################################################### 816 827 817 void G4IonParametrisedLossModel::UpdateDEDXCac 828 void G4IonParametrisedLossModel::UpdateDEDXCache( 818 const G4ParticleDefinitio 829 const G4ParticleDefinition* particle, 819 const G4Material* materia 830 const G4Material* material, 820 G4double cutEnergy) { 831 G4double cutEnergy) { 821 832 822 // ############## Caching ################## 833 // ############## Caching ################################################## 823 // If the ion-material combination is covere 834 // If the ion-material combination is covered by any native ion data 824 // parameterisation (for low energies), a tr 835 // parameterisation (for low energies), a transition factor is computed 825 // which is applied to Bethe-Bloch results a << 836 // which is applied to Bethe-Bloch results at higher energies to 826 // guarantee a smooth transition. 837 // guarantee a smooth transition. 827 // This factor only needs to be calculated f 838 // This factor only needs to be calculated for the first step an ion 828 // performs inside a certain material. 839 // performs inside a certain material. 829 840 830 if(particle == dedxCacheParticle && << 841 if(particle == dedxCacheParticle && 831 material == dedxCacheMaterial && 842 material == dedxCacheMaterial && 832 cutEnergy == dedxCacheEnergyCut) { 843 cutEnergy == dedxCacheEnergyCut) { 833 } 844 } 834 else { 845 else { 835 846 836 dedxCacheParticle = particle; 847 dedxCacheParticle = particle; 837 dedxCacheMaterial = material; 848 dedxCacheMaterial = material; 838 dedxCacheEnergyCut = cutEnergy; 849 dedxCacheEnergyCut = cutEnergy; 839 850 840 G4double massRatio = genericIonPDGMass / 851 G4double massRatio = genericIonPDGMass / particle -> GetPDGMass(); 841 dedxCacheGenIonMassRatio = massRatio; 852 dedxCacheGenIonMassRatio = massRatio; 842 853 843 LossTableList::iterator iter = IsApplicab 854 LossTableList::iterator iter = IsApplicable(particle, material); 844 dedxCacheIter = iter; 855 dedxCacheIter = iter; 845 856 846 // If any table is applicable, the transi 857 // If any table is applicable, the transition factor is computed: 847 if(iter != lossTableList.end()) { 858 if(iter != lossTableList.end()) { 848 859 849 // Retrieving the transition energy fr 860 // Retrieving the transition energy from the parameterisation table 850 G4double transitionEnergy = << 861 G4double transitionEnergy = 851 (*iter) -> GetUpperEnergyEdge << 862 (*iter) -> GetUpperEnergyEdge(particle, material); 852 dedxCacheTransitionEnergy = transition << 863 dedxCacheTransitionEnergy = transitionEnergy; 853 864 854 // Computing dE/dx from low-energy par << 865 // Computing dE/dx from low-energy parameterisation at 855 // transition energy 866 // transition energy 856 G4double dEdxParam = (*iter) -> GetDED << 867 G4double dEdxParam = (*iter) -> GetDEDX(particle, material, 857 tra 868 transitionEnergy); 858 << 869 859 G4double dEdxDeltaRays = DeltaRayMeanE << 870 G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material, 860 par << 871 particle, 861 tra << 872 transitionEnergy, 862 cut 873 cutEnergy); 863 dEdxParam -= dEdxDeltaRays; << 874 dEdxParam -= dEdxDeltaRays; 864 875 865 // Computing dE/dx from Bethe-Bloch fo << 876 // Computing dE/dx from Bethe-Bloch formula at transition 866 // energy 877 // energy 867 G4double transitionChargeSquare = << 878 G4double transitionChargeSquare = 868 GetChargeSquareRatio(particle, m 879 GetChargeSquareRatio(particle, material, transitionEnergy); 869 880 870 G4double scaledTransitionEnergy = tran 881 G4double scaledTransitionEnergy = transitionEnergy * massRatio; 871 882 872 G4double dEdxBetheBloch = << 883 G4double dEdxBetheBloch = 873 betheBlochModel -> 884 betheBlochModel -> ComputeDEDXPerVolume( 874 materi 885 material, genericIon, 875 scaledTransitionEnergy, cu 886 scaledTransitionEnergy, cutEnergy); 876 dEdxBetheBloch *= transitionChargeSqua 887 dEdxBetheBloch *= transitionChargeSquare; 877 888 878 // Additionally, high order correction 889 // Additionally, high order corrections are added 879 dEdxBetheBloch += << 890 dEdxBetheBloch += 880 corrections -> ComputeIonCorrectio << 891 corrections -> ComputeIonCorrections(particle, 881 892 material, transitionEnergy); 882 893 883 // Computing transition factor from bo 894 // Computing transition factor from both dE/dx values 884 dedxCacheTransitionFactor = << 895 dedxCacheTransitionFactor = 885 (dEdxParam - dEdxBetheBlo 896 (dEdxParam - dEdxBetheBloch)/dEdxBetheBloch 886 * transitionEnerg << 897 * transitionEnergy; 887 } 898 } 888 else { 899 else { 889 << 900 890 dedxCacheParticle = particle; 901 dedxCacheParticle = particle; 891 dedxCacheMaterial = material; 902 dedxCacheMaterial = material; 892 dedxCacheEnergyCut = cutEnergy; 903 dedxCacheEnergyCut = cutEnergy; 893 904 894 dedxCacheGenIonMassRatio = << 905 dedxCacheGenIonMassRatio = 895 genericIonPDGMass 906 genericIonPDGMass / particle -> GetPDGMass(); 896 907 897 dedxCacheTransitionEnergy = 0.0; 908 dedxCacheTransitionEnergy = 0.0; 898 dedxCacheTransitionFactor = 0.0; 909 dedxCacheTransitionFactor = 0.0; 899 } 910 } 900 } 911 } 901 } 912 } 902 913 903 // ########################################### 914 // ######################################################################### 904 915 905 void G4IonParametrisedLossModel::CorrectionsAl 916 void G4IonParametrisedLossModel::CorrectionsAlongStep( 906 const G4MaterialC 917 const G4MaterialCutsCouple* couple, 907 const G4DynamicParticle* dynamicPar 918 const G4DynamicParticle* dynamicParticle, 908 const G4double& l << 919 G4double& eloss, 909 G4double& eloss) { << 920 G4double&, >> 921 G4double length) { 910 922 911 // ############## Corrections for along step 923 // ############## Corrections for along step energy loss calculation ###### 912 // The computed energy loss (due to electron << 924 // The computed energy loss (due to electronic stopping) is overwritten 913 // by this function if an ion data parameter << 925 // by this function if an ion data parameterization is available for the 914 // current ion-material pair. 926 // current ion-material pair. 915 // No action on the energy loss (due to elec << 927 // No action on the energy loss (due to electronic stopping) is performed 916 // if no parameterization is available. In t << 928 // if no parameterization is available. In this case the original 917 // generic ion tables (in combination with t 929 // generic ion tables (in combination with the effective charge) are used 918 // in the along step DoIt function. 930 // in the along step DoIt function. 919 // 931 // 920 // 932 // 921 // (Implementation partly adapted from G4Bra 933 // (Implementation partly adapted from G4BraggIonModel/G4BetheBlochModel) 922 934 923 const G4ParticleDefinition* particle = dynam 935 const G4ParticleDefinition* particle = dynamicParticle -> GetDefinition(); 924 const G4Material* material = couple -> GetMa 936 const G4Material* material = couple -> GetMaterial(); 925 937 926 G4double kineticEnergy = dynamicParticle -> 938 G4double kineticEnergy = dynamicParticle -> GetKineticEnergy(); 927 939 928 if(kineticEnergy == eloss) { return; } 940 if(kineticEnergy == eloss) { return; } 929 941 930 G4double cutEnergy = DBL_MAX; 942 G4double cutEnergy = DBL_MAX; 931 std::size_t cutIndex = couple -> GetIndex(); << 943 size_t cutIndex = couple -> GetIndex(); 932 cutEnergy = cutEnergies[cutIndex]; 944 cutEnergy = cutEnergies[cutIndex]; 933 945 934 UpdateDEDXCache(particle, material, cutEnerg 946 UpdateDEDXCache(particle, material, cutEnergy); 935 947 936 LossTableList::iterator iter = dedxCacheIter 948 LossTableList::iterator iter = dedxCacheIter; 937 949 938 // If parameterization for ions is available 950 // If parameterization for ions is available the electronic energy loss 939 // is overwritten 951 // is overwritten 940 if(iter != lossTableList.end()) { 952 if(iter != lossTableList.end()) { 941 953 942 // The energy loss is calculated using th 954 // The energy loss is calculated using the ComputeDEDXPerVolume function 943 // and the step length (it is assumed tha << 955 // and the step length (it is assumed that dE/dx does not change 944 // considerably along the step) 956 // considerably along the step) 945 eloss = << 957 eloss = 946 length * ComputeDEDXPerVolume(material << 958 length * ComputeDEDXPerVolume(material, particle, 947 kineticE 959 kineticEnergy, cutEnergy); 948 960 949 #ifdef PRINT_DEBUG 961 #ifdef PRINT_DEBUG 950 G4cout.precision(6); << 962 G4cout.precision(6); 951 G4cout << "################################# 963 G4cout << "########################################################" 952 << G4endl 964 << G4endl 953 << "# G4IonParametrisedLossModel::Cor 965 << "# G4IonParametrisedLossModel::CorrectionsAlongStep" 954 << G4endl 966 << G4endl 955 << "# cut(MeV) = " << cutEnergy/MeV << 967 << "# cut(MeV) = " << cutEnergy/MeV 956 << G4endl; 968 << G4endl; 957 969 958 G4cout << "#" << 970 G4cout << "#" 959 << std::setw(13) << std::right << "E( 971 << std::setw(13) << std::right << "E(MeV)" 960 << std::setw(14) << "l(um)" 972 << std::setw(14) << "l(um)" 961 << std::setw(14) << "l*dE/dx(MeV)" 973 << std::setw(14) << "l*dE/dx(MeV)" 962 << std::setw(14) << "(l*dE/dx)/E" 974 << std::setw(14) << "(l*dE/dx)/E" 963 << G4endl 975 << G4endl 964 << "# ------------------------------- 976 << "# ------------------------------------------------------" 965 << G4endl; 977 << G4endl; 966 978 967 G4cout << std::setw(14) << std::right << kin 979 G4cout << std::setw(14) << std::right << kineticEnergy / MeV 968 << std::setw(14) << length / um 980 << std::setw(14) << length / um 969 << std::setw(14) << eloss / MeV 981 << std::setw(14) << eloss / MeV 970 << std::setw(14) << eloss / kineticEn 982 << std::setw(14) << eloss / kineticEnergy * 100.0 971 << G4endl; 983 << G4endl; 972 #endif 984 #endif 973 985 974 // If the energy loss exceeds a certain f 986 // If the energy loss exceeds a certain fraction of the kinetic energy 975 // (the fraction is indicated by the para 987 // (the fraction is indicated by the parameter "energyLossLimit") then 976 // the range tables are used to derive a 988 // the range tables are used to derive a more accurate value of the 977 // energy loss 989 // energy loss 978 if(eloss > energyLossLimit * kineticEnerg 990 if(eloss > energyLossLimit * kineticEnergy) { 979 991 980 eloss = ComputeLossForStep(couple, par << 992 eloss = ComputeLossForStep(couple, particle, 981 kineticEner 993 kineticEnergy,length); 982 994 983 #ifdef PRINT_DEBUG 995 #ifdef PRINT_DEBUG 984 G4cout << "# Correction applied:" 996 G4cout << "# Correction applied:" 985 << G4endl; 997 << G4endl; 986 998 987 G4cout << std::setw(14) << std::right << kin 999 G4cout << std::setw(14) << std::right << kineticEnergy / MeV 988 << std::setw(14) << length / um 1000 << std::setw(14) << length / um 989 << std::setw(14) << eloss / MeV 1001 << std::setw(14) << eloss / MeV 990 << std::setw(14) << eloss / kineticEn 1002 << std::setw(14) << eloss / kineticEnergy * 100.0 991 << G4endl; 1003 << G4endl; 992 #endif 1004 #endif 993 1005 994 } 1006 } 995 } 1007 } 996 1008 997 // For all corrections below a kinetic energ << 1009 // For all corrections below a kinetic energy between the Pre- and 998 // Post-step energy values is used 1010 // Post-step energy values is used 999 G4double energy = kineticEnergy - eloss * 0. 1011 G4double energy = kineticEnergy - eloss * 0.5; 1000 if(energy < 0.0) energy = kineticEnergy * 0 1012 if(energy < 0.0) energy = kineticEnergy * 0.5; 1001 1013 1002 G4double chargeSquareRatio = corrections -> 1014 G4double chargeSquareRatio = corrections -> 1003 Effectiv 1015 EffectiveChargeSquareRatio(particle, 1004 1016 material, 1005 1017 energy); 1006 GetModelOfFluctuations() -> SetParticleAndC << 1018 GetModelOfFluctuations() -> SetParticleAndCharge(particle, 1007 1019 chargeSquareRatio); 1008 1020 1009 // A correction is applied considering the << 1021 // A correction is applied considering the change of the effective charge 1010 // along the step (the parameter "corrFacto << 1022 // along the step (the parameter "corrFactor" refers to the effective 1011 // charge at the beginning of the step). No 1023 // charge at the beginning of the step). Note: the correction is not 1012 // applied for energy loss values deriving << 1024 // applied for energy loss values deriving directly from parameterized 1013 // ion stopping power tables 1025 // ion stopping power tables 1014 G4double transitionEnergy = dedxCacheTransi 1026 G4double transitionEnergy = dedxCacheTransitionEnergy; 1015 1027 1016 if(iter != lossTableList.end() && transitio 1028 if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) { 1017 chargeSquareRatio *= corrections -> Effe 1029 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle, 1018 1030 material, 1019 1031 energy); 1020 1032 1021 G4double chargeSquareRatioCorr = chargeS 1033 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor; 1022 eloss *= chargeSquareRatioCorr; 1034 eloss *= chargeSquareRatioCorr; 1023 } 1035 } 1024 else if (iter == lossTableList.end()) { 1036 else if (iter == lossTableList.end()) { 1025 1037 1026 chargeSquareRatio *= corrections -> Effe 1038 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle, 1027 1039 material, 1028 1040 energy); 1029 1041 1030 G4double chargeSquareRatioCorr = chargeS 1042 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor; 1031 eloss *= chargeSquareRatioCorr; 1043 eloss *= chargeSquareRatioCorr; 1032 } 1044 } 1033 1045 1034 // Ion high order corrections are applied i 1046 // Ion high order corrections are applied if the current model does not 1035 // overwrite the energy loss (i.e. when the 1047 // overwrite the energy loss (i.e. when the effective charge approach is 1036 // used) 1048 // used) 1037 if(iter == lossTableList.end()) { 1049 if(iter == lossTableList.end()) { 1038 1050 1039 G4double scaledKineticEnergy = kineticEn 1051 G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio; 1040 G4double lowEnergyLimit = betheBlochMode 1052 G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit(); 1041 1053 1042 // Corrections are only applied in the B 1054 // Corrections are only applied in the Bethe-Bloch energy region 1043 if(scaledKineticEnergy > lowEnergyLimit) << 1055 if(scaledKineticEnergy > lowEnergyLimit) 1044 eloss += length * << 1056 eloss += length * 1045 corrections -> IonHighOrderCorrec 1057 corrections -> IonHighOrderCorrections(particle, couple, energy); 1046 } 1058 } 1047 } 1059 } 1048 1060 1049 // ########################################## 1061 // ######################################################################### 1050 1062 1051 void G4IonParametrisedLossModel::BuildRangeVe 1063 void G4IonParametrisedLossModel::BuildRangeVector( 1052 const G4ParticleDefiniti 1064 const G4ParticleDefinition* particle, 1053 const G4MaterialCutsCoup 1065 const G4MaterialCutsCouple* matCutsCouple) { 1054 1066 1055 G4double cutEnergy = DBL_MAX; 1067 G4double cutEnergy = DBL_MAX; 1056 std::size_t cutIndex = matCutsCouple -> Get << 1068 size_t cutIndex = matCutsCouple -> GetIndex(); 1057 cutEnergy = cutEnergies[cutIndex]; 1069 cutEnergy = cutEnergies[cutIndex]; 1058 1070 1059 const G4Material* material = matCutsCouple 1071 const G4Material* material = matCutsCouple -> GetMaterial(); 1060 1072 1061 G4double massRatio = genericIonPDGMass / pa 1073 G4double massRatio = genericIonPDGMass / particle -> GetPDGMass(); 1062 1074 1063 G4double lowerEnergy = lowerEnergyEdgeInteg 1075 G4double lowerEnergy = lowerEnergyEdgeIntegr / massRatio; 1064 G4double upperEnergy = upperEnergyEdgeInteg 1076 G4double upperEnergy = upperEnergyEdgeIntegr / massRatio; 1065 1077 1066 G4double logLowerEnergyEdge = std::log(lowe 1078 G4double logLowerEnergyEdge = std::log(lowerEnergy); 1067 G4double logUpperEnergyEdge = std::log(uppe 1079 G4double logUpperEnergyEdge = std::log(upperEnergy); 1068 1080 1069 G4double logDeltaEnergy = (logUpperEnergyEd << 1081 G4double logDeltaEnergy = (logUpperEnergyEdge - logLowerEnergyEdge) / 1070 1082 G4double(nmbBins); 1071 << 1083 1072 G4double logDeltaIntegr = logDeltaEnergy / 1084 G4double logDeltaIntegr = logDeltaEnergy / G4double(nmbSubBins); 1073 1085 1074 G4PhysicsFreeVector* energyRangeVector = << 1086 G4LPhysicsFreeVector* energyRangeVector = 1075 new G4PhysicsFreeVector(nmbBins+1, << 1087 new G4LPhysicsFreeVector(nmbBins+1, 1076 lowerEnergy, << 1088 lowerEnergy, 1077 upperEnergy, << 1089 upperEnergy); 1078 /*spline=*/true); << 1090 energyRangeVector -> SetSpline(true); 1079 1091 1080 G4double dedxLow = ComputeDEDXPerVolume(mat << 1092 G4double dedxLow = ComputeDEDXPerVolume(material, 1081 par << 1093 particle, 1082 low 1094 lowerEnergy, 1083 cut 1095 cutEnergy); 1084 1096 1085 G4double range = 2.0 * lowerEnergy / dedxLo 1097 G4double range = 2.0 * lowerEnergy / dedxLow; 1086 << 1098 1087 energyRangeVector -> PutValues(0, lowerEner 1099 energyRangeVector -> PutValues(0, lowerEnergy, range); 1088 1100 1089 G4double logEnergy = std::log(lowerEnergy); 1101 G4double logEnergy = std::log(lowerEnergy); 1090 for(std::size_t i = 1; i < nmbBins+1; ++i) << 1102 for(size_t i = 1; i < nmbBins+1; i++) { 1091 1103 1092 G4double logEnergyIntegr = logEnergy; 1104 G4double logEnergyIntegr = logEnergy; 1093 1105 1094 for(std::size_t j = 0; j < nmbSubBins; << 1106 for(size_t j = 0; j < nmbSubBins; j++) { 1095 1107 1096 G4double binLowerBoundary = G4Exp(l << 1108 G4double binLowerBoundary = std::exp(logEnergyIntegr); 1097 logEnergyIntegr += logDeltaIntegr; 1109 logEnergyIntegr += logDeltaIntegr; 1098 << 1110 1099 G4double binUpperBoundary = G4Exp(l << 1111 G4double binUpperBoundary = std::exp(logEnergyIntegr); 1100 G4double deltaIntegr = binUpperBoun 1112 G4double deltaIntegr = binUpperBoundary - binLowerBoundary; 1101 << 1113 1102 G4double energyIntegr = binLowerBou 1114 G4double energyIntegr = binLowerBoundary + 0.5 * deltaIntegr; 1103 << 1115 1104 G4double dedxValue = ComputeDEDXPer << 1116 G4double dedxValue = ComputeDEDXPerVolume(material, 1105 << 1117 particle, 1106 1118 energyIntegr, 1107 cutEnergy); 1119 cutEnergy); 1108 1120 1109 if(dedxValue > 0.0) range += deltaI << 1121 if(dedxValue > 0.0) range += deltaIntegr / dedxValue; 1110 1122 1111 #ifdef PRINT_DEBUG_DETAILS 1123 #ifdef PRINT_DEBUG_DETAILS 1112 G4cout << " E = "<< energyInt << 1124 G4cout << " E = "<< energyIntegr/MeV 1113 << " MeV -> dE = " << de 1125 << " MeV -> dE = " << deltaIntegr/MeV 1114 << " MeV -> dE/dx = " << << 1126 << " MeV -> dE/dx = " << dedxValue/MeV*mm 1115 << " MeV/mm -> dE/(dE/dx << 1127 << " MeV/mm -> dE/(dE/dx) = " << deltaIntegr / 1116 1128 dedxValue / mm 1117 << " mm -> range = " << 1129 << " mm -> range = " << range / mm 1118 << " mm " << G4endl; 1130 << " mm " << G4endl; 1119 #endif 1131 #endif 1120 } 1132 } 1121 << 1133 1122 logEnergy += logDeltaEnergy; 1134 logEnergy += logDeltaEnergy; 1123 << 1135 1124 G4double energy = G4Exp(logEnergy); << 1136 G4double energy = std::exp(logEnergy); 1125 1137 1126 energyRangeVector -> PutValues(i, energ 1138 energyRangeVector -> PutValues(i, energy, range); 1127 1139 1128 #ifdef PRINT_DEBUG_DETAILS 1140 #ifdef PRINT_DEBUG_DETAILS 1129 G4cout << "G4IonParametrisedLossModel:: << 1141 G4cout << "G4IonParametrisedLossModel::BuildRangeVector() bin = " 1130 << i <<", E = " 1142 << i <<", E = " 1131 << energy / MeV << " MeV, R = " << 1143 << energy / MeV << " MeV, R = " 1132 << range / mm << " mm" 1144 << range / mm << " mm" 1133 << G4endl; << 1145 << G4endl; 1134 #endif << 1146 #endif 1135 1147 1136 } 1148 } 1137 //vector is filled: activate spline << 1138 energyRangeVector -> FillSecondDerivatives( << 1139 1149 1140 G4double lowerRangeEdge = << 1150 G4bool b; 1141 energyRangeVector -> Valu << 1151 1142 G4double upperRangeEdge = << 1152 G4double lowerRangeEdge = 1143 energyRangeVector -> Valu << 1153 energyRangeVector -> GetValue(lowerEnergy, b); 1144 << 1154 G4double upperRangeEdge = 1145 G4PhysicsFreeVector* rangeEnergyVector << 1155 energyRangeVector -> GetValue(upperEnergy, b); 1146 = new G4PhysicsFreeVector(nmbBins+1, << 1147 lowerRangeEdge, << 1148 upperRangeEdge, << 1149 /*spline=*/true); << 1150 << 1151 for(std::size_t i = 0; i < nmbBins+1; ++i) << 1152 G4double energy = energyRangeVector -> << 1153 rangeEnergyVector -> << 1154 PutValues(i, energyRangeVector - << 1155 } << 1156 1156 1157 rangeEnergyVector -> FillSecondDerivatives( << 1157 G4LPhysicsFreeVector* rangeEnergyVector >> 1158 = new G4LPhysicsFreeVector(nmbBins+1, >> 1159 lowerRangeEdge, >> 1160 upperRangeEdge); >> 1161 rangeEnergyVector -> SetSpline(true); >> 1162 >> 1163 for(size_t i = 0; i < nmbBins+1; i++) { >> 1164 G4double energy = energyRangeVector -> GetLowEdgeEnergy(i); >> 1165 rangeEnergyVector -> >> 1166 PutValues(i, energyRangeVector -> GetValue(energy, b), energy); >> 1167 } 1158 1168 1159 #ifdef PRINT_DEBUG_TABLES 1169 #ifdef PRINT_DEBUG_TABLES 1160 G4cout << *energyLossVector << 1170 G4cout << *energyLossVector 1161 << *energyRangeVector << 1171 << *energyRangeVector 1162 << *rangeEnergyVector << G4endl; << 1172 << *rangeEnergyVector << G4endl; 1163 #endif << 1173 #endif 1164 1174 1165 IonMatCouple ionMatCouple = std::make_pair( << 1175 IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple); 1166 1176 1167 E[ionMatCouple] = energyRangeVector; 1177 E[ionMatCouple] = energyRangeVector; 1168 r[ionMatCouple] = rangeEnergyVector; 1178 r[ionMatCouple] = rangeEnergyVector; 1169 } 1179 } 1170 1180 1171 // ########################################## 1181 // ######################################################################### 1172 1182 1173 G4double G4IonParametrisedLossModel::ComputeL 1183 G4double G4IonParametrisedLossModel::ComputeLossForStep( 1174 const G4MaterialCutsCoup 1184 const G4MaterialCutsCouple* matCutsCouple, 1175 const G4ParticleDefiniti 1185 const G4ParticleDefinition* particle, 1176 G4double kineticEnergy, 1186 G4double kineticEnergy, 1177 G4double stepLength) { 1187 G4double stepLength) { 1178 1188 1179 G4double loss = 0.0; 1189 G4double loss = 0.0; 1180 1190 1181 UpdateRangeCache(particle, matCutsCouple); 1191 UpdateRangeCache(particle, matCutsCouple); 1182 1192 1183 G4PhysicsVector* energyRange = rangeCacheEn 1193 G4PhysicsVector* energyRange = rangeCacheEnergyRange; 1184 G4PhysicsVector* rangeEnergy = rangeCacheRa 1194 G4PhysicsVector* rangeEnergy = rangeCacheRangeEnergy; 1185 1195 1186 if(energyRange != 0 && rangeEnergy != 0) { 1196 if(energyRange != 0 && rangeEnergy != 0) { >> 1197 G4bool b; 1187 1198 1188 G4double lowerEnEdge = energyRange -> En << 1199 G4double lowerEnEdge = energyRange -> GetLowEdgeEnergy( 0 ); 1189 G4double lowerRangeEdge = rangeEnergy -> << 1200 G4double lowerRangeEdge = rangeEnergy -> GetLowEdgeEnergy( 0 ); 1190 1201 1191 // Computing range for pre-step kinetic 1202 // Computing range for pre-step kinetic energy: 1192 G4double range = energyRange -> Value(ki << 1203 G4double range = energyRange -> GetValue(kineticEnergy, b); 1193 1204 1194 // Energy below vector boundary: 1205 // Energy below vector boundary: 1195 if(kineticEnergy < lowerEnEdge) { 1206 if(kineticEnergy < lowerEnEdge) { 1196 1207 1197 range = energyRange -> Value(lowerEn << 1208 range = energyRange -> GetValue(lowerEnEdge, b); 1198 range *= std::sqrt(kineticEnergy / lo 1209 range *= std::sqrt(kineticEnergy / lowerEnEdge); 1199 } 1210 } 1200 1211 1201 #ifdef PRINT_DEBUG 1212 #ifdef PRINT_DEBUG 1202 G4cout << "G4IonParametrisedLossModel::C << 1213 G4cout << "G4IonParametrisedLossModel::ComputeLossForStep() range = " 1203 << range / mm << " mm, step = " < << 1214 << range / mm << " mm, step = " << stepLength / mm << " mm" 1204 << G4endl; 1215 << G4endl; 1205 #endif 1216 #endif 1206 1217 1207 // Remaining range: 1218 // Remaining range: 1208 G4double remRange = range - stepLength; 1219 G4double remRange = range - stepLength; 1209 1220 1210 // If range is smaller than step length, << 1221 // If range is smaller than step length, the loss is set to kinetic 1211 // energy 1222 // energy 1212 if(remRange < 0.0) loss = kineticEnergy; 1223 if(remRange < 0.0) loss = kineticEnergy; 1213 else if(remRange < lowerRangeEdge) { 1224 else if(remRange < lowerRangeEdge) { 1214 << 1225 1215 G4double ratio = remRange / lowerRang 1226 G4double ratio = remRange / lowerRangeEdge; 1216 loss = kineticEnergy - ratio * ratio 1227 loss = kineticEnergy - ratio * ratio * lowerEnEdge; 1217 } 1228 } 1218 else { 1229 else { 1219 1230 1220 G4double energy = rangeEnergy -> Valu << 1231 G4double energy = rangeEnergy -> GetValue(range - stepLength, b); 1221 loss = kineticEnergy - energy; << 1232 loss = kineticEnergy - energy; 1222 } 1233 } 1223 } 1234 } 1224 1235 1225 if(loss < 0.0) loss = 0.0; 1236 if(loss < 0.0) loss = 0.0; 1226 1237 1227 return loss; 1238 return loss; 1228 } 1239 } 1229 1240 1230 // ########################################## 1241 // ######################################################################### 1231 1242 1232 G4bool G4IonParametrisedLossModel::AddDEDXTab 1243 G4bool G4IonParametrisedLossModel::AddDEDXTable( 1233 const G4Strin << 1244 const G4String& name, 1234 G4VIonDEDXTab << 1245 G4VIonDEDXTable* table, 1235 G4VIonDEDXSca 1246 G4VIonDEDXScalingAlgorithm* algorithm) { 1236 1247 1237 if(table == 0) { 1248 if(table == 0) { 1238 G4cout << "G4IonParametrisedLossModel::A << 1249 G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot " 1239 << " add table: Invalid pointer." 1250 << " add table: Invalid pointer." 1240 << G4endl; << 1251 << G4endl; 1241 1252 1242 return false; 1253 return false; 1243 } 1254 } 1244 1255 1245 // Checking uniqueness of name 1256 // Checking uniqueness of name 1246 LossTableList::iterator iter = lossTableLis 1257 LossTableList::iterator iter = lossTableList.begin(); 1247 LossTableList::iterator iter_end = lossTabl 1258 LossTableList::iterator iter_end = lossTableList.end(); >> 1259 >> 1260 for(;iter != iter_end; iter++) { >> 1261 G4String tableName = (*iter) -> GetName(); 1248 1262 1249 for(;iter != iter_end; ++iter) { << 1263 if(tableName == name) { 1250 const G4String& tableName = (*iter)->Get << 1264 G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot " 1251 << 1265 << " add table: Name already exists." 1252 if(tableName == nam) { << 1253 G4cout << "G4IonParametrisedLossModel << 1254 << " add table: Name already e << 1255 << G4endl; 1266 << G4endl; 1256 1267 1257 return false; 1268 return false; 1258 } 1269 } 1259 } 1270 } 1260 1271 1261 G4VIonDEDXScalingAlgorithm* scalingAlgorith 1272 G4VIonDEDXScalingAlgorithm* scalingAlgorithm = algorithm; 1262 if(scalingAlgorithm == 0) << 1273 if(scalingAlgorithm == 0) 1263 scalingAlgorithm = new G4VIonDEDXScaling << 1274 scalingAlgorithm = new G4VIonDEDXScalingAlgorithm; 1264 << 1275 1265 G4IonDEDXHandler* handler = << 1276 G4IonDEDXHandler* handler = 1266 new G4IonDEDXHandler(ta << 1277 new G4IonDEDXHandler(table, scalingAlgorithm, name); 1267 1278 1268 lossTableList.push_front(handler); 1279 lossTableList.push_front(handler); 1269 1280 1270 return true; 1281 return true; 1271 } 1282 } 1272 1283 1273 // ########################################## 1284 // ######################################################################### 1274 1285 1275 G4bool G4IonParametrisedLossModel::RemoveDEDX 1286 G4bool G4IonParametrisedLossModel::RemoveDEDXTable( 1276 const G4String& nam) { << 1287 const G4String& name) { 1277 1288 1278 LossTableList::iterator iter = lossTableLis 1289 LossTableList::iterator iter = lossTableList.begin(); 1279 LossTableList::iterator iter_end = lossTabl 1290 LossTableList::iterator iter_end = lossTableList.end(); 1280 << 1291 1281 for(;iter != iter_end; iter++) { 1292 for(;iter != iter_end; iter++) { 1282 G4String tableName = (*iter) -> GetName( 1293 G4String tableName = (*iter) -> GetName(); 1283 1294 1284 if(tableName == nam) { << 1295 if(tableName == name) { 1285 delete (*iter); 1296 delete (*iter); 1286 1297 1287 // Remove from table list 1298 // Remove from table list 1288 lossTableList.erase(iter); 1299 lossTableList.erase(iter); 1289 1300 1290 // Range vs energy and energy vs rang 1301 // Range vs energy and energy vs range vectors are cleared 1291 RangeEnergyTable::iterator iterRange 1302 RangeEnergyTable::iterator iterRange = r.begin(); 1292 RangeEnergyTable::iterator iterRange_ 1303 RangeEnergyTable::iterator iterRange_end = r.end(); 1293 1304 1294 for(;iterRange != iterRange_end; iter << 1305 for(;iterRange != iterRange_end; iterRange++) 1295 d 1306 delete iterRange -> second; 1296 r.clear(); 1307 r.clear(); 1297 1308 1298 EnergyRangeTable::iterator iterEnergy 1309 EnergyRangeTable::iterator iterEnergy = E.begin(); 1299 EnergyRangeTable::iterator iterEnergy 1310 EnergyRangeTable::iterator iterEnergy_end = E.end(); 1300 1311 1301 for(;iterEnergy != iterEnergy_end; it << 1312 for(;iterEnergy != iterEnergy_end; iterEnergy++) 1302 d 1313 delete iterEnergy -> second; 1303 E.clear(); 1314 E.clear(); 1304 1315 1305 return true; 1316 return true; 1306 } 1317 } 1307 } 1318 } 1308 1319 1309 return false; 1320 return false; 1310 } 1321 } >> 1322 >> 1323 // ######################################################################### >> 1324 >> 1325 void G4IonParametrisedLossModel::DeactivateICRU73Scaling() { >> 1326 >> 1327 RemoveDEDXTable("ICRU73"); >> 1328 AddDEDXTable("ICRU73", new G4IonStoppingData("ion_stopping_data/icru73")); >> 1329 } >> 1330 >> 1331 // ######################################################################### 1311 1332