00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065 #include "G4EmCalculator.hh"
00066 #include "G4SystemOfUnits.hh"
00067 #include "G4LossTableManager.hh"
00068 #include "G4VEmProcess.hh"
00069 #include "G4VEnergyLossProcess.hh"
00070 #include "G4VMultipleScattering.hh"
00071 #include "G4Material.hh"
00072 #include "G4MaterialCutsCouple.hh"
00073 #include "G4ParticleDefinition.hh"
00074 #include "G4ParticleTable.hh"
00075 #include "G4PhysicsTable.hh"
00076 #include "G4ProductionCutsTable.hh"
00077 #include "G4ProcessManager.hh"
00078 #include "G4ionEffectiveCharge.hh"
00079 #include "G4RegionStore.hh"
00080 #include "G4Element.hh"
00081 #include "G4EmCorrections.hh"
00082 #include "G4GenericIon.hh"
00083 #include "G4ProcessVector.hh"
00084 #include "G4Gamma.hh"
00085
00086
00087
00088 G4EmCalculator::G4EmCalculator()
00089 {
00090 manager = G4LossTableManager::Instance();
00091 corr = manager->EmCorrections();
00092 nLocalMaterials = 0;
00093 verbose = 0;
00094 currentCoupleIndex = 0;
00095 currentCouple = 0;
00096 currentMaterial = 0;
00097 currentParticle = 0;
00098 lambdaParticle = 0;
00099 baseParticle = 0;
00100 currentLambda = 0;
00101 currentModel = 0;
00102 currentProcess = 0;
00103 loweModel = 0;
00104 chargeSquare = 1.0;
00105 massRatio = 1.0;
00106 mass = 0.0;
00107 currentCut = 0.0;
00108 currentParticleName= "";
00109 currentMaterialName= "";
00110 currentName = "";
00111 lambdaName = "";
00112 theGenericIon = G4GenericIon::GenericIon();
00113 ionEffCharge = new G4ionEffectiveCharge();
00114 isIon = false;
00115 isApplicable = false;
00116 }
00117
00118
00119
00120 G4EmCalculator::~G4EmCalculator()
00121 {
00122 delete ionEffCharge;
00123 for (G4int i=0; i<nLocalMaterials; ++i) {
00124 delete localCouples[i];
00125 }
00126 }
00127
00128
00129
00130 G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4ParticleDefinition* p,
00131 const G4Material* mat, const G4Region* region)
00132 {
00133 G4double res = 0.0;
00134 const G4MaterialCutsCouple* couple = FindCouple(mat, region);
00135 if(couple && UpdateParticle(p, kinEnergy) ) {
00136 res = manager->GetDEDX(p, kinEnergy, couple);
00137
00138 if(isIon) {
00139 if(FindEmModel(p, currentProcessName, kinEnergy)) {
00140 G4double length = CLHEP::nm;
00141 G4double eloss = res*length;
00142
00143
00144 G4double niel = 0.0;
00145 dynParticle.SetKineticEnergy(kinEnergy);
00146 currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
00147 currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
00148 res = eloss/length;
00149
00150
00151 }
00152 }
00153
00154 if(verbose>0) {
00155 G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV
00156 << " DEDX(MeV/mm)= " << res*mm/MeV
00157 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
00158 << " " << p->GetParticleName()
00159 << " in " << mat->GetName()
00160 << " isIon= " << isIon
00161 << G4endl;
00162 }
00163 }
00164 return res;
00165 }
00166
00167
00168
00169 G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4String& particle,
00170 const G4String& material, const G4String& reg)
00171 {
00172 return GetDEDX(kinEnergy,FindParticle(particle),
00173 FindMaterial(material),FindRegion(reg));
00174 }
00175
00176
00177
00178 G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy,
00179 const G4ParticleDefinition* p,
00180 const G4Material* mat,
00181 const G4Region* region)
00182 {
00183 G4double res = 0.0;
00184 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
00185 if(couple && UpdateParticle(p, kinEnergy)) {
00186 res = manager->GetRangeFromRestricteDEDX(p, kinEnergy, couple);
00187 if(verbose>0) {
00188 G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
00189 << " range(mm)= " << res/mm
00190 << " " << p->GetParticleName()
00191 << " in " << mat->GetName()
00192 << G4endl;
00193 }
00194 }
00195 return res;
00196 }
00197
00198
00199
00200 G4double G4EmCalculator::GetCSDARange(G4double kinEnergy,
00201 const G4ParticleDefinition* p,
00202 const G4Material* mat,
00203 const G4Region* region)
00204 {
00205 G4double res = 0.0;
00206 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
00207 if(couple && UpdateParticle(p, kinEnergy)) {
00208 res = manager->GetCSDARange(p, kinEnergy, couple);
00209 if(verbose>0) {
00210 G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
00211 << " range(mm)= " << res/mm
00212 << " " << p->GetParticleName()
00213 << " in " << mat->GetName()
00214 << G4endl;
00215 }
00216 }
00217 return res;
00218 }
00219
00220
00221
00222 G4double G4EmCalculator::GetRange(G4double kinEnergy,
00223 const G4ParticleDefinition* p,
00224 const G4Material* mat,
00225 const G4Region* region)
00226 {
00227 G4double res = 0.0;
00228 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
00229 if(couple && UpdateParticle(p, kinEnergy)) {
00230 res = manager->GetRange(p, kinEnergy, couple);
00231 if(verbose>0) {
00232 G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
00233 << " range(mm)= " << res/mm
00234 << " " << p->GetParticleName()
00235 << " in " << mat->GetName()
00236 << G4endl;
00237 }
00238 }
00239 return res;
00240 }
00241
00242
00243
00244 G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy,
00245 const G4String& particle,
00246 const G4String& material,
00247 const G4String& reg)
00248 {
00249 return GetRangeFromRestricteDEDX(kinEnergy,FindParticle(particle),
00250 FindMaterial(material),FindRegion(reg));
00251 }
00252
00253
00254
00255 G4double G4EmCalculator::GetCSDARange(G4double kinEnergy,
00256 const G4String& particle,
00257 const G4String& material,
00258 const G4String& reg)
00259 {
00260 return GetCSDARange(kinEnergy,FindParticle(particle),
00261 FindMaterial(material),FindRegion(reg));
00262 }
00263
00264
00265
00266 G4double G4EmCalculator::GetRange(G4double kinEnergy,
00267 const G4String& particle,
00268 const G4String& material,
00269 const G4String& reg)
00270 {
00271 return GetRange(kinEnergy,FindParticle(particle),
00272 FindMaterial(material),FindRegion(reg));
00273 }
00274
00275
00276
00277 G4double G4EmCalculator::GetKinEnergy(G4double range,
00278 const G4ParticleDefinition* p,
00279 const G4Material* mat,
00280 const G4Region* region)
00281 {
00282 G4double res = 0.0;
00283 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
00284 if(couple && UpdateParticle(p, 1.0*GeV)) {
00285 res = manager->GetEnergy(p, range, couple);
00286 if(verbose>0) {
00287 G4cout << "G4EmCalculator::GetKinEnergy: Range(mm)= " << range/mm
00288 << " KinE(MeV)= " << res/MeV
00289 << " " << p->GetParticleName()
00290 << " in " << mat->GetName()
00291 << G4endl;
00292 }
00293 }
00294 return res;
00295 }
00296
00297
00298
00299 G4double G4EmCalculator::GetKinEnergy(G4double range, const G4String& particle,
00300 const G4String& material, const G4String& reg)
00301 {
00302 return GetKinEnergy(range,FindParticle(particle),
00303 FindMaterial(material),FindRegion(reg));
00304 }
00305
00306
00307
00308 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
00309 const G4ParticleDefinition* p,
00310 const G4String& processName,
00311 const G4Material* mat,
00312 const G4Region* region)
00313 {
00314 G4double res = 0.0;
00315 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
00316
00317 if(couple && UpdateParticle(p, kinEnergy)) {
00318 G4int idx = couple->GetIndex();
00319 FindLambdaTable(p, processName, kinEnergy);
00320
00321 if(currentLambda) {
00322 G4double e = kinEnergy*massRatio;
00323 res = (((*currentLambda)[idx])->Value(e))*chargeSquare;
00324 } else {
00325 res = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat,
00326 kinEnergy);
00327 }
00328 if(verbose>0) {
00329 G4cout << "G4EmCalculator::GetXSPerVolume: E(MeV)= " << kinEnergy/MeV
00330 << " cross(cm-1)= " << res*cm
00331 << " " << p->GetParticleName()
00332 << " in " << mat->GetName();
00333 if(verbose>1)
00334 G4cout << " idx= " << idx << " Escaled((MeV)= "
00335 << kinEnergy*massRatio
00336 << " q2= " << chargeSquare;
00337 G4cout << G4endl;
00338 }
00339 }
00340 return res;
00341 }
00342
00343
00344
00345 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
00346 const G4String& particle,
00347 const G4String& processName,
00348 const G4String& material,
00349 const G4String& reg)
00350 {
00351 return GetCrossSectionPerVolume(kinEnergy,FindParticle(particle),processName,
00352 FindMaterial(material),FindRegion(reg));
00353 }
00354
00355
00356
00357 G4double G4EmCalculator::GetShellIonisationCrossSectionPerAtom(
00358 const G4String& particle,
00359 G4int Z,
00360 G4AtomicShellEnumerator shell,
00361 G4double kinEnergy)
00362 {
00363 G4double res = 0.0;
00364 const G4ParticleDefinition* p = FindParticle(particle);
00365 G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
00366 if(p && ad) {
00367 res = ad->GetShellIonisationCrossSectionPerAtom(p, Z, shell, kinEnergy);
00368 }
00369 return res;
00370 }
00371
00372
00373
00374 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
00375 const G4ParticleDefinition* p,
00376 const G4String& processName,
00377 const G4Material* mat,
00378 const G4Region* region)
00379 {
00380 G4double res = DBL_MAX;
00381 G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region);
00382 if(x > 0.0) { res = 1.0/x; }
00383 if(verbose>1) {
00384 G4cout << "G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/MeV
00385 << " MFP(mm)= " << res/mm
00386 << " " << p->GetParticleName()
00387 << " in " << mat->GetName()
00388 << G4endl;
00389 }
00390 return res;
00391 }
00392
00393
00394
00395 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
00396 const G4String& particle,
00397 const G4String& processName,
00398 const G4String& material,
00399 const G4String& reg)
00400 {
00401 return GetMeanFreePath(kinEnergy,FindParticle(particle),processName,
00402 FindMaterial(material),FindRegion(reg));
00403 }
00404
00405
00406
00407 void G4EmCalculator::PrintDEDXTable(const G4ParticleDefinition* p)
00408 {
00409 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
00410 G4cout << "##### DEDX Table for " << p->GetParticleName() << G4endl;
00411 if(elp) G4cout << *(elp->DEDXTable()) << G4endl;
00412 }
00413
00414
00415
00416 void G4EmCalculator::PrintRangeTable(const G4ParticleDefinition* p)
00417 {
00418 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
00419 G4cout << "##### Range Table for " << p->GetParticleName() << G4endl;
00420 if(elp) G4cout << *(elp->RangeTableForLoss()) << G4endl;
00421 }
00422
00423
00424
00425 void G4EmCalculator::PrintInverseRangeTable(const G4ParticleDefinition* p)
00426 {
00427 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
00428 G4cout << "### G4EmCalculator: Inverse Range Table for "
00429 << p->GetParticleName() << G4endl;
00430 if(elp) G4cout << *(elp->InverseRangeTable()) << G4endl;
00431 }
00432
00433
00434
00435 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
00436 const G4ParticleDefinition* p,
00437 const G4String& processName,
00438 const G4Material* mat,
00439 G4double cut)
00440 {
00441 currentMaterial = mat;
00442 currentMaterialName = mat->GetName();
00443 G4double res = 0.0;
00444 if(verbose > 1) {
00445 G4cout << "### G4EmCalculator::ComputeDEDX: " << p->GetParticleName()
00446 << " in " << currentMaterialName
00447 << " e(MeV)= " << kinEnergy/MeV << " cut(MeV)= " << cut/MeV
00448 << G4endl;
00449 }
00450 if(UpdateParticle(p, kinEnergy)) {
00451 if(FindEmModel(p, processName, kinEnergy)) {
00452 G4double escaled = kinEnergy*massRatio;
00453 if(baseParticle) {
00454 res = currentModel->ComputeDEDXPerVolume(
00455 mat, baseParticle, escaled, cut) * chargeSquare;
00456 if(verbose > 1) {
00457 G4cout << baseParticle->GetParticleName()
00458 << " Escaled(MeV)= " << escaled;
00459 }
00460 } else {
00461 res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
00462 if(verbose > 1) { G4cout << " no basePart E(MeV)= " << kinEnergy << " "; }
00463 }
00464 if(verbose > 1) {
00465 G4cout << currentModel->GetName() << ": DEDX(MeV/mm)= " << res*mm/MeV
00466 << " DEDX(MeV*cm^2/g)= "
00467 << res*gram/(MeV*cm2*mat->GetDensity())
00468 << G4endl;
00469 }
00470
00471
00472 G4double eth = currentModel->LowEnergyLimit();
00473
00474 if(loweModel) {
00475 G4double res0 = 0.0;
00476 G4double res1 = 0.0;
00477 if(baseParticle) {
00478 res1 = currentModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
00479 * chargeSquare;
00480 res0 = loweModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
00481 * chargeSquare;
00482 } else {
00483 res1 = currentModel->ComputeDEDXPerVolume(mat, p, eth, cut);
00484 res0 = loweModel->ComputeDEDXPerVolume(mat, p, eth, cut);
00485 }
00486 if(verbose > 1) {
00487 G4cout << "At boundary energy(MeV)= " << eth/MeV
00488 << " DEDX(MeV/mm)= " << res1*mm/MeV
00489 << G4endl;
00490 }
00491
00492
00493
00494
00495
00496 if(res1 > 0.0 && escaled > 0.0) {
00497 res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
00498 }
00499 }
00500
00501
00502 if(isIon) {
00503 G4double length = CLHEP::nm;
00504 const G4Region* r = 0;
00505 const G4MaterialCutsCouple* couple = FindCouple(mat, r);
00506 G4double eloss = res*length;
00507 G4double niel = 0.0;
00508 dynParticle.SetKineticEnergy(kinEnergy);
00509 currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
00510 currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
00511 res = eloss/length;
00512
00513 if(verbose > 1) {
00514 G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
00515 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
00516 << G4endl;
00517 }
00518 }
00519 }
00520
00521 if(verbose > 0) {
00522 G4cout << "Sum: E(MeV)= " << kinEnergy/MeV
00523 << " DEDX(MeV/mm)= " << res*mm/MeV
00524 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
00525 << " cut(MeV)= " << cut/MeV
00526 << " " << p->GetParticleName()
00527 << " in " << currentMaterialName
00528 << " Zi^2= " << chargeSquare
00529 << " isIon=" << isIon
00530 << G4endl;
00531 }
00532 }
00533 return res;
00534 }
00535
00536
00537
00538 G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy,
00539 const G4ParticleDefinition* part,
00540 const G4Material* mat,
00541 G4double cut)
00542 {
00543 currentMaterial = mat;
00544 currentMaterialName = mat->GetName();
00545 G4double dedx = 0.0;
00546 if(UpdateParticle(part, kinEnergy)) {
00547
00548 G4LossTableManager* lManager = G4LossTableManager::Instance();
00549 const std::vector<G4VEnergyLossProcess*> vel =
00550 lManager->GetEnergyLossProcessVector();
00551 G4int n = vel.size();
00552
00553
00554
00555
00556 for(G4int i=0; i<n; ++i) {
00557 if(vel[i]) {
00558 G4VProcess* p = reinterpret_cast<G4VProcess*>(vel[i]);
00559 if(ActiveForParticle(part, p)) {
00560
00561
00562 dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),mat,cut);
00563 }
00564 }
00565 }
00566 }
00567 return dedx;
00568 }
00569
00570
00571
00572 G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy, const G4String& part,
00573 const G4String& mat, G4double cut)
00574 {
00575 return ComputeElectronicDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut);
00576 }
00577
00578
00579
00580 G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy,
00581 const G4ParticleDefinition* part,
00582 const G4Material* mat,
00583 G4double cut)
00584 {
00585 G4double dedx = ComputeElectronicDEDX(kinEnergy,part,mat,cut);
00586 if(mass > 700.*MeV) { dedx += ComputeNuclearDEDX(kinEnergy,part,mat); }
00587 return dedx;
00588 }
00589
00590
00591
00592 G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy,
00593 const G4String& part,
00594 const G4String& mat,
00595 G4double cut)
00596 {
00597 return ComputeTotalDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut);
00598 }
00599
00600
00601
00602 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
00603 const G4String& particle,
00604 const G4String& processName,
00605 const G4String& material,
00606 G4double cut)
00607 {
00608 return ComputeDEDX(kinEnergy,FindParticle(particle),processName,
00609 FindMaterial(material),cut);
00610 }
00611
00612
00613
00614 G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
00615 const G4ParticleDefinition* p,
00616 const G4Material* mat)
00617 {
00618
00619 G4double res = corr->NuclearDEDX(p, mat, kinEnergy, false);
00620
00621 if(verbose > 1) {
00622 G4cout << p->GetParticleName() << " E(MeV)= " << kinEnergy/MeV
00623 << " NuclearDEDX(MeV/mm)= " << res*mm/MeV
00624 << " NuclearDEDX(MeV*cm^2/g)= "
00625 << res*gram/(MeV*cm2*mat->GetDensity())
00626 << G4endl;
00627 }
00628 return res;
00629 }
00630
00631
00632
00633 G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
00634 const G4String& particle,
00635 const G4String& material)
00636 {
00637 return ComputeNuclearDEDX(kinEnergy,FindParticle(particle),
00638 FindMaterial(material));
00639 }
00640
00641
00642
00643 G4double G4EmCalculator::ComputeCrossSectionPerVolume(
00644 G4double kinEnergy,
00645 const G4ParticleDefinition* p,
00646 const G4String& processName,
00647 const G4Material* mat,
00648 G4double cut)
00649 {
00650 currentMaterial = mat;
00651 currentMaterialName = mat->GetName();
00652 G4double res = 0.0;
00653 if(UpdateParticle(p, kinEnergy)) {
00654 if(FindEmModel(p, processName, kinEnergy)) {
00655 G4double e = kinEnergy;
00656 if(baseParticle) {
00657 e *= kinEnergy*massRatio;
00658 res = currentModel->CrossSectionPerVolume(
00659 mat, baseParticle, e, cut, e) * chargeSquare;
00660 } else {
00661 res = currentModel->CrossSectionPerVolume(mat, p, e, cut, e);
00662 }
00663 if(verbose>0) {
00664 G4cout << "G4EmCalculator::ComputeXSPerVolume: E(MeV)= " << kinEnergy/MeV
00665 << " cross(cm-1)= " << res*cm
00666 << " cut(keV)= " << cut/keV
00667 << " " << p->GetParticleName()
00668 << " in " << mat->GetName()
00669 << G4endl;
00670 }
00671 }
00672 }
00673 return res;
00674 }
00675
00676
00677
00678 G4double G4EmCalculator::ComputeCrossSectionPerVolume(
00679 G4double kinEnergy,
00680 const G4String& particle,
00681 const G4String& processName,
00682 const G4String& material,
00683 G4double cut)
00684 {
00685 return ComputeCrossSectionPerVolume(kinEnergy,FindParticle(particle),
00686 processName,
00687 FindMaterial(material),cut);
00688 }
00689
00690
00691
00692 G4double G4EmCalculator::ComputeCrossSectionPerAtom(
00693 G4double kinEnergy,
00694 const G4ParticleDefinition* p,
00695 const G4String& processName,
00696 G4double Z, G4double A,
00697 G4double cut)
00698 {
00699 G4double res = 0.0;
00700 if(UpdateParticle(p, kinEnergy)) {
00701 if(FindEmModel(p, processName, kinEnergy)) {
00702 G4double e = kinEnergy;
00703 if(baseParticle) {
00704 e *= kinEnergy*massRatio;
00705 res = currentModel->ComputeCrossSectionPerAtom(
00706 baseParticle, e, Z, A, cut) * chargeSquare;
00707 } else {
00708 res = currentModel->ComputeCrossSectionPerAtom(p, e, Z, A, cut);
00709 }
00710 if(verbose>0) {
00711 G4cout << "E(MeV)= " << kinEnergy/MeV
00712 << " cross(barn)= " << res/barn
00713 << " " << p->GetParticleName()
00714 << " Z= " << Z << " A= " << A/(g/mole) << " g/mole"
00715 << G4endl;
00716 }
00717 }
00718 }
00719 return res;
00720 }
00721
00722
00723
00724 G4double G4EmCalculator::ComputeCrossSectionPerAtom(G4double kinEnergy,
00725 const G4String& particle,
00726 const G4String& processName,
00727 const G4Element* elm,
00728 G4double cut)
00729 {
00730 return ComputeCrossSectionPerAtom(kinEnergy,FindParticle(particle),
00731 processName,
00732 elm->GetZ(),elm->GetN(),cut);
00733 }
00734
00735
00736
00737 G4double
00738 G4EmCalculator::ComputeGammaAttenuationLength(G4double kinEnergy,
00739 const G4Material* mat)
00740 {
00741 G4double res = 0.0;
00742 const G4ParticleDefinition* gamma = G4Gamma::Gamma();
00743 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "conv", mat, 0.0);
00744 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "compt", mat, 0.0);
00745 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "phot", mat, 0.0);
00746 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "Rayl", mat, 0.0);
00747 if(res > 0.0) { res = 1.0/res; }
00748 return res;
00749 }
00750
00751
00752
00753 G4double G4EmCalculator::ComputeShellIonisationCrossSectionPerAtom(
00754 const G4String& particle,
00755 G4int Z,
00756 G4AtomicShellEnumerator shell,
00757 G4double kinEnergy,
00758 const G4Material* mat)
00759 {
00760 G4double res = 0.0;
00761 const G4ParticleDefinition* p = FindParticle(particle);
00762 G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
00763 if(p && ad) {
00764 res = ad->ComputeShellIonisationCrossSectionPerAtom(p, Z, shell,
00765 kinEnergy, mat);
00766 }
00767 return res;
00768 }
00769
00770
00771
00772 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
00773 const G4ParticleDefinition* p,
00774 const G4String& processName,
00775 const G4Material* mat,
00776 G4double cut)
00777 {
00778 G4double mfp = DBL_MAX;
00779 G4double x = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, cut);
00780 if(x > 0.0) { mfp = 1.0/x; }
00781 if(verbose>1) {
00782 G4cout << "E(MeV)= " << kinEnergy/MeV
00783 << " MFP(mm)= " << mfp/mm
00784 << " " << p->GetParticleName()
00785 << " in " << mat->GetName()
00786 << G4endl;
00787 }
00788 return mfp;
00789 }
00790
00791
00792
00793 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
00794 const G4String& particle,
00795 const G4String& processName,
00796 const G4String& material,
00797 G4double cut)
00798 {
00799 return ComputeMeanFreePath(kinEnergy,FindParticle(particle),processName,
00800 FindMaterial(material),cut);
00801 }
00802
00803
00804
00805 G4double G4EmCalculator::ComputeEnergyCutFromRangeCut(
00806 G4double range,
00807 const G4ParticleDefinition* part,
00808 const G4Material* mat)
00809 {
00810 return G4ProductionCutsTable::GetProductionCutsTable()->
00811 ConvertRangeToEnergy(part, mat, range);
00812 }
00813
00814
00815
00816 G4double G4EmCalculator::ComputeEnergyCutFromRangeCut(
00817 G4double range,
00818 const G4String& particle,
00819 const G4String& material)
00820 {
00821 return ComputeEnergyCutFromRangeCut(range,FindParticle(particle),
00822 FindMaterial(material));
00823 }
00824
00825
00826
00827
00828 G4bool G4EmCalculator::UpdateParticle(const G4ParticleDefinition* p,
00829 G4double kinEnergy)
00830 {
00831 if(p != currentParticle) {
00832
00833
00834 currentParticle = p;
00835 dynParticle.SetDefinition(const_cast<G4ParticleDefinition*>(p));
00836 dynParticle.SetKineticEnergy(kinEnergy);
00837 baseParticle = 0;
00838 currentParticleName = p->GetParticleName();
00839 massRatio = 1.0;
00840 mass = p->GetPDGMass();
00841 chargeSquare = 1.0;
00842 currentProcess = FindEnergyLossProcess(p);
00843 currentProcessName = "";
00844 isIon = false;
00845
00846
00847 if(currentProcess) {
00848 currentProcessName = currentProcess->GetProcessName();
00849 baseParticle = currentProcess->BaseParticle();
00850
00851
00852 if(baseParticle) {
00853 massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
00854 G4double q = p->GetPDGCharge()/baseParticle->GetPDGCharge();
00855 chargeSquare = q*q;
00856 }
00857
00858 if(p->GetParticleType() == "nucleus"
00859 && currentParticleName != "deuteron"
00860 && currentParticleName != "triton"
00861 && currentParticleName != "alpha+"
00862 && currentParticleName != "helium"
00863 && currentParticleName != "hydrogen"
00864 ) {
00865 isIon = true;
00866 massRatio = theGenericIon->GetPDGMass()/p->GetPDGMass();
00867 baseParticle = theGenericIon;
00868
00869
00870
00871 }
00872 }
00873 }
00874
00875
00876 if(isIon) {
00877 chargeSquare =
00878 corr->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy)
00879 * corr->EffectiveChargeCorrection(p,currentMaterial,kinEnergy);
00880 if(currentProcess) {
00881 currentProcess->SetDynamicMassCharge(massRatio,chargeSquare);
00882
00883 }
00884 }
00885 return true;
00886 }
00887
00888
00889
00890 const G4ParticleDefinition* G4EmCalculator::FindParticle(const G4String& name)
00891 {
00892 const G4ParticleDefinition* p = 0;
00893 if(name != currentParticleName) {
00894 p = G4ParticleTable::GetParticleTable()->FindParticle(name);
00895 if(!p) {
00896 G4cout << "### WARNING: G4EmCalculator::FindParticle fails to find "
00897 << name << G4endl;
00898 }
00899 } else {
00900 p = currentParticle;
00901 }
00902 return p;
00903 }
00904
00905
00906
00907 const G4ParticleDefinition* G4EmCalculator::FindIon(G4int Z, G4int A)
00908 {
00909 const G4ParticleDefinition* p =
00910 G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
00911 return p;
00912 }
00913
00914
00915
00916 const G4Material* G4EmCalculator::FindMaterial(const G4String& name)
00917 {
00918 if(name != currentMaterialName) {
00919 currentMaterial = G4Material::GetMaterial(name);
00920 currentMaterialName = name;
00921 if(!currentMaterial)
00922 G4cout << "### WARNING: G4EmCalculator::FindMaterial fails to find "
00923 << name << G4endl;
00924 }
00925 return currentMaterial;
00926 }
00927
00928
00929
00930 const G4Region* G4EmCalculator::FindRegion(const G4String& reg)
00931 {
00932 const G4Region* r = 0;
00933 if(reg != "" && reg != "world") {
00934 r = G4RegionStore::GetInstance()->GetRegion(reg);
00935 } else {
00936 r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
00937 }
00938 return r;
00939 }
00940
00941
00942
00943 const G4MaterialCutsCouple* G4EmCalculator::FindCouple(
00944 const G4Material* material,
00945 const G4Region* region)
00946 {
00947 if(!material) return 0;
00948 currentMaterial = material;
00949 currentMaterialName = material->GetName();
00950
00951 const G4ProductionCutsTable* theCoupleTable=
00952 G4ProductionCutsTable::GetProductionCutsTable();
00953 const G4Region* r = region;
00954 if(!r) {
00955 r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
00956 }
00957 return theCoupleTable->GetMaterialCutsCouple(material,r->GetProductionCuts());
00958
00959 }
00960
00961
00962
00963 G4bool G4EmCalculator::UpdateCouple(const G4Material* material, G4double cut)
00964 {
00965 if(!material) return false;
00966 currentMaterial = material;
00967 currentMaterialName = material->GetName();
00968 for (G4int i=0; i<nLocalMaterials; ++i) {
00969 if(material == localMaterials[i] && cut == localCuts[i]) {
00970 currentCouple = localCouples[i];
00971 currentCoupleIndex = currentCouple->GetIndex();
00972 currentCut = cut;
00973 return true;
00974 }
00975 }
00976 const G4MaterialCutsCouple* cc = new G4MaterialCutsCouple(material);
00977 localMaterials.push_back(material);
00978 localCouples.push_back(cc);
00979 localCuts.push_back(cut);
00980 nLocalMaterials++;
00981 currentCouple = cc;
00982 currentCoupleIndex = currentCouple->GetIndex();
00983 currentCut = cut;
00984 return true;
00985 }
00986
00987
00988
00989 void G4EmCalculator::FindLambdaTable(const G4ParticleDefinition* p,
00990 const G4String& processName,
00991 G4double kinEnergy)
00992 {
00993
00994 if (!currentLambda || p != lambdaParticle || processName != lambdaName) {
00995 lambdaName = processName;
00996 currentLambda = 0;
00997 lambdaParticle = p;
00998
00999 const G4ParticleDefinition* part = p;
01000 if(isIon) { part = theGenericIon; }
01001
01002
01003 currentName = processName;
01004 currentModel = 0;
01005 loweModel = 0;
01006
01007 G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
01008 if(elproc) {
01009 currentLambda = elproc->LambdaTable();
01010 if(currentLambda) {
01011 isApplicable = true;
01012 if(verbose>1) {
01013 G4cout << "G4VEnergyLossProcess is found out: " << currentName
01014 << G4endl;
01015 }
01016 }
01017 return;
01018 }
01019
01020
01021 G4VEmProcess* proc = FindDiscreteProcess(part, processName);
01022 if(proc) {
01023 currentLambda = proc->LambdaTable();
01024 if(currentLambda) {
01025 isApplicable = true;
01026 if(verbose>1) {
01027 G4cout << "G4VEmProcess is found out: " << currentName << G4endl;
01028 }
01029 }
01030 return;
01031 }
01032
01033
01034 G4VMultipleScattering* msc = FindMscProcess(part, processName);
01035 if(msc) {
01036 currentModel = msc->SelectModel(kinEnergy,0);
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049 }
01050 }
01051 }
01052
01053
01054
01055 G4bool G4EmCalculator::FindEmModel(const G4ParticleDefinition* p,
01056 const G4String& processName,
01057 G4double kinEnergy)
01058 {
01059 isApplicable = false;
01060 if(!p) {
01061 G4cout << "G4EmCalculator::FindEmModel WARNING: no particle defined"
01062 << G4endl;
01063 return isApplicable;
01064 }
01065 G4String partname = p->GetParticleName();
01066 const G4ParticleDefinition* part = p;
01067 G4double scaledEnergy = kinEnergy*massRatio;
01068 if(isIon) { part = theGenericIon; }
01069
01070 if(verbose > 1) {
01071 G4cout << "## G4EmCalculator::FindEmModel for " << partname
01072 << " (type= " << p->GetParticleType()
01073 << ") and " << processName << " at E(MeV)= " << scaledEnergy
01074 << G4endl;
01075 if(p != part) { G4cout << " GenericIon is the base particle" << G4endl; }
01076 }
01077
01078
01079 currentName = processName;
01080 currentModel = 0;
01081 loweModel = 0;
01082 size_t idx = 0;
01083
01084 G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
01085 if(elproc) {
01086 currentModel = elproc->SelectModelForMaterial(scaledEnergy, idx);
01087 G4double eth = currentModel->LowEnergyLimit();
01088 if(eth > 0.0) {
01089 loweModel = elproc->SelectModelForMaterial(eth - CLHEP::eV, idx);
01090 if(loweModel == currentModel) { loweModel = 0; }
01091 }
01092 }
01093
01094
01095 if(!currentModel) {
01096 G4VEmProcess* proc = FindDiscreteProcess(part, processName);
01097 if(proc) {
01098 currentModel = proc->SelectModelForMaterial(kinEnergy, idx);
01099 G4double eth = currentModel->LowEnergyLimit();
01100 if(eth > 0.0) {
01101 loweModel = proc->SelectModelForMaterial(eth - CLHEP::eV, idx);
01102 if(loweModel == currentModel) { loweModel = 0; }
01103 }
01104 }
01105 }
01106
01107
01108 if(!currentModel) {
01109 G4VMultipleScattering* proc = FindMscProcess(part, processName);
01110 if(proc) {
01111 currentModel = proc->SelectModel(kinEnergy, idx);
01112 loweModel = 0;
01113 }
01114 }
01115 if(currentModel) {
01116 if(loweModel == currentModel) { loweModel = 0; }
01117 isApplicable = true;
01118 if(verbose > 1) {
01119 G4cout << " Model <" << currentModel->GetName()
01120 << "> Emin(MeV)= " << currentModel->LowEnergyLimit()/MeV
01121 << " for " << part->GetParticleName();
01122 if(elproc) {
01123 G4cout << " and " << elproc->GetProcessName() << " " << elproc
01124 << G4endl;
01125 }
01126 if(loweModel) {
01127 G4cout << " LowEnergy model <" << loweModel->GetName() << ">";
01128 }
01129 G4cout << G4endl;
01130 }
01131 }
01132 return isApplicable;
01133 }
01134
01135
01136
01137 G4VEnergyLossProcess* G4EmCalculator::FindEnergyLossProcess(
01138 const G4ParticleDefinition* p)
01139 {
01140 G4VEnergyLossProcess* elp = 0;
01141 G4String partname = p->GetParticleName();
01142 const G4ParticleDefinition* part = p;
01143
01144 if(p->GetParticleType() == "nucleus"
01145 && currentParticleName != "deuteron"
01146 && currentParticleName != "triton"
01147 && currentParticleName != "alpha+"
01148 && currentParticleName != "helium"
01149 && currentParticleName != "hydrogen"
01150 ) { part = theGenericIon; }
01151
01152 elp = G4LossTableManager::Instance()->GetEnergyLossProcess(part);
01153 return elp;
01154 }
01155
01156
01157
01158 G4VEnergyLossProcess*
01159 G4EmCalculator::FindEnLossProcess(const G4ParticleDefinition* part,
01160 const G4String& processName)
01161 {
01162 G4VEnergyLossProcess* proc = 0;
01163 const std::vector<G4VEnergyLossProcess*> v =
01164 G4LossTableManager::Instance()->GetEnergyLossProcessVector();
01165 G4int n = v.size();
01166 for(G4int i=0; i<n; ++i) {
01167 if((v[i])->GetProcessName() == processName) {
01168 G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
01169 if(ActiveForParticle(part, p)) {
01170 proc = v[i];
01171 break;
01172 }
01173 }
01174 }
01175 return proc;
01176 }
01177
01178
01179
01180 G4VEmProcess*
01181 G4EmCalculator::FindDiscreteProcess(const G4ParticleDefinition* part,
01182 const G4String& processName)
01183 {
01184 G4VEmProcess* proc = 0;
01185 const std::vector<G4VEmProcess*> v =
01186 G4LossTableManager::Instance()->GetEmProcessVector();
01187 G4int n = v.size();
01188 for(G4int i=0; i<n; ++i) {
01189 if((v[i])->GetProcessName() == processName) {
01190 G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
01191 if(ActiveForParticle(part, p)) {
01192 proc = v[i];
01193 break;
01194 }
01195 }
01196 }
01197 return proc;
01198 }
01199
01200
01201
01202 G4VMultipleScattering*
01203 G4EmCalculator::FindMscProcess(const G4ParticleDefinition* part,
01204 const G4String& processName)
01205 {
01206 G4VMultipleScattering* proc = 0;
01207 const std::vector<G4VMultipleScattering*> v =
01208 G4LossTableManager::Instance()->GetMultipleScatteringVector();
01209 G4int n = v.size();
01210 for(G4int i=0; i<n; ++i) {
01211 if((v[i])->GetProcessName() == processName) {
01212 G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
01213 if(ActiveForParticle(part, p)) {
01214 proc = v[i];
01215 break;
01216 }
01217 }
01218 }
01219 return proc;
01220 }
01221
01222
01223
01224 G4bool G4EmCalculator::ActiveForParticle(const G4ParticleDefinition* part,
01225 G4VProcess* proc)
01226 {
01227 G4ProcessManager* pm = part->GetProcessManager();
01228 G4ProcessVector* pv = pm->GetProcessList();
01229 G4int n = pv->size();
01230 G4bool res = false;
01231 for(G4int i=0; i<n; ++i) {
01232 if((*pv)[i] == proc) {
01233 if(pm->GetProcessActivation(i)) { res = true; }
01234 break;
01235 }
01236 }
01237 return res;
01238 }
01239
01240
01241
01242 void G4EmCalculator::SetVerbose(G4int verb)
01243 {
01244 verbose = verb;
01245 }
01246
01247
01248