Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 //--------------------------------------------------------------------------- 27 // 28 // ClassName: G4TablesForExtrapolator 29 // 30 // Description: This class provide calculation of energy loss, fluctuation, 31 // and msc angle 32 // 33 // Author: 09.12.04 V.Ivanchenko 34 // 35 // Modification: 36 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko) 37 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko) 38 // 21-03-06 Add verbosity defined in the constructor and Initialisation 39 // start only when first public method is called (V.Ivanchenko) 40 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI) 41 // 12-05-06 SEt linLossLimit=0.001 (VI) 42 // 43 //---------------------------------------------------------------------------- 44 // 45 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 47 48 #include "G4TablesForExtrapolator.hh" 49 #include "G4PhysicalConstants.hh" 50 #include "G4SystemOfUnits.hh" 51 #include "G4LossTableManager.hh" 52 #include "G4PhysicsLogVector.hh" 53 #include "G4ParticleDefinition.hh" 54 #include "G4Material.hh" 55 #include "G4MaterialCutsCouple.hh" 56 #include "G4Electron.hh" 57 #include "G4Positron.hh" 58 #include "G4Proton.hh" 59 #include "G4MuonPlus.hh" 60 #include "G4MuonMinus.hh" 61 #include "G4EmParameters.hh" 62 #include "G4MollerBhabhaModel.hh" 63 #include "G4BetheBlochModel.hh" 64 #include "G4eBremsstrahlungRelModel.hh" 65 #include "G4MuPairProductionModel.hh" 66 #include "G4MuBremsstrahlungModel.hh" 67 #include "G4ProductionCuts.hh" 68 #include "G4LossTableBuilder.hh" 69 #include "G4WentzelVIModel.hh" 70 #include "G4ios.hh" 71 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 73 74 G4TablesForExtrapolator::G4TablesForExtrapolator(G4int verb, G4int bins, 75 G4double e1, G4double e2) 76 : emin(e1), emax(e2), verbose(verb), nbins(bins) 77 { 78 electron = G4Electron::Electron(); 79 positron = G4Positron::Positron(); 80 proton = G4Proton::Proton(); 81 muonPlus = G4MuonPlus::MuonPlus(); 82 muonMinus= G4MuonMinus::MuonMinus(); 83 } 84 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 86 87 G4TablesForExtrapolator::~G4TablesForExtrapolator() 88 { 89 if(nullptr != dedxElectron) { 90 dedxElectron->clearAndDestroy(); 91 delete dedxElectron; 92 } 93 if(nullptr != dedxPositron) { 94 dedxPositron->clearAndDestroy(); 95 delete dedxPositron; 96 } 97 if(nullptr != dedxProton) { 98 dedxProton->clearAndDestroy(); 99 delete dedxProton; 100 } 101 if(nullptr != dedxMuon) { 102 dedxMuon->clearAndDestroy(); 103 delete dedxMuon; 104 } 105 if(nullptr != rangeElectron) { 106 rangeElectron->clearAndDestroy(); 107 delete rangeElectron; 108 } 109 if(nullptr != rangePositron) { 110 rangePositron->clearAndDestroy(); 111 delete rangePositron; 112 } 113 if(nullptr != rangeProton) { 114 rangeProton->clearAndDestroy(); 115 delete rangeProton; 116 } 117 if(nullptr != rangeMuon) { 118 rangeMuon->clearAndDestroy(); 119 delete rangeMuon; 120 } 121 if(nullptr != invRangeElectron) { 122 invRangeElectron->clearAndDestroy(); 123 delete invRangeElectron; 124 } 125 if(nullptr != invRangePositron) { 126 invRangePositron->clearAndDestroy(); 127 delete invRangePositron; 128 } 129 if(nullptr != invRangeProton) { 130 invRangeProton->clearAndDestroy(); 131 delete invRangeProton; 132 } 133 if(nullptr != invRangeMuon) { 134 invRangeMuon->clearAndDestroy(); 135 delete invRangeMuon; 136 } 137 if(nullptr != mscElectron) { 138 mscElectron->clearAndDestroy(); 139 delete mscElectron; 140 } 141 delete pcuts; 142 delete builder; 143 } 144 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 146 147 const G4PhysicsTable* 148 G4TablesForExtrapolator::GetPhysicsTable(ExtTableType type) const 149 { 150 const G4PhysicsTable* table = nullptr; 151 switch (type) 152 { 153 case fDedxElectron: 154 table = dedxElectron; 155 break; 156 case fDedxPositron: 157 table = dedxPositron; 158 break; 159 case fDedxProton: 160 table = dedxProton; 161 break; 162 case fDedxMuon: 163 table = dedxMuon; 164 break; 165 case fRangeElectron: 166 table = rangeElectron; 167 break; 168 case fRangePositron: 169 table = rangePositron; 170 break; 171 case fRangeProton: 172 table = rangeProton; 173 break; 174 case fRangeMuon: 175 table = rangeMuon; 176 break; 177 case fInvRangeElectron: 178 table = invRangeElectron; 179 break; 180 case fInvRangePositron: 181 table = invRangePositron; 182 break; 183 case fInvRangeProton: 184 table = invRangeProton; 185 break; 186 case fInvRangeMuon: 187 table = invRangeMuon; 188 break; 189 case fMscElectron: 190 table = mscElectron; 191 } 192 return table; 193 } 194 195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 196 197 void G4TablesForExtrapolator::Initialisation() 198 { 199 if(verbose>1) { 200 G4cout << "### G4TablesForExtrapolator::Initialisation" << G4endl; 201 } 202 G4int num = (G4int)G4Material::GetNumberOfMaterials(); 203 if(nmat == num) { return; } 204 nmat = num; 205 cuts.resize(nmat, DBL_MAX); 206 couples.resize(nmat, nullptr); 207 208 const G4MaterialTable* mtable = G4Material::GetMaterialTable(); 209 if(!pcuts) { pcuts = new G4ProductionCuts(); } 210 211 for(G4int i=0; i<nmat; ++i) { 212 couples[i] = new G4MaterialCutsCouple((*mtable)[i],pcuts); 213 } 214 215 dedxElectron = PrepareTable(dedxElectron); 216 dedxPositron = PrepareTable(dedxPositron); 217 dedxMuon = PrepareTable(dedxMuon); 218 dedxProton = PrepareTable(dedxProton); 219 rangeElectron = PrepareTable(rangeElectron); 220 rangePositron = PrepareTable(rangePositron); 221 rangeMuon = PrepareTable(rangeMuon); 222 rangeProton = PrepareTable(rangeProton); 223 invRangeElectron = PrepareTable(invRangeElectron); 224 invRangePositron = PrepareTable(invRangePositron); 225 invRangeMuon = PrepareTable(invRangeMuon); 226 invRangeProton = PrepareTable(invRangeProton); 227 mscElectron = PrepareTable(mscElectron); 228 229 builder = new G4LossTableBuilder(true); 230 builder->SetBaseMaterialActive(false); 231 232 if(verbose>1) { 233 G4cout << "### G4TablesForExtrapolator Builds electron tables" 234 << G4endl; 235 } 236 ComputeElectronDEDX(electron, dedxElectron); 237 builder->BuildRangeTable(dedxElectron,rangeElectron); 238 builder->BuildInverseRangeTable(rangeElectron, invRangeElectron); 239 240 if(verbose>1) { 241 G4cout << "### G4TablesForExtrapolator Builds positron tables" 242 << G4endl; 243 } 244 ComputeElectronDEDX(positron, dedxPositron); 245 builder->BuildRangeTable(dedxPositron, rangePositron); 246 builder->BuildInverseRangeTable(rangePositron, invRangePositron); 247 248 if(verbose>1) { 249 G4cout << "### G4TablesForExtrapolator Builds muon tables" << G4endl; 250 } 251 ComputeMuonDEDX(muonPlus, dedxMuon); 252 builder->BuildRangeTable(dedxMuon, rangeMuon); 253 builder->BuildInverseRangeTable(rangeMuon, invRangeMuon); 254 255 if(verbose>2) { 256 G4cout << "DEDX MUON" << G4endl; 257 G4cout << *dedxMuon << G4endl; 258 G4cout << "RANGE MUON" << G4endl; 259 G4cout << *rangeMuon << G4endl; 260 G4cout << "INVRANGE MUON" << G4endl; 261 G4cout << *invRangeMuon << G4endl; 262 } 263 if(verbose>1) { 264 G4cout << "### G4TablesForExtrapolator Builds proton tables" 265 << G4endl; 266 } 267 ComputeProtonDEDX(proton, dedxProton); 268 builder->BuildRangeTable(dedxProton, rangeProton); 269 builder->BuildInverseRangeTable(rangeProton, invRangeProton); 270 271 ComputeTrasportXS(electron, mscElectron); 272 } 273 274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 275 276 G4PhysicsTable* G4TablesForExtrapolator::PrepareTable(G4PhysicsTable* ptr) 277 { 278 G4PhysicsTable* table = ptr; 279 if(nullptr == ptr) { table = new G4PhysicsTable(); } 280 G4int n = (G4int)table->length(); 281 for(G4int i=n; i<nmat; ++i) { 282 G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins, splineFlag); 283 table->push_back(v); 284 } 285 return table; 286 } 287 288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 289 290 void G4TablesForExtrapolator::ComputeElectronDEDX( 291 const G4ParticleDefinition* part, 292 G4PhysicsTable* table) 293 { 294 G4MollerBhabhaModel* ioni = new G4MollerBhabhaModel(); 295 G4eBremsstrahlungRelModel* brem = new G4eBremsstrahlungRelModel(); 296 ioni->Initialise(part, cuts); 297 brem->Initialise(part, cuts); 298 ioni->SetUseBaseMaterials(false); 299 brem->SetUseBaseMaterials(false); 300 301 mass = electron_mass_c2; 302 charge2 = 1.0; 303 currentParticle = part; 304 305 const G4MaterialTable* mtable = G4Material::GetMaterialTable(); 306 if(0<verbose) { 307 G4cout << "G4TablesForExtrapolator::ComputeElectronDEDX for " 308 << part->GetParticleName() 309 << G4endl; 310 } 311 for(G4int i=0; i<nmat; ++i) { 312 313 const G4Material* mat = (*mtable)[i]; 314 if(1<verbose) { 315 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl; 316 } 317 G4PhysicsVector* aVector = (*table)[i]; 318 319 for(G4int j=0; j<=nbins; ++j) { 320 321 G4double e = aVector->Energy(j); 322 G4double dedx = ioni->ComputeDEDXPerVolume(mat,part,e,e) 323 + brem->ComputeDEDXPerVolume(mat,part,e,e); 324 if(1<verbose) { 325 G4cout << "j= " << j 326 << " e(MeV)= " << e/MeV 327 << " dedx(Mev/cm)= " << dedx*cm/MeV 328 << " dedx(Mev.cm2/g)= " 329 << dedx/((MeV*mat->GetDensity())/(g/cm2)) 330 << G4endl; 331 } 332 aVector->PutValue(j,dedx); 333 } 334 if(splineFlag) { aVector->FillSecondDerivatives(); } 335 } 336 } 337 338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 339 340 void 341 G4TablesForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part, 342 G4PhysicsTable* table) 343 { 344 G4BetheBlochModel* ioni = new G4BetheBlochModel(); 345 G4MuPairProductionModel* pair = new G4MuPairProductionModel(part); 346 G4MuBremsstrahlungModel* brem = new G4MuBremsstrahlungModel(part); 347 ioni->Initialise(part, cuts); 348 pair->Initialise(part, cuts); 349 brem->Initialise(part, cuts); 350 ioni->SetUseBaseMaterials(false); 351 pair->SetUseBaseMaterials(false); 352 brem->SetUseBaseMaterials(false); 353 354 mass = part->GetPDGMass(); 355 charge2 = 1.0; 356 currentParticle = part; 357 358 const G4MaterialTable* mtable = G4Material::GetMaterialTable(); 359 360 if(0<verbose) { 361 G4cout << "G4TablesForExtrapolator::ComputeMuonDEDX for " 362 << part->GetParticleName() 363 << G4endl; 364 } 365 366 for(G4int i=0; i<nmat; ++i) { 367 368 const G4Material* mat = (*mtable)[i]; 369 if(1<verbose) { 370 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl; 371 } 372 G4PhysicsVector* aVector = (*table)[i]; 373 for(G4int j=0; j<=nbins; ++j) { 374 375 G4double e = aVector->Energy(j); 376 G4double dedx = ioni->ComputeDEDXPerVolume(mat,part,e,e) + 377 pair->ComputeDEDXPerVolume(mat,part,e,e) + 378 brem->ComputeDEDXPerVolume(mat,part,e,e); 379 aVector->PutValue(j,dedx); 380 if(1<verbose) { 381 G4cout << "j= " << j 382 << " e(MeV)= " << e/MeV 383 << " dedx(Mev/cm)= " << dedx*cm/MeV 384 << " dedx(Mev/(g/cm2)= " 385 << dedx/((MeV*mat->GetDensity())/(g/cm2)) 386 << G4endl; 387 } 388 } 389 if(splineFlag) { aVector->FillSecondDerivatives(); } 390 } 391 } 392 393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 394 395 void 396 G4TablesForExtrapolator::ComputeProtonDEDX(const G4ParticleDefinition* part, 397 G4PhysicsTable* table) 398 { 399 G4BetheBlochModel* ioni = new G4BetheBlochModel(); 400 ioni->Initialise(part, cuts); 401 ioni->SetUseBaseMaterials(false); 402 403 mass = part->GetPDGMass(); 404 charge2 = 1.0; 405 currentParticle = part; 406 407 const G4MaterialTable* mtable = G4Material::GetMaterialTable(); 408 409 if(0<verbose) { 410 G4cout << "G4TablesForExtrapolator::ComputeProtonDEDX for " 411 << part->GetParticleName() 412 << G4endl; 413 } 414 415 for(G4int i=0; i<nmat; ++i) { 416 417 const G4Material* mat = (*mtable)[i]; 418 if(1<verbose) 419 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl; 420 G4PhysicsVector* aVector = (*table)[i]; 421 for(G4int j=0; j<=nbins; ++j) { 422 423 G4double e = aVector->Energy(j); 424 G4double dedx = ioni->ComputeDEDXPerVolume(mat,part,e,e); 425 aVector->PutValue(j,dedx); 426 if(1<verbose) { 427 G4cout << "j= " << j 428 << " e(MeV)= " << e/MeV 429 << " dedx(Mev/cm)= " << dedx*cm/MeV 430 << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2)) 431 << G4endl; 432 } 433 } 434 if(splineFlag) { aVector->FillSecondDerivatives(); } 435 } 436 } 437 438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 439 440 void 441 G4TablesForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part, 442 G4PhysicsTable* table) 443 { 444 G4WentzelVIModel* msc = new G4WentzelVIModel(); 445 msc->SetPolarAngleLimit(CLHEP::pi); 446 msc->Initialise(part, cuts); 447 msc->SetUseBaseMaterials(false); 448 449 mass = part->GetPDGMass(); 450 charge2 = 1.0; 451 currentParticle = part; 452 453 const G4MaterialTable* mtable = G4Material::GetMaterialTable(); 454 455 if(0<verbose) { 456 G4cout << "G4TablesForExtrapolator::ComputeTransportXS for " 457 << part->GetParticleName() 458 << G4endl; 459 } 460 461 for(G4int i=0; i<nmat; ++i) { 462 463 const G4Material* mat = (*mtable)[i]; 464 msc->SetCurrentCouple(couples[i]); 465 if(1<verbose) 466 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl; 467 G4PhysicsVector* aVector = (*table)[i]; 468 for(G4int j=0; j<=nbins; ++j) { 469 470 G4double e = aVector->Energy(j); 471 G4double xs = msc->CrossSectionPerVolume(mat,part,e); 472 aVector->PutValue(j,xs); 473 if(1<verbose) { 474 G4cout << "j= " << j << " e(MeV)= " << e/MeV 475 << " xs(1/mm)= " << xs*mm << G4endl; 476 } 477 } 478 if(splineFlag) { aVector->FillSecondDerivatives(); } 479 } 480 } 481 482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 483 484