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: G4EnergyLossForExtrapolator.cc,v 1.21 2010-11-04 12:40:29 vnivanch Exp $ >> 27 // GEANT4 tag $Name: not supported by cvs2svn $ 26 // 28 // 27 //-------------------------------------------- 29 //--------------------------------------------------------------------------- 28 // 30 // 29 // ClassName: G4EnergyLossForExtrapolator 31 // ClassName: G4EnergyLossForExtrapolator 30 // 32 // 31 // Description: This class provide calculatio 33 // Description: This class provide calculation of energy loss, fluctuation, 32 // and msc angle 34 // and msc angle 33 // 35 // 34 // Author: 09.12.04 V.Ivanchenko 36 // Author: 09.12.04 V.Ivanchenko 35 // 37 // 36 // Modification: 38 // Modification: 37 // 08-04-05 Rename Propogator -> Extrapolator 39 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko) 38 // 16-03-06 Add muon tables and fix bug in uni 40 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko) 39 // 21-03-06 Add verbosity defined in the const 41 // 21-03-06 Add verbosity defined in the constructor and Initialisation 40 // start only when first public metho 42 // start only when first public method is called (V.Ivanchenko) 41 // 03-05-06 Remove unused pointer G4Material* 43 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI) 42 // 12-05-06 SEt linLossLimit=0.001 (VI) 44 // 12-05-06 SEt linLossLimit=0.001 (VI) 43 // 45 // 44 //-------------------------------------------- 46 //---------------------------------------------------------------------------- 45 // 47 // 46 48 47 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 48 50 49 #include "G4EnergyLossForExtrapolator.hh" 51 #include "G4EnergyLossForExtrapolator.hh" 50 #include "G4PhysicalConstants.hh" << 52 #include "G4PhysicsLogVector.hh" 51 #include "G4SystemOfUnits.hh" << 52 #include "G4ParticleDefinition.hh" 53 #include "G4ParticleDefinition.hh" 53 #include "G4Material.hh" 54 #include "G4Material.hh" 54 #include "G4MaterialCutsCouple.hh" 55 #include "G4MaterialCutsCouple.hh" 55 #include "G4Electron.hh" 56 #include "G4Electron.hh" 56 #include "G4Positron.hh" 57 #include "G4Positron.hh" 57 #include "G4Proton.hh" 58 #include "G4Proton.hh" 58 #include "G4MuonPlus.hh" 59 #include "G4MuonPlus.hh" 59 #include "G4MuonMinus.hh" 60 #include "G4MuonMinus.hh" 60 #include "G4ParticleTable.hh" 61 #include "G4ParticleTable.hh" >> 62 #include "G4LossTableBuilder.hh" >> 63 #include "G4MollerBhabhaModel.hh" >> 64 #include "G4BetheBlochModel.hh" >> 65 #include "G4eBremsstrahlungModel.hh" >> 66 #include "G4MuPairProductionModel.hh" >> 67 #include "G4MuBremsstrahlungModel.hh" >> 68 #include "G4ProductionCuts.hh" >> 69 #include "G4LossTableManager.hh" >> 70 #include "G4WentzelVIModel.hh" 61 71 62 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 63 73 64 #ifdef G4MULTITHREADED << 65 G4Mutex G4EnergyLossForExtrapolator::extrMutex << 66 #endif << 67 << 68 G4TablesForExtrapolator* G4EnergyLossForExtrap << 69 << 70 G4EnergyLossForExtrapolator::G4EnergyLossForEx 74 G4EnergyLossForExtrapolator::G4EnergyLossForExtrapolator(G4int verb) 71 : maxEnergyTransfer(DBL_MAX), verbose(verb) << 75 :maxEnergyTransfer(DBL_MAX),verbose(verb),isInitialised(false) 72 { 76 { 73 emin = 1.*CLHEP::MeV; << 77 currentParticle = 0; 74 emax = 100.*CLHEP::TeV; << 78 currentMaterial = 0; >> 79 >> 80 linLossLimit = 0.01; >> 81 emin = 1.*MeV; >> 82 emax = 10.*TeV; >> 83 nbins = 70; >> 84 >> 85 nmat = index = 0; >> 86 cuts = 0; >> 87 >> 88 mass = charge2 = electronDensity = radLength = bg2 = beta2 >> 89 = kineticEnergy = tmax = 0; >> 90 gam = 1.0; >> 91 >> 92 dedxElectron = dedxPositron = dedxProton = rangeElectron >> 93 = rangePositron = rangeProton = invRangeElectron = invRangePositron >> 94 = invRangeProton = mscElectron = dedxMuon = rangeMuon = invRangeMuon = 0; >> 95 cuts = 0; >> 96 electron = positron = proton = muonPlus = muonMinus = 0; 75 } 97 } 76 98 77 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 78 100 79 G4EnergyLossForExtrapolator::~G4EnergyLossForE << 101 G4EnergyLossForExtrapolator:: ~G4EnergyLossForExtrapolator() 80 { 102 { 81 if(isMaster) { << 103 for(G4int i=0; i<nmat; i++) {delete couples[i];} 82 delete tables; << 104 delete dedxElectron; 83 tables = nullptr; << 105 delete dedxPositron; 84 } << 106 delete dedxProton; >> 107 delete dedxMuon; >> 108 delete rangeElectron; >> 109 delete rangePositron; >> 110 delete rangeProton; >> 111 delete rangeMuon; >> 112 delete invRangeElectron; >> 113 delete invRangePositron; >> 114 delete invRangeProton; >> 115 delete invRangeMuon; >> 116 delete mscElectron; >> 117 delete cuts; 85 } 118 } 86 119 87 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 88 121 89 G4double << 122 G4double G4EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy, 90 G4EnergyLossForExtrapolator::EnergyAfterStep(G << 123 G4double stepLength, 91 G4double stepLength, << 124 const G4Material* mat, 92 const G4Material* mat, << 125 const G4ParticleDefinition* part) 93 const G4ParticleDefinition* par << 94 { 126 { >> 127 if(!isInitialised) Initialisation(); 95 G4double kinEnergyFinal = kinEnergy; 128 G4double kinEnergyFinal = kinEnergy; 96 if(SetupKinematics(part, mat, kinEnergy)) { 129 if(SetupKinematics(part, mat, kinEnergy)) { 97 G4double step = TrueStepLength(kinEnergy,s 130 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part); 98 G4double r = ComputeRange(kinEnergy,part, << 131 G4double r = ComputeRange(kinEnergy,part); 99 if(r <= step) { 132 if(r <= step) { 100 kinEnergyFinal = 0.0; 133 kinEnergyFinal = 0.0; 101 } else if(step < linLossLimit*r) { 134 } else if(step < linLossLimit*r) { 102 kinEnergyFinal -= step*ComputeDEDX(kinEn << 135 kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part); 103 } else { 136 } else { 104 G4double r1 = r - step; 137 G4double r1 = r - step; 105 kinEnergyFinal = ComputeEnergy(r1,part,m << 138 kinEnergyFinal = ComputeEnergy(r1,part); 106 } 139 } 107 } 140 } 108 return kinEnergyFinal; 141 return kinEnergyFinal; 109 } 142 } 110 143 111 //....oooOO0OOooo........oooOO0OOooo........oo 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 112 145 113 G4double << 146 G4double G4EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy, 114 G4EnergyLossForExtrapolator::EnergyBeforeStep( << 147 G4double stepLength, 115 G4double stepLength, << 148 const G4Material* mat, 116 const G4Material* mat, << 149 const G4ParticleDefinition* part) 117 const G4ParticleDefinition* pa << 118 { 150 { 119 //G4cout << "G4EnergyLossForExtrapolator::En << 151 if(!isInitialised) Initialisation(); 120 G4double kinEnergyFinal = kinEnergy; 152 G4double kinEnergyFinal = kinEnergy; 121 153 122 if(SetupKinematics(part, mat, kinEnergy)) { 154 if(SetupKinematics(part, mat, kinEnergy)) { 123 G4double step = TrueStepLength(kinEnergy,s 155 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part); 124 G4double r = ComputeRange(kinEnergy,part,m << 156 G4double r = ComputeRange(kinEnergy,part); 125 157 126 if(step < linLossLimit*r) { 158 if(step < linLossLimit*r) { 127 kinEnergyFinal += step*ComputeDEDX(kinEn << 159 kinEnergyFinal += step*ComputeDEDX(kinEnergy,part); 128 } else { 160 } else { 129 G4double r1 = r + step; 161 G4double r1 = r + step; 130 kinEnergyFinal = ComputeEnergy(r1,part,m << 162 kinEnergyFinal = ComputeEnergy(r1,part); 131 } 163 } 132 } 164 } 133 return kinEnergyFinal; 165 return kinEnergyFinal; 134 } 166 } 135 167 136 //....oooOO0OOooo........oooOO0OOooo........oo 168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 137 169 138 G4double << 170 G4double G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy, 139 G4EnergyLossForExtrapolator::TrueStepLength(G4 << 171 G4double stepLength, 140 G4double stepLength, << 172 const G4Material* mat, 141 const G4Material* mat, << 173 const G4ParticleDefinition* part) 142 const G4ParticleDefinition* part << 143 { 174 { 144 G4double res = stepLength; 175 G4double res = stepLength; 145 //G4cout << "## G4EnergyLossForExtrapolator: << 176 if(!isInitialised) Initialisation(); 146 // << " " << part->GetParticleName() << << 147 if(SetupKinematics(part, mat, kinEnergy)) { 177 if(SetupKinematics(part, mat, kinEnergy)) { 148 if(part == electron || part == positron) { 178 if(part == electron || part == positron) { 149 const G4double x = stepLength* << 179 G4double x = stepLength*ComputeValue(kinEnergy, mscElectron); 150 ComputeValue(kinEnergy, GetPhysicsTable(fMsc << 180 if(x < 0.2) res *= (1.0 + 0.5*x + x*x/3.0); 151 //G4cout << " x= " << x << G4endl; << 181 else if(x < 0.9999) res = -std::log(1.0 - x)*stepLength/x; 152 if(x < 0.2) { res *= (1.0 + 0.5* << 182 else res = ComputeRange(kinEnergy,part); 153 else if(x < 0.9999) { res = -G4Log(1.0 - << 183 154 else { res = ComputeRange(kinEnergy, par << 155 } else { 184 } else { 156 res = ComputeTrueStep(mat,part,kinEnergy 185 res = ComputeTrueStep(mat,part,kinEnergy,stepLength); 157 } 186 } 158 } 187 } 159 return res; 188 return res; 160 } 189 } 161 190 162 //....oooOO0OOooo........oooOO0OOooo........oo 191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 163 192 164 G4bool << 193 G4bool G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part, 165 G4EnergyLossForExtrapolator::SetupKinematics(c << 194 const G4Material* mat, 166 const G4Material* mat, << 195 G4double kinEnergy) 167 G4double kinEnergy) << 196 { 168 { << 197 if(!part || !mat || kinEnergy < keV) return false; 169 if(mat->GetNumberOfMaterials() != nmat) { In << 198 if(!isInitialised) Initialisation(); 170 if(nullptr == part || nullptr == mat || kinE << 199 G4bool flag = false; 171 { return false; } << 172 if(part != currentParticle) { 200 if(part != currentParticle) { >> 201 flag = true; 173 currentParticle = part; 202 currentParticle = part; >> 203 mass = part->GetPDGMass(); 174 G4double q = part->GetPDGCharge()/eplus; 204 G4double q = part->GetPDGCharge()/eplus; 175 charge2 = q*q; 205 charge2 = q*q; 176 } 206 } 177 if(mat != currentMaterial) { 207 if(mat != currentMaterial) { 178 size_t i = mat->GetIndex(); << 208 G4int i = mat->GetIndex(); 179 if(i >= nmat) { 209 if(i >= nmat) { 180 G4cout << "### G4EnergyLossForExtrapolat << 210 G4cout << "### G4EnergyLossForExtrapolator WARNING:index i= " 181 << i << " above number of materials " < << 211 << i << " is out of table - NO extrapolation" << G4endl; 182 return false; << 183 } else { 212 } else { >> 213 flag = true; 184 currentMaterial = mat; 214 currentMaterial = mat; 185 electronDensity = mat->GetElectronDensit 215 electronDensity = mat->GetElectronDensity(); 186 radLength = mat->GetRadlen(); 216 radLength = mat->GetRadlen(); >> 217 index = i; 187 } 218 } 188 } 219 } 189 if(kinEnergy != kineticEnergy) { << 220 if(flag || kinEnergy != kineticEnergy) { 190 kineticEnergy = kinEnergy; 221 kineticEnergy = kinEnergy; 191 G4double mass = part->GetPDGMass(); << 192 G4double tau = kinEnergy/mass; 222 G4double tau = kinEnergy/mass; 193 223 194 gam = tau + 1.0; 224 gam = tau + 1.0; 195 bg2 = tau * (tau + 2.0); 225 bg2 = tau * (tau + 2.0); 196 beta2 = bg2/(gam*gam); 226 beta2 = bg2/(gam*gam); 197 tmax = kinEnergy; 227 tmax = kinEnergy; 198 if(part == electron) tmax *= 0.5; 228 if(part == electron) tmax *= 0.5; 199 else if(part != positron) { 229 else if(part != positron) { 200 G4double r = CLHEP::electron_mass_c2/mas << 230 G4double r = electron_mass_c2/mass; 201 tmax = 2.0*bg2*CLHEP::electron_mass_c2/( << 231 tmax = 2.0*bg2*electron_mass_c2/(1.0 + 2.0*gam*r + r*r); 202 } 232 } 203 tmax = std::min(tmax, maxEnergyTransfer); << 233 if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer; 204 } 234 } 205 return true; 235 return true; 206 } 236 } 207 237 208 //....oooOO0OOooo........oooOO0OOooo........oo 238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 209 239 210 const G4ParticleDefinition* << 240 void G4EnergyLossForExtrapolator::Initialisation() 211 G4EnergyLossForExtrapolator::FindParticle(cons << 212 { 241 { 213 currentParticle = G4ParticleTable::GetPartic << 242 isInitialised = true; 214 if(nullptr == currentParticle) { << 243 if(verbose>1) { 215 G4cout << "### G4EnergyLossForExtrapolator << 244 G4cout << "### G4EnergyLossForExtrapolator::Initialisation" << G4endl; 216 << "FindParticle() fails to find " << nam << 245 } >> 246 currentParticle = 0; >> 247 currentMaterial = 0; >> 248 kineticEnergy = 0.0; >> 249 >> 250 electron = G4Electron::Electron(); >> 251 positron = G4Positron::Positron(); >> 252 proton = G4Proton::Proton(); >> 253 muonPlus = G4MuonPlus::MuonPlus(); >> 254 muonMinus= G4MuonMinus::MuonMinus(); >> 255 >> 256 currentParticleName = ""; >> 257 >> 258 nmat = G4Material::GetNumberOfMaterials(); >> 259 const G4MaterialTable* mtable = G4Material::GetMaterialTable(); >> 260 cuts = new G4ProductionCuts(); >> 261 >> 262 const G4MaterialCutsCouple* couple; >> 263 for(G4int i=0; i<nmat; i++) { >> 264 couple = new G4MaterialCutsCouple((*mtable)[i],cuts); >> 265 couples.push_back(couple); 217 } 266 } 218 return currentParticle; << 267 >> 268 dedxElectron = PrepareTable(); >> 269 dedxPositron = PrepareTable(); >> 270 dedxMuon = PrepareTable(); >> 271 dedxProton = PrepareTable(); >> 272 rangeElectron = PrepareTable(); >> 273 rangePositron = PrepareTable(); >> 274 rangeMuon = PrepareTable(); >> 275 rangeProton = PrepareTable(); >> 276 invRangeElectron = PrepareTable(); >> 277 invRangePositron = PrepareTable(); >> 278 invRangeMuon = PrepareTable(); >> 279 invRangeProton = PrepareTable(); >> 280 mscElectron = PrepareTable(); >> 281 >> 282 G4LossTableBuilder builder; >> 283 >> 284 if(verbose>1) >> 285 G4cout << "### G4EnergyLossForExtrapolator Builds electron tables" << G4endl; >> 286 >> 287 ComputeElectronDEDX(electron, dedxElectron); >> 288 builder.BuildRangeTable(dedxElectron,rangeElectron); >> 289 builder.BuildInverseRangeTable(rangeElectron, invRangeElectron); >> 290 >> 291 if(verbose>1) >> 292 G4cout << "### G4EnergyLossForExtrapolator Builds positron tables" << G4endl; >> 293 >> 294 ComputeElectronDEDX(positron, dedxPositron); >> 295 builder.BuildRangeTable(dedxPositron, rangePositron); >> 296 builder.BuildInverseRangeTable(rangePositron, invRangePositron); >> 297 >> 298 if(verbose>1) >> 299 G4cout << "### G4EnergyLossForExtrapolator Builds muon tables" << G4endl; >> 300 >> 301 ComputeMuonDEDX(muonPlus, dedxMuon); >> 302 builder.BuildRangeTable(dedxMuon, rangeMuon); >> 303 builder.BuildInverseRangeTable(rangeMuon, invRangeMuon); >> 304 >> 305 if(verbose>1) >> 306 G4cout << "### G4EnergyLossForExtrapolator Builds proton tables" << G4endl; >> 307 >> 308 ComputeProtonDEDX(proton, dedxProton); >> 309 builder.BuildRangeTable(dedxProton, rangeProton); >> 310 builder.BuildInverseRangeTable(rangeProton, invRangeProton); >> 311 >> 312 ComputeTrasportXS(electron, mscElectron); >> 313 } >> 314 >> 315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 316 >> 317 G4PhysicsTable* G4EnergyLossForExtrapolator::PrepareTable() >> 318 { >> 319 G4PhysicsTable* table = new G4PhysicsTable(); >> 320 >> 321 for(G4int i=0; i<nmat; i++) { >> 322 >> 323 G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins); >> 324 v->SetSpline(G4LossTableManager::Instance()->SplineFlag()); >> 325 table->push_back(v); >> 326 } >> 327 return table; 219 } 328 } 220 329 221 //....oooOO0OOooo........oooOO0OOooo........oo 330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 222 331 223 G4double << 332 const G4ParticleDefinition* G4EnergyLossForExtrapolator::FindParticle(const G4String& name) 224 G4EnergyLossForExtrapolator::ComputeDEDX(G4dou << 333 { 225 const G4ParticleDefinition* part, << 334 const G4ParticleDefinition* p = 0; 226 const << 335 if(name != currentParticleName) { >> 336 p = G4ParticleTable::GetParticleTable()->FindParticle(name); >> 337 if(!p) { >> 338 G4cout << "### G4EnergyLossForExtrapolator WARNING:FindParticle fails to find " >> 339 << name << G4endl; >> 340 } >> 341 } else { >> 342 p = currentParticle; >> 343 } >> 344 return p; >> 345 } >> 346 >> 347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 348 >> 349 G4double G4EnergyLossForExtrapolator::ComputeDEDX(G4double kinEnergy, >> 350 const G4ParticleDefinition* part) 227 { 351 { 228 if(mat->GetNumberOfMaterials() != nmat) { In << 229 G4double x = 0.0; 352 G4double x = 0.0; 230 if(part == electron) { << 353 if(part == electron) x = ComputeValue(kinEnergy, dedxElectron); 231 x = ComputeValue(ekin, GetPhysicsTable(fDe << 354 else if(part == positron) x = ComputeValue(kinEnergy, dedxPositron); 232 } else if(part == positron) { << 355 else if(part == muonPlus || part == muonMinus) { 233 x = ComputeValue(ekin, GetPhysicsTable(fDe << 356 x = ComputeValue(kinEnergy, dedxMuon); 234 } else if(part == muonPlus || part == muonMi << 235 x = ComputeValue(ekin, GetPhysicsTable(fDe << 236 } else { 357 } else { 237 G4double e = ekin*CLHEP::proton_mass_c2/pa << 358 G4double e = kinEnergy*proton_mass_c2/mass; 238 G4double q = part->GetPDGCharge()/CLHEP::e << 359 x = ComputeValue(e, dedxProton)*charge2; 239 x = ComputeValue(e, GetPhysicsTable(fDedxP << 240 } 360 } 241 return x; 361 return x; 242 } 362 } 243 363 244 //....oooOO0OOooo........oooOO0OOooo........oo 364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 245 365 246 G4double << 366 G4double G4EnergyLossForExtrapolator::ComputeRange(G4double kinEnergy, 247 G4EnergyLossForExtrapolator::ComputeRange(G4do << 367 const G4ParticleDefinition* part) 248 const G4ParticleDefinition* part, << 249 const G4Material* mat) << 250 { 368 { 251 if(mat->GetNumberOfMaterials() != nmat) { In << 252 G4double x = 0.0; 369 G4double x = 0.0; 253 if(part == electron) { << 370 if(part == electron) x = ComputeValue(kinEnergy, rangeElectron); 254 x = ComputeValue(ekin, GetPhysicsTable(fRa << 371 else if(part == positron) x = ComputeValue(kinEnergy, rangePositron); 255 } else if(part == positron) { << 372 else if(part == muonPlus || part == muonMinus) 256 x = ComputeValue(ekin, GetPhysicsTable(fRa << 373 x = ComputeValue(kinEnergy, rangeMuon); 257 } else if(part == muonPlus || part == muonMi << 374 else { 258 x = ComputeValue(ekin, GetPhysicsTable(fRa << 375 G4double massratio = proton_mass_c2/mass; 259 } else { << 376 G4double e = kinEnergy*massratio; 260 G4double massratio = CLHEP::proton_mass_c2 << 377 x = ComputeValue(e, rangeProton)/(charge2*massratio); 261 G4double e = ekin*massratio; << 262 G4double q = part->GetPDGCharge()/CLHEP::e << 263 x = ComputeValue(e, GetPhysicsTable(fRange << 264 /(q*q*massratio); << 265 } 378 } 266 return x; 379 return x; 267 } 380 } 268 381 269 //....oooOO0OOooo........oooOO0OOooo........oo 382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 270 383 271 G4double << 384 G4double G4EnergyLossForExtrapolator::ComputeEnergy(G4double range, 272 G4EnergyLossForExtrapolator::ComputeEnergy(G4d << 385 const G4ParticleDefinition* part) 273 const G4ParticleDefinition* part, << 274 const G4Material* mat) << 275 { 386 { 276 if(mat->GetNumberOfMaterials() != nmat) { In << 277 G4double x = 0.0; 387 G4double x = 0.0; 278 if(part == electron) { << 388 if(part == electron) x = ComputeValue(range, invRangeElectron); 279 x = ComputeValue(range,GetPhysicsTable(fIn << 389 else if(part == positron) x = ComputeValue(range, invRangePositron); 280 } else if(part == positron) { << 390 else if(part == muonPlus || part == muonMinus) 281 x = ComputeValue(range,GetPhysicsTable(fIn << 391 x = ComputeValue(range, invRangeMuon); 282 } else if(part == muonPlus || part == muonMi << 392 else { 283 x = ComputeValue(range, GetPhysicsTable(fI << 393 G4double massratio = proton_mass_c2/mass; 284 } else { << 394 G4double r = range*massratio*charge2; 285 G4double massratio = CLHEP::proton_mass_c2 << 395 x = ComputeValue(r, invRangeProton)/massratio; 286 G4double q = part->GetPDGCharge()/CLHEP::e << 287 G4double r = range*massratio*q*q; << 288 x = ComputeValue(r, GetPhysicsTable(fInvRa << 289 } 396 } 290 return x; 397 return x; 291 } 398 } 292 399 293 //....oooOO0OOooo........oooOO0OOooo........oo 400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 294 << 401 295 G4double << 402 void G4EnergyLossForExtrapolator::ComputeElectronDEDX(const G4ParticleDefinition* part, 296 G4EnergyLossForExtrapolator::EnergyDispersion( << 403 G4PhysicsTable* table) 297 G4double stepLength, << 298 const G4Material* mat, << 299 const G4ParticleDefinition* pa << 300 { 404 { 301 G4double sig2 = 0.0; << 405 G4DataVector v; 302 if(SetupKinematics(part, mat, kinEnergy)) { << 406 G4MollerBhabhaModel* ioni = new G4MollerBhabhaModel(); 303 G4double step = ComputeTrueStep(mat,part,k << 407 G4eBremsstrahlungModel* brem = new G4eBremsstrahlungModel(); 304 sig2 = (1.0/beta2 - 0.5) << 408 ioni->Initialise(part, v); 305 *CLHEP::twopi_mc2_rcl2*tmax*step*electro << 409 brem->Initialise(part, v); >> 410 >> 411 mass = electron_mass_c2; >> 412 charge2 = 1.0; >> 413 currentParticle = part; >> 414 >> 415 const G4MaterialTable* mtable = G4Material::GetMaterialTable(); >> 416 if(0<verbose) { >> 417 G4cout << "G4EnergyLossForExtrapolator::ComputeElectronDEDX for " >> 418 << part->GetParticleName() >> 419 << G4endl; >> 420 } >> 421 for(G4int i=0; i<nmat; i++) { >> 422 >> 423 const G4Material* mat = (*mtable)[i]; >> 424 if(1<verbose) >> 425 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl; >> 426 const G4MaterialCutsCouple* couple = couples[i]; >> 427 G4PhysicsVector* aVector = (*table)[i]; >> 428 >> 429 for(G4int j=0; j<=nbins; j++) { >> 430 >> 431 G4double e = aVector->Energy(j); >> 432 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) + brem->ComputeDEDX(couple,part,e,e); >> 433 if(1<verbose) { >> 434 G4cout << "j= " << j >> 435 << " e(MeV)= " << e/MeV >> 436 << " dedx(Mev/cm)= " << dedx*cm/MeV >> 437 << " dedx(Mev.cm2/g)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl; >> 438 } >> 439 aVector->PutValue(j,dedx); >> 440 } >> 441 } >> 442 delete ioni; >> 443 delete brem; >> 444 } >> 445 >> 446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 447 >> 448 void G4EnergyLossForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part, >> 449 G4PhysicsTable* table) >> 450 { >> 451 G4DataVector v; >> 452 G4BetheBlochModel* ioni = new G4BetheBlochModel(); >> 453 G4MuPairProductionModel* pair = new G4MuPairProductionModel(); >> 454 G4MuBremsstrahlungModel* brem = new G4MuBremsstrahlungModel(); >> 455 ioni->Initialise(part, v); >> 456 pair->Initialise(part, v); >> 457 brem->Initialise(part, v); >> 458 >> 459 mass = part->GetPDGMass(); >> 460 charge2 = 1.0; >> 461 currentParticle = part; >> 462 >> 463 const G4MaterialTable* mtable = G4Material::GetMaterialTable(); >> 464 >> 465 if(0<verbose) { >> 466 G4cout << "G4EnergyLossForExtrapolator::ComputeMuonDEDX for " << part->GetParticleName() >> 467 << G4endl; >> 468 } >> 469 >> 470 for(G4int i=0; i<nmat; i++) { >> 471 >> 472 const G4Material* mat = (*mtable)[i]; >> 473 if(1<verbose) >> 474 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl; >> 475 const G4MaterialCutsCouple* couple = couples[i]; >> 476 G4PhysicsVector* aVector = (*table)[i]; >> 477 for(G4int j=0; j<=nbins; j++) { >> 478 >> 479 G4double e = aVector->Energy(j); >> 480 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) + >> 481 pair->ComputeDEDX(couple,part,e,e) + >> 482 brem->ComputeDEDX(couple,part,e,e); >> 483 aVector->PutValue(j,dedx); >> 484 if(1<verbose) { >> 485 G4cout << "j= " << j >> 486 << " e(MeV)= " << e/MeV >> 487 << " dedx(Mev/cm)= " << dedx*cm/MeV >> 488 << " dedx(Mev/(g/cm2)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl; >> 489 } >> 490 } 306 } 491 } 307 return sig2; << 492 delete ioni; 308 } 493 } 309 494 310 //....oooOO0OOooo........oooOO0OOooo........oo 495 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 311 496 312 G4double G4EnergyLossForExtrapolator::AverageS << 497 void G4EnergyLossForExtrapolator::ComputeProtonDEDX(const G4ParticleDefinition* part, 313 G4double kinEnergy, << 498 G4PhysicsTable* table) 314 G4double stepLength, << 315 const G4Material* mat, << 316 const G4ParticleDefinition* part) << 317 { 499 { 318 G4double theta = 0.0; << 500 G4DataVector v; 319 if(SetupKinematics(part, mat, kinEnergy)) { << 501 G4BetheBlochModel* ioni = new G4BetheBlochModel(); 320 G4double t = stepLength/radLength; << 502 ioni->Initialise(part, v); 321 G4double y = std::max(0.001, t); << 503 322 theta = 19.23*CLHEP::MeV*std::sqrt(charge2 << 504 mass = part->GetPDGMass(); 323 /(beta2*gam*part->GetPDGMass()); << 505 charge2 = 1.0; >> 506 currentParticle = part; >> 507 >> 508 const G4MaterialTable* mtable = G4Material::GetMaterialTable(); >> 509 >> 510 if(0<verbose) { >> 511 G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName() >> 512 << G4endl; >> 513 } >> 514 >> 515 for(G4int i=0; i<nmat; i++) { >> 516 >> 517 const G4Material* mat = (*mtable)[i]; >> 518 if(1<verbose) >> 519 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl; >> 520 const G4MaterialCutsCouple* couple = couples[i]; >> 521 G4PhysicsVector* aVector = (*table)[i]; >> 522 for(G4int j=0; j<=nbins; j++) { >> 523 >> 524 G4double e = aVector->Energy(j); >> 525 G4double dedx = ioni->ComputeDEDX(couple,part,e,e); >> 526 aVector->PutValue(j,dedx); >> 527 if(1<verbose) { >> 528 G4cout << "j= " << j >> 529 << " e(MeV)= " << e/MeV >> 530 << " dedx(Mev/cm)= " << dedx*cm/MeV >> 531 << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2)) << G4endl; >> 532 } >> 533 } 324 } 534 } 325 return theta; << 535 delete ioni; 326 } 536 } 327 537 328 //....oooOO0OOooo........oooOO0OOooo........oo 538 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 329 539 330 void G4EnergyLossForExtrapolator::Initialisati << 540 void G4EnergyLossForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part, >> 541 G4PhysicsTable* table) 331 { 542 { 332 if(verbose>0) { << 543 G4DataVector v; 333 G4cout << "### G4EnergyLossForExtrapolator << 544 G4WentzelVIModel* msc = new G4WentzelVIModel(); 334 << tables << G4endl; << 545 msc->SetPolarAngleLimit(CLHEP::pi); >> 546 msc->Initialise(part, v); >> 547 >> 548 mass = part->GetPDGMass(); >> 549 charge2 = 1.0; >> 550 currentParticle = part; >> 551 >> 552 const G4MaterialTable* mtable = G4Material::GetMaterialTable(); >> 553 >> 554 if(0<verbose) { >> 555 G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName() >> 556 << G4endl; 335 } 557 } 336 electron = G4Electron::Electron(); << 558 337 positron = G4Positron::Positron(); << 559 for(G4int i=0; i<nmat; i++) { 338 proton = G4Proton::Proton(); << 339 muonPlus = G4MuonPlus::MuonPlus(); << 340 muonMinus= G4MuonMinus::MuonMinus(); << 341 560 342 // initialisation for the 1st run << 561 const G4Material* mat = (*mtable)[i]; 343 if(nullptr == tables) { << 562 if(1<verbose) 344 #ifdef G4MULTITHREADED << 563 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl; 345 G4MUTEXLOCK(&extrMutex); << 564 G4PhysicsVector* aVector = (*table)[i]; 346 if(nullptr == tables) { << 565 for(G4int j=0; j<=nbins; j++) { 347 #endif << 566 348 isMaster = true; << 567 G4double e = aVector->Energy(j); 349 tables = new G4TablesForExtrapolator(ver << 568 G4double xs = msc->CrossSectionPerVolume(mat,part,e); 350 tables->Initialisation(); << 569 aVector->PutValue(j,xs); 351 nmat = G4Material::GetNumberOfMaterials( << 570 if(1<verbose) { 352 if(verbose > 0) { << 571 G4cout << "j= " << j << " e(MeV)= " << e/MeV 353 G4cout << "### G4EnergyLossForExtrapol << 572 << " xs(1/mm)= " << xs*mm << G4endl; 354 << nmat << " materials Nbins= " << 573 } 355 << nbins << " Emin(MeV)= " << e << 356 << G4endl; << 357 } << 358 #ifdef G4MULTITHREADED << 359 } 574 } 360 G4MUTEXUNLOCK(&extrMutex); << 361 #endif << 362 } << 363 << 364 // initialisation for the next run << 365 if(isMaster && G4Material::GetNumberOfMateri << 366 #ifdef G4MULTITHREADED << 367 G4MUTEXLOCK(&extrMutex); << 368 #endif << 369 tables->Initialisation(); << 370 #ifdef G4MULTITHREADED << 371 G4MUTEXUNLOCK(&extrMutex); << 372 #endif << 373 } 575 } 374 nmat = G4Material::GetNumberOfMaterials(); << 576 delete msc; 375 } 577 } 376 578 377 //....oooOO0OOooo........oooOO0OOooo........oo 579 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 580 378 581