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
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 #include "G4VEnergyLossProcess.hh"
00126 #include "G4PhysicalConstants.hh"
00127 #include "G4SystemOfUnits.hh"
00128 #include "G4ProcessManager.hh"
00129 #include "G4LossTableManager.hh"
00130 #include "G4LossTableBuilder.hh"
00131 #include "G4Step.hh"
00132 #include "G4ParticleDefinition.hh"
00133 #include "G4VEmModel.hh"
00134 #include "G4VEmFluctuationModel.hh"
00135 #include "G4DataVector.hh"
00136 #include "G4PhysicsLogVector.hh"
00137 #include "G4VParticleChange.hh"
00138 #include "G4Gamma.hh"
00139 #include "G4Electron.hh"
00140 #include "G4Positron.hh"
00141 #include "G4Proton.hh"
00142 #include "G4ProcessManager.hh"
00143 #include "G4UnitsTable.hh"
00144 #include "G4GenericIon.hh"
00145 #include "G4ProductionCutsTable.hh"
00146 #include "G4Region.hh"
00147 #include "G4RegionStore.hh"
00148 #include "G4PhysicsTableHelper.hh"
00149 #include "G4SafetyHelper.hh"
00150 #include "G4TransportationManager.hh"
00151 #include "G4EmConfigurator.hh"
00152 #include "G4VAtomDeexcitation.hh"
00153 #include "G4EmBiasingManager.hh"
00154
00155
00156
00157 G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name,
00158 G4ProcessType type):
00159 G4VContinuousDiscreteProcess(name, type),
00160 secondaryParticle(0),
00161 nSCoffRegions(0),
00162 idxSCoffRegions(0),
00163 nProcesses(0),
00164 theDEDXTable(0),
00165 theDEDXSubTable(0),
00166 theDEDXunRestrictedTable(0),
00167 theIonisationTable(0),
00168 theIonisationSubTable(0),
00169 theRangeTableForLoss(0),
00170 theCSDARangeTable(0),
00171 theSecondaryRangeTable(0),
00172 theInverseRangeTable(0),
00173 theLambdaTable(0),
00174 theSubLambdaTable(0),
00175 theDensityFactor(0),
00176 theDensityIdx(0),
00177 baseParticle(0),
00178 minSubRange(0.1),
00179 lossFluctuationFlag(true),
00180 rndmStepFlag(false),
00181 tablesAreBuilt(false),
00182 integral(true),
00183 isIon(false),
00184 isIonisation(true),
00185 useSubCutoff(false),
00186 useDeexcitation(false),
00187 particle(0),
00188 currentCouple(0),
00189 nWarnings(0),
00190 mfpKinEnergy(0.0)
00191 {
00192 SetVerboseLevel(1);
00193
00194
00195 lowestKinEnergy = 1.*eV;
00196
00197
00198 minKinEnergy = 0.1*keV;
00199 maxKinEnergy = 10.0*TeV;
00200 nBins = 77;
00201 maxKinEnergyCSDA = 1.0*GeV;
00202 nBinsCSDA = 35;
00203
00204
00205 linLossLimit = 0.01;
00206
00207
00208 SetStepFunction(0.2, 1.0*mm);
00209
00210
00211 lambdaFactor = 0.8;
00212
00213
00214 biasFactor = 1.0;
00215
00216
00217 theElectron = G4Electron::Electron();
00218 thePositron = G4Positron::Positron();
00219 theGamma = G4Gamma::Gamma();
00220 theGenericIon = 0;
00221
00222
00223 pParticleChange = &fParticleChange;
00224 fParticleChange.SetSecondaryWeightByProcess(true);
00225 modelManager = new G4EmModelManager();
00226 safetyHelper = G4TransportationManager::GetTransportationManager()
00227 ->GetSafetyHelper();
00228 aGPILSelection = CandidateForSelection;
00229
00230
00231 (G4LossTableManager::Instance())->Register(this);
00232 fluctModel = 0;
00233 atomDeexcitation = 0;
00234
00235 biasManager = 0;
00236 biasFlag = false;
00237 weightFlag = false;
00238
00239 scTracks.reserve(5);
00240 secParticles.reserve(5);
00241 }
00242
00243
00244
00245 G4VEnergyLossProcess::~G4VEnergyLossProcess()
00246 {
00247 if(1 < verboseLevel) {
00248 G4cout << "G4VEnergyLossProcess destruct " << GetProcessName()
00249 << " " << this << " " << baseParticle << G4endl;
00250 }
00251 Clean();
00252
00253 if ( !baseParticle ) {
00254 if(theDEDXTable) {
00255 if(theIonisationTable == theDEDXTable) { theIonisationTable = 0; }
00256 theDEDXTable->clearAndDestroy();
00257 delete theDEDXTable;
00258 if(theDEDXSubTable) {
00259 if(theIonisationSubTable == theDEDXSubTable)
00260 theIonisationSubTable = 0;
00261 theDEDXSubTable->clearAndDestroy();
00262 delete theDEDXSubTable;
00263 }
00264 }
00265 if(theIonisationTable) {
00266 theIonisationTable->clearAndDestroy();
00267 delete theIonisationTable;
00268 }
00269 if(theIonisationSubTable) {
00270 theIonisationSubTable->clearAndDestroy();
00271 delete theIonisationSubTable;
00272 }
00273 if(theDEDXunRestrictedTable && theCSDARangeTable) {
00274 theDEDXunRestrictedTable->clearAndDestroy();
00275 delete theDEDXunRestrictedTable;
00276 }
00277 if(theCSDARangeTable) {
00278 theCSDARangeTable->clearAndDestroy();
00279 delete theCSDARangeTable;
00280 }
00281 if(theRangeTableForLoss) {
00282 theRangeTableForLoss->clearAndDestroy();
00283 delete theRangeTableForLoss;
00284 }
00285 if(theInverseRangeTable) {
00286 theInverseRangeTable->clearAndDestroy();
00287 delete theInverseRangeTable;
00288 }
00289 if(theLambdaTable) {
00290 theLambdaTable->clearAndDestroy();
00291 delete theLambdaTable;
00292 }
00293 if(theSubLambdaTable) {
00294 theSubLambdaTable->clearAndDestroy();
00295 delete theSubLambdaTable;
00296 }
00297 }
00298
00299 delete modelManager;
00300 delete biasManager;
00301 (G4LossTableManager::Instance())->DeRegister(this);
00302 }
00303
00304
00305
00306 void G4VEnergyLossProcess::Clean()
00307 {
00308 if(1 < verboseLevel) {
00309 G4cout << "G4VEnergyLossProcess::Clear() for " << GetProcessName()
00310 << G4endl;
00311 }
00312 delete [] idxSCoffRegions;
00313
00314 tablesAreBuilt = false;
00315
00316 scProcesses.clear();
00317 nProcesses = 0;
00318 }
00319
00320
00321
00322 G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
00323 const G4Material*,
00324 G4double cut)
00325 {
00326 return cut;
00327 }
00328
00329
00330
00331 void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* p,
00332 G4VEmFluctuationModel* fluc,
00333 const G4Region* region)
00334 {
00335 modelManager->AddEmModel(order, p, fluc, region);
00336 if(p) { p->SetParticleChange(pParticleChange, fluc); }
00337 }
00338
00339
00340
00341 void G4VEnergyLossProcess::UpdateEmModel(const G4String& nam,
00342 G4double emin, G4double emax)
00343 {
00344 modelManager->UpdateEmModel(nam, emin, emax);
00345 }
00346
00347
00348
00349 void G4VEnergyLossProcess::SetEmModel(G4VEmModel* p, G4int index)
00350 {
00351 G4int n = emModels.size();
00352 if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
00353 emModels[index] = p;
00354 }
00355
00356
00357
00358 G4VEmModel* G4VEnergyLossProcess::EmModel(G4int index)
00359 {
00360 G4VEmModel* p = 0;
00361 if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; }
00362 return p;
00363 }
00364
00365
00366
00367 G4VEmModel* G4VEnergyLossProcess::GetModelByIndex(G4int idx, G4bool ver)
00368 {
00369 return modelManager->GetModel(idx, ver);
00370 }
00371
00372
00373
00374 G4int G4VEnergyLossProcess::NumberOfModels()
00375 {
00376 return modelManager->NumberOfModels();
00377 }
00378
00379
00380
00381 void
00382 G4VEnergyLossProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
00383 {
00384 if(1 < verboseLevel) {
00385 G4cout << "G4VEnergyLossProcess::PreparePhysicsTable for "
00386 << GetProcessName() << " for " << part.GetParticleName()
00387 << " " << this << G4endl;
00388 }
00389
00390 currentCouple = 0;
00391 preStepLambda = 0.0;
00392 mfpKinEnergy = DBL_MAX;
00393 fRange = DBL_MAX;
00394 preStepKinEnergy = 0.0;
00395 chargeSqRatio = 1.0;
00396 massRatio = 1.0;
00397 reduceFactor = 1.0;
00398 fFactor = 1.0;
00399
00400 G4LossTableManager* lManager = G4LossTableManager::Instance();
00401
00402
00403 if( !particle ) { particle = ∂ }
00404
00405 if(part.GetParticleType() == "nucleus") {
00406
00407 G4String pname = part.GetParticleName();
00408 if(pname != "deuteron" && pname != "triton" &&
00409 pname != "alpha+" && pname != "helium" &&
00410 pname != "hydrogen") {
00411
00412 theGenericIon = G4GenericIon::GenericIon();
00413 isIon = true;
00414
00415
00416 if(pname != "He3" && pname != "alpha") { particle = theGenericIon; }
00417 }
00418 }
00419
00420 if( particle != &part ) {
00421 if(isIon) {
00422 lManager->RegisterIon(&part, this);
00423 } else {
00424 lManager->RegisterExtraParticle(&part, this);
00425 }
00426 if(1 < verboseLevel) {
00427 G4cout << "### G4VEnergyLossProcess::PreparePhysicsTable() interrupted for "
00428 << part.GetParticleName() << " isIon= " << isIon
00429 << " particle " << particle << " GenericIon " << theGenericIon
00430 << G4endl;
00431 }
00432 return;
00433 }
00434
00435 Clean();
00436 lManager->PreparePhysicsTable(&part, this);
00437 G4LossTableBuilder* bld = lManager->GetTableBuilder();
00438
00439
00440 InitialiseEnergyLossProcess(particle, baseParticle);
00441
00442 const G4ProductionCutsTable* theCoupleTable=
00443 G4ProductionCutsTable::GetProductionCutsTable();
00444 size_t n = theCoupleTable->GetTableSize();
00445
00446 theDEDXAtMaxEnergy.resize(n, 0.0);
00447 theRangeAtMaxEnergy.resize(n, 0.0);
00448 theEnergyOfCrossSectionMax.resize(n, 0.0);
00449 theCrossSectionMax.resize(n, DBL_MAX);
00450
00451
00452 if (!baseParticle) {
00453
00454 theDEDXTable = G4PhysicsTableHelper::PreparePhysicsTable(theDEDXTable);
00455 bld->InitialiseBaseMaterials(theDEDXTable);
00456
00457 if (lManager->BuildCSDARange()) {
00458 theDEDXunRestrictedTable =
00459 G4PhysicsTableHelper::PreparePhysicsTable(theDEDXunRestrictedTable);
00460 theCSDARangeTable =
00461 G4PhysicsTableHelper::PreparePhysicsTable(theCSDARangeTable);
00462 }
00463
00464 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
00465 bld->InitialiseBaseMaterials(theLambdaTable);
00466
00467 if(isIonisation) {
00468 theRangeTableForLoss =
00469 G4PhysicsTableHelper::PreparePhysicsTable(theRangeTableForLoss);
00470 theInverseRangeTable =
00471 G4PhysicsTableHelper::PreparePhysicsTable(theInverseRangeTable);
00472 }
00473
00474 if (nSCoffRegions) {
00475 theDEDXSubTable =
00476 G4PhysicsTableHelper::PreparePhysicsTable(theDEDXSubTable);
00477 theSubLambdaTable =
00478 G4PhysicsTableHelper::PreparePhysicsTable(theSubLambdaTable);
00479 }
00480 }
00481
00482 theDensityFactor = bld->GetDensityFactors();
00483 theDensityIdx = bld->GetCoupleIndexes();
00484
00485
00486 if(biasManager) {
00487 biasManager->Initialise(part,GetProcessName(),verboseLevel);
00488 biasFlag = false;
00489 }
00490
00491 G4double initialCharge = particle->GetPDGCharge();
00492 G4double initialMass = particle->GetPDGMass();
00493
00494 if (baseParticle) {
00495 massRatio = (baseParticle->GetPDGMass())/initialMass;
00496 G4double q = initialCharge/baseParticle->GetPDGCharge();
00497 chargeSqRatio = q*q;
00498 if(chargeSqRatio > 0.0) { reduceFactor = 1.0/(chargeSqRatio*massRatio); }
00499 }
00500
00501
00502 G4int nmod = modelManager->NumberOfModels();
00503 for(G4int i=0; i<nmod; ++i) {
00504 G4VEmModel* mod = modelManager->GetModel(i);
00505 if(mod->HighEnergyLimit() > maxKinEnergy) {
00506 mod->SetHighEnergyLimit(maxKinEnergy);
00507 }
00508 }
00509
00510 theCuts = modelManager->Initialise(particle, secondaryParticle,
00511 minSubRange, verboseLevel);
00512
00513
00514 if (nSCoffRegions>0) {
00515 theSubCuts = modelManager->SubCutoff();
00516
00517 if(nSCoffRegions>0) { idxSCoffRegions = new G4bool[n]; }
00518 for (size_t j=0; j<n; ++j) {
00519
00520 const G4MaterialCutsCouple* couple =
00521 theCoupleTable->GetMaterialCutsCouple(j);
00522 const G4ProductionCuts* pcuts = couple->GetProductionCuts();
00523
00524 if(nSCoffRegions>0) {
00525 G4bool reg = false;
00526 for(G4int i=0; i<nSCoffRegions; ++i) {
00527 if( pcuts == scoffRegions[i]->GetProductionCuts()) { reg = true; }
00528 }
00529 idxSCoffRegions[j] = reg;
00530 }
00531 }
00532 }
00533
00534 if(1 < verboseLevel) {
00535 G4cout << "G4VEnergyLossProcess::Initialise() is done "
00536 << " for local " << particle->GetParticleName()
00537 << " isIon= " << isIon;
00538 if(baseParticle) { G4cout << "; base: " << baseParticle->GetParticleName(); }
00539 G4cout << " chargeSqRatio= " << chargeSqRatio
00540 << " massRatio= " << massRatio
00541 << " reduceFactor= " << reduceFactor << G4endl;
00542 if (nSCoffRegions) {
00543 G4cout << " SubCutoff Regime is ON for regions: " << G4endl;
00544 for (G4int i=0; i<nSCoffRegions; ++i) {
00545 const G4Region* r = scoffRegions[i];
00546 G4cout << " " << r->GetName() << G4endl;
00547 }
00548 }
00549 }
00550 }
00551
00552
00553
00554 void G4VEnergyLossProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
00555 {
00556 if(1 < verboseLevel) {
00557 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for "
00558 << GetProcessName()
00559 << " and particle " << part.GetParticleName()
00560 << "; local: " << particle->GetParticleName();
00561 if(baseParticle) { G4cout << "; base: " << baseParticle->GetParticleName(); }
00562 G4cout << " TablesAreBuilt= " << tablesAreBuilt
00563 << " isIon= " << isIon << " " << this << G4endl;
00564 }
00565
00566 if(&part == particle) {
00567 if(!tablesAreBuilt) {
00568 G4LossTableManager::Instance()->BuildPhysicsTable(particle, this);
00569 }
00570 if(!baseParticle) {
00571
00572 safetyHelper->InitialiseHelper();
00573 }
00574 }
00575
00576
00577 G4String num = part.GetParticleName();
00578 if(1 < verboseLevel ||
00579 (0 < verboseLevel && (num == "e-" ||
00580 num == "e+" || num == "mu+" ||
00581 num == "mu-" || num == "proton"||
00582 num == "pi+" || num == "pi-" ||
00583 num == "kaon+" || num == "kaon-" ||
00584 num == "alpha" || num == "anti_proton" ||
00585 num == "GenericIon")))
00586 {
00587 particle = ∂
00588 PrintInfoDefinition();
00589 }
00590
00591
00592
00593 if(isIonisation) {
00594 fParticleChange.SetLowEnergyLimit(lowestKinEnergy);
00595 atomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00596 if(atomDeexcitation) {
00597 if(atomDeexcitation->IsPIXEActive()) { useDeexcitation = true; }
00598 }
00599 }
00600
00601 if(1 < verboseLevel) {
00602 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for "
00603 << GetProcessName()
00604 << " and particle " << part.GetParticleName();
00605 if(isIonisation) { G4cout << " isIonisation flag = 1"; }
00606 G4cout << G4endl;
00607 }
00608 }
00609
00610
00611
00612 G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTable(G4EmTableType tType)
00613 {
00614 if(1 < verboseLevel) {
00615 G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType
00616 << " for " << GetProcessName()
00617 << " and particle " << particle->GetParticleName()
00618 << G4endl;
00619 }
00620 G4PhysicsTable* table = 0;
00621 G4double emax = maxKinEnergy;
00622 G4int bin = nBins;
00623
00624 if(fTotal == tType) {
00625 emax = maxKinEnergyCSDA;
00626 bin = nBinsCSDA;
00627 table = theDEDXunRestrictedTable;
00628 } else if(fRestricted == tType) {
00629 table = theDEDXTable;
00630 if(theIonisationTable)
00631 table = G4PhysicsTableHelper::PreparePhysicsTable(theIonisationTable);
00632 } else if(fSubRestricted == tType) {
00633 table = theDEDXSubTable;
00634 if(theIonisationSubTable)
00635 table = G4PhysicsTableHelper::PreparePhysicsTable(theIonisationSubTable);
00636 } else {
00637 G4cout << "G4VEnergyLossProcess::BuildDEDXTable WARNING: wrong type "
00638 << tType << G4endl;
00639 }
00640
00641
00642 const G4ProductionCutsTable* theCoupleTable=
00643 G4ProductionCutsTable::GetProductionCutsTable();
00644 size_t numOfCouples = theCoupleTable->GetTableSize();
00645
00646 if(1 < verboseLevel) {
00647 G4cout << numOfCouples << " materials"
00648 << " minKinEnergy= " << minKinEnergy
00649 << " maxKinEnergy= " << emax
00650 << " nbin= " << bin
00651 << " EmTableType= " << tType
00652 << " table= " << table << " " << this
00653 << G4endl;
00654 }
00655 if(!table) return table;
00656
00657 G4LossTableBuilder* bld = (G4LossTableManager::Instance())->GetTableBuilder();
00658 G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
00659 G4PhysicsLogVector* aVector = 0;
00660 G4PhysicsLogVector* bVector = 0;
00661
00662 for(size_t i=0; i<numOfCouples; ++i) {
00663
00664 if(1 < verboseLevel) {
00665 G4cout << "G4VEnergyLossProcess::BuildDEDXVector flagTable= "
00666 << table->GetFlag(i) << " Flag= " << bld->GetFlag(i) << G4endl;
00667 }
00668 if(bld->GetFlag(i)) {
00669
00670
00671 const G4MaterialCutsCouple* couple =
00672 theCoupleTable->GetMaterialCutsCouple(i);
00673 delete (*table)[i];
00674 if(!bVector) {
00675 aVector = new G4PhysicsLogVector(minKinEnergy, emax, bin);
00676 bVector = aVector;
00677 } else {
00678 aVector = new G4PhysicsLogVector(*bVector);
00679 }
00680 aVector->SetSpline(splineFlag);
00681
00682 modelManager->FillDEDXVector(aVector, couple, tType);
00683 if(splineFlag) { aVector->FillSecondDerivatives(); }
00684
00685
00686 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
00687 }
00688 }
00689
00690 if(1 < verboseLevel) {
00691 G4cout << "G4VEnergyLossProcess::BuildDEDXTable(): table is built for "
00692 << particle->GetParticleName()
00693 << " and process " << GetProcessName()
00694 << G4endl;
00695
00696 }
00697
00698 return table;
00699 }
00700
00701
00702
00703 G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaTable(G4EmTableType tType)
00704 {
00705 G4PhysicsTable* table = 0;
00706
00707 if(fRestricted == tType) {
00708 table = theLambdaTable;
00709 } else if(fSubRestricted == tType) {
00710 table = theSubLambdaTable;
00711 } else {
00712 G4cout << "G4VEnergyLossProcess::BuildLambdaTable WARNING: wrong type "
00713 << tType << G4endl;
00714 }
00715
00716 if(1 < verboseLevel) {
00717 G4cout << "G4VEnergyLossProcess::BuildLambdaTable() of type "
00718 << tType << " for process "
00719 << GetProcessName() << " and particle "
00720 << particle->GetParticleName()
00721 << " EmTableType= " << tType
00722 << " table= " << table
00723 << G4endl;
00724 }
00725 if(!table) {return table;}
00726
00727
00728 const G4ProductionCutsTable* theCoupleTable=
00729 G4ProductionCutsTable::GetProductionCutsTable();
00730 size_t numOfCouples = theCoupleTable->GetTableSize();
00731
00732 G4LossTableBuilder* bld = (G4LossTableManager::Instance())->GetTableBuilder();
00733 G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
00734 G4PhysicsLogVector* aVector = 0;
00735 G4double scale = std::log(maxKinEnergy/minKinEnergy);
00736
00737 for(size_t i=0; i<numOfCouples; ++i) {
00738
00739 if (bld->GetFlag(i)) {
00740
00741
00742 const G4MaterialCutsCouple* couple =
00743 theCoupleTable->GetMaterialCutsCouple(i);
00744 delete (*table)[i];
00745 G4double emin = MinPrimaryEnergy(particle,couple->GetMaterial(),(*theCuts)[i]);
00746 if(0.0 >= emin) { emin = eV; }
00747 else if(maxKinEnergy <= emin) { emin = 0.5*maxKinEnergy; }
00748 G4int bin = G4int(nBins*std::log(maxKinEnergy/emin)/scale + 0.5);
00749 if(bin < 3) { bin = 3; }
00750 aVector = new G4PhysicsLogVector(emin, maxKinEnergy, bin);
00751 aVector->SetSpline(splineFlag);
00752
00753 modelManager->FillLambdaVector(aVector, couple, true, tType);
00754 if(splineFlag) { aVector->FillSecondDerivatives(); }
00755
00756
00757 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
00758 }
00759 }
00760
00761 if(1 < verboseLevel) {
00762 G4cout << "Lambda table is built for "
00763 << particle->GetParticleName()
00764 << G4endl;
00765 }
00766
00767 return table;
00768 }
00769
00770
00771
00772 void G4VEnergyLossProcess::PrintInfoDefinition()
00773 {
00774 if(0 < verboseLevel) {
00775 G4cout << G4endl << GetProcessName() << ": for "
00776 << particle->GetParticleName()
00777 << " SubType= " << GetProcessSubType()
00778 << G4endl
00779 << " dE/dx and range tables from "
00780 << G4BestUnit(minKinEnergy,"Energy")
00781 << " to " << G4BestUnit(maxKinEnergy,"Energy")
00782 << " in " << nBins << " bins" << G4endl
00783 << " Lambda tables from threshold to "
00784 << G4BestUnit(maxKinEnergy,"Energy")
00785 << " in " << nBins << " bins, spline: "
00786 << (G4LossTableManager::Instance())->SplineFlag()
00787 << G4endl;
00788 if(theRangeTableForLoss && isIonisation) {
00789 G4cout << " finalRange(mm)= " << finalRange/mm
00790 << ", dRoverRange= " << dRoverRange
00791 << ", integral: " << integral
00792 << ", fluct: " << lossFluctuationFlag
00793 << ", linLossLimit= " << linLossLimit
00794 << G4endl;
00795 }
00796 PrintInfo();
00797 modelManager->DumpModelList(verboseLevel);
00798 if(theCSDARangeTable && isIonisation) {
00799 G4cout << " CSDA range table up"
00800 << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
00801 << " in " << nBinsCSDA << " bins" << G4endl;
00802 }
00803 if(nSCoffRegions>0 && isIonisation) {
00804 G4cout << " Subcutoff sampling in " << nSCoffRegions
00805 << " regions" << G4endl;
00806 }
00807 if(2 < verboseLevel) {
00808 G4cout << " DEDXTable address= " << theDEDXTable << G4endl;
00809 if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
00810 G4cout << "non restricted DEDXTable address= "
00811 << theDEDXunRestrictedTable << G4endl;
00812 if(theDEDXunRestrictedTable && isIonisation) {
00813 G4cout << (*theDEDXunRestrictedTable) << G4endl;
00814 }
00815 if(theDEDXSubTable && isIonisation) {
00816 G4cout << (*theDEDXSubTable) << G4endl;
00817 }
00818 G4cout << " CSDARangeTable address= " << theCSDARangeTable
00819 << G4endl;
00820 if(theCSDARangeTable && isIonisation) {
00821 G4cout << (*theCSDARangeTable) << G4endl;
00822 }
00823 G4cout << " RangeTableForLoss address= " << theRangeTableForLoss
00824 << G4endl;
00825 if(theRangeTableForLoss && isIonisation) {
00826 G4cout << (*theRangeTableForLoss) << G4endl;
00827 }
00828 G4cout << " InverseRangeTable address= " << theInverseRangeTable
00829 << G4endl;
00830 if(theInverseRangeTable && isIonisation) {
00831 G4cout << (*theInverseRangeTable) << G4endl;
00832 }
00833 G4cout << " LambdaTable address= " << theLambdaTable << G4endl;
00834 if(theLambdaTable && isIonisation) {
00835 G4cout << (*theLambdaTable) << G4endl;
00836 }
00837 G4cout << " SubLambdaTable address= " << theSubLambdaTable << G4endl;
00838 if(theSubLambdaTable && isIonisation) {
00839 G4cout << (*theSubLambdaTable) << G4endl;
00840 }
00841 }
00842 }
00843 }
00844
00845
00846
00847 void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
00848 {
00849 G4RegionStore* regionStore = G4RegionStore::GetInstance();
00850 const G4Region* reg = r;
00851 if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
00852
00853
00854 if (nSCoffRegions) {
00855 for (G4int i=0; i<nSCoffRegions; ++i) {
00856 if (reg == scoffRegions[i]) {
00857 return;
00858 }
00859 }
00860 }
00861
00862
00863 if(val) {
00864 useSubCutoff = true;
00865 scoffRegions.push_back(reg);
00866 ++nSCoffRegions;
00867 } else {
00868 useSubCutoff = false;
00869 }
00870 }
00871
00872
00873
00874 void G4VEnergyLossProcess::StartTracking(G4Track* track)
00875 {
00876
00877 theNumberOfInteractionLengthLeft = -1.0;
00878 mfpKinEnergy = DBL_MAX;
00879
00880
00881 if(isIon) {
00882 chargeSqRatio = 0.5;
00883
00884 G4double newmass = track->GetDefinition()->GetPDGMass();
00885 if(baseParticle) {
00886 massRatio = baseParticle->GetPDGMass()/newmass;
00887 } else {
00888 massRatio = proton_mass_c2/newmass;
00889 }
00890 }
00891
00892 if(biasManager) {
00893 if(0 == track->GetParentID()) {
00894
00895 biasFlag = true;
00896 biasManager->ResetForcedInteraction();
00897 }
00898 }
00899 }
00900
00901
00902
00903 G4double G4VEnergyLossProcess::AlongStepGetPhysicalInteractionLength(
00904 const G4Track&,G4double,G4double,G4double&,
00905 G4GPILSelection* selection)
00906 {
00907 G4double x = DBL_MAX;
00908 *selection = aGPILSelection;
00909 if(isIonisation) {
00910 fRange = GetScaledRangeForScaledEnergy(preStepScaledEnergy)*reduceFactor;
00911
00912 x = fRange;
00913 G4double y = x*dRoverRange;
00914 G4double finR = finalRange;
00915 if(rndmStepFlag) {
00916 finR = std::min(finR,currentCouple->GetProductionCuts()->GetProductionCut(1));
00917 }
00918 if(x > finR) { x = y + finR*(1.0 - dRoverRange)*(2.0 - finR/fRange); }
00919
00920
00921
00922
00923
00924
00925
00926 }
00927
00928 return x;
00929 }
00930
00931
00932
00933 G4double G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength(
00934 const G4Track& track,
00935 G4double previousStepSize,
00936 G4ForceCondition* condition)
00937 {
00938
00939 *condition = NotForced;
00940 G4double x = DBL_MAX;
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954 DefineMaterial(track.GetMaterialCutsCouple());
00955 preStepKinEnergy = track.GetKineticEnergy();
00956 preStepScaledEnergy = preStepKinEnergy*massRatio;
00957 SelectModel(preStepScaledEnergy);
00958
00959 if(!currentModel->IsActive(preStepScaledEnergy)) { return x; }
00960
00961
00962 if(isIon) {
00963 G4double q2 = currentModel->ChargeSquareRatio(track);
00964 if(q2 != chargeSqRatio) {
00965 chargeSqRatio = q2;
00966 fFactor = q2*biasFactor*(*theDensityFactor)[currentCoupleIndex];
00967 reduceFactor = 1.0/(fFactor*massRatio);
00968 }
00969 }
00970
00971
00972
00973
00974
00975
00976
00977 if(biasManager) {
00978 if(0 == track.GetParentID()) {
00979 if(biasFlag && biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
00980 return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize);
00981 }
00982 }
00983 }
00984
00985
00986 if(preStepScaledEnergy < mfpKinEnergy) {
00987 if (integral) { ComputeLambdaForScaledEnergy(preStepScaledEnergy); }
00988 else { preStepLambda = GetLambdaForScaledEnergy(preStepScaledEnergy); }
00989
00990
00991 if(preStepLambda <= 0.0) {
00992 theNumberOfInteractionLengthLeft = -1.0;
00993 currentInteractionLength = DBL_MAX;
00994 }
00995 }
00996
00997
00998 if(preStepLambda > 0.0) {
00999 if (theNumberOfInteractionLengthLeft < 0.0) {
01000
01001
01002 ResetNumberOfInteractionLengthLeft();
01003
01004 } else if(currentInteractionLength < DBL_MAX) {
01005
01006
01007 theNumberOfInteractionLengthLeft -= previousStepSize/currentInteractionLength;
01008
01009 if(theNumberOfInteractionLengthLeft < 0.) {
01010 theNumberOfInteractionLengthLeft = 0.0;
01011
01012 }
01013 }
01014
01015
01016 currentInteractionLength = 1.0/preStepLambda;
01017 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
01018
01019 #ifdef G4VERBOSE
01020 if (verboseLevel>2){
01021
01022 G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength ";
01023 G4cout << "[ " << GetProcessName() << "]" << G4endl;
01024 G4cout << " for " << track.GetDefinition()->GetParticleName()
01025 << " in Material " << currentMaterial->GetName()
01026 << " Ekin(MeV)= " << preStepKinEnergy/MeV
01027 <<G4endl;
01028 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]"
01029 << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
01030 }
01031 #endif
01032 }
01033 return x;
01034 }
01035
01036
01037
01038 G4VParticleChange* G4VEnergyLossProcess::AlongStepDoIt(const G4Track& track,
01039 const G4Step& step)
01040 {
01041 fParticleChange.InitializeForAlongStep(track);
01042
01043 if(!isIonisation || !currentModel->IsActive(preStepScaledEnergy)) {
01044 return &fParticleChange;
01045 }
01046
01047
01048 G4double length = step.GetStepLength();
01049 if(length <= 0.0) { return &fParticleChange; }
01050 G4double eloss = 0.0;
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069 const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
01070
01071
01072 G4double weight = fParticleChange.GetParentWeight();
01073 if(weightFlag) {
01074 weight /= biasFactor;
01075 fParticleChange.ProposeWeight(weight);
01076 }
01077
01078
01079 if (length >= fRange) {
01080 eloss = preStepKinEnergy;
01081 if (useDeexcitation) {
01082 atomDeexcitation->AlongStepDeexcitation(scTracks, step,
01083 eloss, currentCoupleIndex);
01084 if(scTracks.size() > 0) { FillSecondariesAlongStep(eloss, weight); }
01085 if(eloss < 0.0) { eloss = 0.0; }
01086 }
01087 fParticleChange.SetProposedKineticEnergy(0.0);
01088 fParticleChange.ProposeLocalEnergyDeposit(eloss);
01089 return &fParticleChange;
01090 }
01091
01092
01093 eloss = GetDEDXForScaledEnergy(preStepScaledEnergy)*length;
01094
01095
01096 if(eloss > preStepKinEnergy*linLossLimit) {
01097
01098
01099
01100 G4double x = (fRange - length)/reduceFactor;
01101 eloss = preStepKinEnergy - ScaledKinEnergyForLoss(x)/massRatio;
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115 }
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129 G4double cut = (*theCuts)[currentCoupleIndex];
01130 G4double esec = 0.0;
01131
01132
01133 if(useSubCutoff) {
01134 if(idxSCoffRegions[currentCoupleIndex]) {
01135
01136 G4bool yes = false;
01137 G4StepPoint* prePoint = step.GetPreStepPoint();
01138
01139
01140 if(prePoint->GetStepStatus() == fGeomBoundary) { yes = true; }
01141
01142
01143 else {
01144 G4double preSafety = prePoint->GetSafety();
01145 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
01146
01147
01148 if(preSafety < rcut) {
01149 preSafety = safetyHelper->ComputeSafety(prePoint->GetPosition());
01150 }
01151
01152 if(preSafety < rcut) { yes = true; }
01153
01154
01155 else {
01156 G4double postSafety = preSafety - length;
01157 if(postSafety < rcut) {
01158 postSafety =
01159 safetyHelper->ComputeSafety(step.GetPostStepPoint()->GetPosition());
01160 if(postSafety < rcut) { yes = true; }
01161 }
01162 }
01163 }
01164
01165
01166 if(yes) {
01167
01168 cut = (*theSubCuts)[currentCoupleIndex];
01169 eloss -= GetSubDEDXForScaledEnergy(preStepScaledEnergy)*length;
01170 esec = SampleSubCutSecondaries(scTracks, step,
01171 currentModel,currentCoupleIndex);
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183 }
01184 }
01185 }
01186
01187
01188 if(isIon) {
01189 G4double eadd = 0.0;
01190 G4double eloss_before = eloss;
01191 currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
01192 eloss, eadd, length);
01193 if(eloss < 0.0) { eloss = 0.5*eloss_before; }
01194 }
01195
01196
01197 if (lossFluctuationFlag) {
01198 G4VEmFluctuationModel* fluc = currentModel->GetModelOfFluctuations();
01199 if(fluc &&
01200 (eloss + esec + lowestKinEnergy) < preStepKinEnergy) {
01201
01202 G4double tmax =
01203 std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
01204 G4double emean = eloss;
01205 eloss = fluc->SampleFluctuations(currentMaterial,dynParticle,
01206 tmax,length,emean);
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216 }
01217 }
01218
01219
01220 if (useDeexcitation) {
01221 G4double esecfluo = preStepKinEnergy - esec;
01222 G4double de = esecfluo;
01223
01224
01225
01226
01227
01228
01229 atomDeexcitation->AlongStepDeexcitation(scTracks, step,
01230 de, currentCoupleIndex);
01231
01232
01233 esecfluo -= de;
01234
01235
01236 if(eloss >= esecfluo) {
01237 esec += esecfluo;
01238 eloss -= esecfluo;
01239 } else {
01240 esec += esecfluo;
01241 eloss = 0.0;
01242 }
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253 }
01254 if(scTracks.size() > 0) { FillSecondariesAlongStep(eloss, weight); }
01255
01256
01257 G4double finalT = preStepKinEnergy - eloss - esec;
01258 if (finalT <= lowestKinEnergy) {
01259 eloss += finalT;
01260 finalT = 0.0;
01261 } else if(isIon) {
01262 fParticleChange.SetProposedCharge(
01263 currentModel->GetParticleCharge(track.GetParticleDefinition(),
01264 currentMaterial,finalT));
01265 }
01266
01267 if(eloss < 0.0) { eloss = 0.0; }
01268 fParticleChange.SetProposedKineticEnergy(finalT);
01269 fParticleChange.ProposeLocalEnergyDeposit(eloss);
01270
01271 if(1 < verboseLevel) {
01272 G4double del = finalT + eloss + esec - preStepKinEnergy;
01273 G4cout << "Final value eloss(MeV)= " << eloss/MeV
01274 << " preStepKinEnergy= " << preStepKinEnergy
01275 << " postStepKinEnergy= " << finalT
01276 << " de(keV)= " << del/keV
01277 << " lossFlag= " << lossFluctuationFlag
01278 << " status= " << track.GetTrackStatus()
01279 << G4endl;
01280 }
01281
01282 return &fParticleChange;
01283 }
01284
01285
01286
01287 void
01288 G4VEnergyLossProcess::FillSecondariesAlongStep(G4double&, G4double& weight)
01289 {
01290
01291 if(biasManager) {
01292 if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) {
01293 weight *=
01294 biasManager->ApplySecondaryBiasing(scTracks, currentCoupleIndex);
01295 }
01296 }
01297
01298
01299 G4int n = scTracks.size();
01300 fParticleChange.SetNumberOfSecondaries(n);
01301
01302 for(G4int i=0; i<n; ++i) {
01303 G4Track* t = scTracks[i];
01304 if(t) {
01305
01306 t->SetWeight(weight);
01307 pParticleChange->AddSecondary(t);
01308
01309
01310 }
01311 }
01312 scTracks.clear();
01313 }
01314
01315
01316
01317 G4double
01318 G4VEnergyLossProcess::SampleSubCutSecondaries(std::vector<G4Track*>& tracks,
01319 const G4Step& step,
01320 G4VEmModel* model,
01321 G4int idx)
01322 {
01323
01324 G4double esec = 0.0;
01325 G4double subcut = (*theSubCuts)[idx];
01326 G4double cut = (*theCuts)[idx];
01327 if(cut <= subcut) { return esec; }
01328
01329 const G4Track* track = step.GetTrack();
01330 const G4DynamicParticle* dp = track->GetDynamicParticle();
01331 G4double e = dp->GetKineticEnergy()*massRatio;
01332 G4double cross = (*theDensityFactor)[idx]*chargeSqRatio
01333 *(((*theSubLambdaTable)[(*theDensityIdx)[idx]])->Value(e));
01334 G4double length = step.GetStepLength();
01335
01336
01337 if(length*cross < perMillion) { return esec; }
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348 G4StepPoint* preStepPoint = step.GetPreStepPoint();
01349 G4StepPoint* postStepPoint = step.GetPostStepPoint();
01350 G4ThreeVector prepoint = preStepPoint->GetPosition();
01351 G4ThreeVector dr = postStepPoint->GetPosition() - prepoint;
01352 G4double pretime = preStepPoint->GetGlobalTime();
01353 G4double dt = postStepPoint->GetGlobalTime() - pretime;
01354
01355 G4double fragment = 0.0;
01356
01357 do {
01358 G4double del = -std::log(G4UniformRand())/cross;
01359 fragment += del/length;
01360 if (fragment > 1.0) break;
01361
01362
01363 secParticles.clear();
01364 model->SampleSecondaries(&secParticles,track->GetMaterialCutsCouple(),
01365 dp,subcut,cut);
01366
01367
01368 G4ThreeVector r = prepoint + fragment*dr;
01369 std::vector<G4DynamicParticle*>::iterator it;
01370 for(it=secParticles.begin(); it!=secParticles.end(); ++it) {
01371
01372 G4bool addSec = true;
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386 if(addSec) {
01387 G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
01388 t->SetTouchableHandle(track->GetTouchableHandle());
01389 tracks.push_back(t);
01390 esec += t->GetKineticEnergy();
01391 if (t->GetParticleDefinition() == thePositron) {
01392 esec += 2.0*electron_mass_c2;
01393 }
01394
01395
01396
01397
01398
01399
01400
01401
01402 }
01403 }
01404 } while (fragment <= 1.0);
01405 return esec;
01406 }
01407
01408
01409
01410 G4VParticleChange* G4VEnergyLossProcess::PostStepDoIt(const G4Track& track,
01411 const G4Step& step)
01412 {
01413
01414 theNumberOfInteractionLengthLeft = -1.0;
01415 mfpKinEnergy = DBL_MAX;
01416
01417 fParticleChange.InitializeForPostStep(track);
01418 G4double finalT = track.GetKineticEnergy();
01419 if(finalT <= lowestKinEnergy) { return &fParticleChange; }
01420
01421 G4double postStepScaledEnergy = finalT*massRatio;
01422
01423 if(!currentModel->IsActive(postStepScaledEnergy)) { return &fParticleChange; }
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433 if(biasFlag) {
01434 if(biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
01435 biasFlag = false;
01436 }
01437 }
01438
01439
01440 if (integral) {
01441 G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy);
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453 if(lx <= 0.0) {
01454 return &fParticleChange;
01455 } else if(preStepLambda*G4UniformRand() > lx) {
01456 return &fParticleChange;
01457 }
01458 }
01459
01460 SelectModel(postStepScaledEnergy);
01461
01462
01463 G4double weight = fParticleChange.GetParentWeight();
01464 if(weightFlag) {
01465 weight /= biasFactor;
01466 fParticleChange.ProposeWeight(weight);
01467 }
01468
01469 const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
01470 G4double tcut = (*theCuts)[currentCoupleIndex];
01471
01472
01473 secParticles.clear();
01474
01475 currentModel->SampleSecondaries(&secParticles, currentCouple, dynParticle, tcut);
01476
01477
01478 if(biasManager) {
01479 if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) {
01480 G4double eloss = 0.0;
01481 weight *= biasManager->ApplySecondaryBiasing(secParticles,
01482 track, currentModel,
01483 &fParticleChange, eloss,
01484 currentCoupleIndex, tcut,
01485 step.GetPostStepPoint()->GetSafety());
01486 if(eloss > 0.0) {
01487 eloss += fParticleChange.GetLocalEnergyDeposit();
01488 fParticleChange.ProposeLocalEnergyDeposit(eloss);
01489 }
01490 }
01491 }
01492
01493
01494 G4int num = secParticles.size();
01495 if(num > 0) {
01496
01497 fParticleChange.SetNumberOfSecondaries(num);
01498
01499 for (G4int i=0; i<num; ++i) {
01500 if(secParticles[i]) {
01501 G4Track* t = new G4Track(secParticles[i], track.GetGlobalTime(),
01502 track.GetPosition());
01503 t->SetTouchableHandle(track.GetTouchableHandle());
01504 t->SetWeight(weight);
01505
01506
01507 pParticleChange->AddSecondary(t);
01508 }
01509 }
01510 }
01511
01512 if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
01513 fAlive == fParticleChange.GetTrackStatus()) {
01514 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
01515 { fParticleChange.ProposeTrackStatus(fStopButAlive); }
01516 else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
01517 }
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531 return &fParticleChange;
01532 }
01533
01534
01535
01536 G4bool G4VEnergyLossProcess::StorePhysicsTable(
01537 const G4ParticleDefinition* part, const G4String& directory,
01538 G4bool ascii)
01539 {
01540 G4bool res = true;
01541 if ( baseParticle || part != particle ) return res;
01542
01543 if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
01544 {res = false;}
01545
01546 if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
01547 {res = false;}
01548
01549 if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
01550 {res = false;}
01551
01552 if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
01553 {res = false;}
01554
01555 if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
01556 {res = false;}
01557
01558 if(isIonisation &&
01559 !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
01560 {res = false;}
01561
01562 if(isIonisation &&
01563 !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
01564 {res = false;}
01565
01566 if(isIonisation &&
01567 !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
01568 {res = false;}
01569
01570 if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
01571 {res = false;}
01572
01573 if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
01574 {res = false;}
01575
01576 if ( res ) {
01577 if(0 < verboseLevel) {
01578 G4cout << "Physics tables are stored for " << particle->GetParticleName()
01579 << " and process " << GetProcessName()
01580 << " in the directory <" << directory
01581 << "> " << G4endl;
01582 }
01583 } else {
01584 G4cout << "Fail to store Physics Tables for "
01585 << particle->GetParticleName()
01586 << " and process " << GetProcessName()
01587 << " in the directory <" << directory
01588 << "> " << G4endl;
01589 }
01590 return res;
01591 }
01592
01593
01594
01595 G4bool
01596 G4VEnergyLossProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
01597 const G4String& directory,
01598 G4bool ascii)
01599 {
01600 G4bool res = true;
01601 const G4String particleName = part->GetParticleName();
01602
01603 if(1 < verboseLevel) {
01604 G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
01605 << particleName << " and process " << GetProcessName()
01606 << "; tables_are_built= " << tablesAreBuilt
01607 << G4endl;
01608 }
01609 if(particle == part) {
01610
01611 if ( !baseParticle ) {
01612
01613 G4bool fpi = true;
01614 if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
01615 {fpi = false;}
01616
01617
01618 if(!RetrieveTable(part,theDEDXTable,ascii,directory,"Ionisation",false))
01619 {fpi = false;}
01620
01621 if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
01622 {res = false;}
01623
01624 if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false))
01625 {res = false;}
01626
01627 if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false))
01628 {res = false;}
01629
01630 if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi))
01631 {res = false;}
01632
01633 if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
01634 {res = false;}
01635
01636 G4bool yes = false;
01637 if(nSCoffRegions > 0) {yes = true;}
01638
01639 if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
01640 {res = false;}
01641
01642 if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes))
01643 {res = false;}
01644
01645 if(!fpi) yes = false;
01646 if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes))
01647 {res = false;}
01648 }
01649 }
01650
01651 return res;
01652 }
01653
01654
01655
01656 G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part,
01657 G4PhysicsTable* aTable, G4bool ascii,
01658 const G4String& directory,
01659 const G4String& tname)
01660 {
01661 G4bool res = true;
01662 if ( aTable ) {
01663 const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
01664 if( !aTable->StorePhysicsTable(name,ascii)) res = false;
01665 }
01666 return res;
01667 }
01668
01669
01670
01671 G4bool
01672 G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part,
01673 G4PhysicsTable* aTable,
01674 G4bool ascii,
01675 const G4String& directory,
01676 const G4String& tname,
01677 G4bool mandatory)
01678 {
01679 G4bool isRetrieved = false;
01680 G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
01681 if(aTable) {
01682 if(aTable->ExistPhysicsTable(filename)) {
01683 if(G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii)) {
01684 isRetrieved = true;
01685 if((G4LossTableManager::Instance())->SplineFlag()) {
01686 size_t n = aTable->length();
01687 for(size_t i=0; i<n; ++i) {
01688 if((*aTable)[i]) { (*aTable)[i]->SetSpline(true); }
01689 }
01690 }
01691 if (0 < verboseLevel) {
01692 G4cout << tname << " table for " << part->GetParticleName()
01693 << " is Retrieved from <" << filename << ">"
01694 << G4endl;
01695 }
01696 }
01697 }
01698 }
01699 if(mandatory && !isRetrieved) {
01700 if(0 < verboseLevel) {
01701 G4cout << tname << " table for " << part->GetParticleName()
01702 << " from file <"
01703 << filename << "> is not Retrieved"
01704 << G4endl;
01705 }
01706 return false;
01707 }
01708 return true;
01709 }
01710
01711
01712
01713 G4double G4VEnergyLossProcess::GetDEDXDispersion(
01714 const G4MaterialCutsCouple *couple,
01715 const G4DynamicParticle* dp,
01716 G4double length)
01717 {
01718 DefineMaterial(couple);
01719 G4double ekin = dp->GetKineticEnergy();
01720 SelectModel(ekin*massRatio);
01721 G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
01722 tmax = std::min(tmax,(*theCuts)[currentCoupleIndex]);
01723 G4double d = 0.0;
01724 G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
01725 if(fm) { d = fm->Dispersion(currentMaterial,dp,tmax,length); }
01726 return d;
01727 }
01728
01729
01730
01731 G4double G4VEnergyLossProcess::CrossSectionPerVolume(
01732 G4double kineticEnergy, const G4MaterialCutsCouple* couple)
01733 {
01734
01735 DefineMaterial(couple);
01736 G4double cross = 0.0;
01737 if(theLambdaTable) {
01738 cross = (*theDensityFactor)[currentCoupleIndex]*
01739 ((*theLambdaTable)[basedCoupleIndex])->Value(kineticEnergy);
01740 } else {
01741 SelectModel(kineticEnergy);
01742 cross = currentModel->CrossSectionPerVolume(currentMaterial,
01743 particle, kineticEnergy,
01744 (*theCuts)[currentCoupleIndex]);
01745 }
01746 if(cross < 0.0) { cross = 0.0; }
01747 return cross;
01748 }
01749
01750
01751
01752 G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
01753 {
01754 DefineMaterial(track.GetMaterialCutsCouple());
01755 preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio);
01756 G4double x = DBL_MAX;
01757 if(0.0 < preStepLambda) { x = 1.0/preStepLambda; }
01758 return x;
01759 }
01760
01761
01762
01763 G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track,
01764 G4double x, G4double y,
01765 G4double& z)
01766 {
01767 G4GPILSelection sel;
01768 return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel);
01769 }
01770
01771
01772
01773 G4double G4VEnergyLossProcess::GetMeanFreePath(
01774 const G4Track& track,
01775 G4double,
01776 G4ForceCondition* condition)
01777
01778 {
01779 *condition = NotForced;
01780 return MeanFreePath(track);
01781 }
01782
01783
01784
01785 G4double G4VEnergyLossProcess::GetContinuousStepLimit(
01786 const G4Track&,
01787 G4double, G4double, G4double&)
01788 {
01789 return DBL_MAX;
01790 }
01791
01792
01793
01794 G4PhysicsVector*
01795 G4VEnergyLossProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* ,
01796 G4double )
01797 {
01798 G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nBins);
01799 v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
01800 return v;
01801 }
01802
01803
01804
01805 void G4VEnergyLossProcess::AddCollaborativeProcess(
01806 G4VEnergyLossProcess* p)
01807 {
01808 G4bool add = true;
01809 if(p->GetProcessName() != "eBrem") { add = false; }
01810 if(add && nProcesses > 0) {
01811 for(G4int i=0; i<nProcesses; ++i) {
01812 if(p == scProcesses[i]) {
01813 add = false;
01814 break;
01815 }
01816 }
01817 }
01818 if(add) {
01819 scProcesses.push_back(p);
01820 ++nProcesses;
01821 if (1 < verboseLevel) {
01822 G4cout << "### The process " << p->GetProcessName()
01823 << " is added to the list of collaborative processes of "
01824 << GetProcessName() << G4endl;
01825 }
01826 }
01827 }
01828
01829
01830
01831 void G4VEnergyLossProcess::SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType)
01832 {
01833 if(fTotal == tType && theDEDXunRestrictedTable != p && !baseParticle) {
01834 if(theDEDXunRestrictedTable) {
01835 theDEDXunRestrictedTable->clearAndDestroy();
01836 delete theDEDXunRestrictedTable;
01837 }
01838 theDEDXunRestrictedTable = p;
01839 if(p) {
01840 size_t n = p->length();
01841 G4PhysicsVector* pv = (*p)[0];
01842 G4double emax = maxKinEnergyCSDA;
01843
01844 for (size_t i=0; i<n; ++i) {
01845 G4double dedx = 0.0;
01846 pv = (*p)[i];
01847 if(pv) { dedx = pv->Value(emax); }
01848 else {
01849 pv = (*p)[(*theDensityIdx)[i]];
01850 if(pv) { dedx = pv->Value(emax)*(*theDensityFactor)[i]; }
01851 }
01852 theDEDXAtMaxEnergy[i] = dedx;
01853
01854
01855 }
01856 }
01857
01858 } else if(fRestricted == tType && theDEDXTable != p) {
01859
01860
01861
01862 if(theDEDXTable && !baseParticle) {
01863 if(theDEDXTable == theIonisationTable) { theIonisationTable = 0; }
01864 theDEDXTable->clearAndDestroy();
01865 delete theDEDXTable;
01866 }
01867 theDEDXTable = p;
01868 } else if(fSubRestricted == tType && theDEDXSubTable != p) {
01869 if(theDEDXSubTable && !baseParticle) {
01870 if(theDEDXSubTable == theIonisationSubTable) { theIonisationSubTable = 0; }
01871 theDEDXSubTable->clearAndDestroy();
01872 delete theDEDXSubTable;
01873 }
01874 theDEDXSubTable = p;
01875 } else if(fIsIonisation == tType && theIonisationTable != p) {
01876 if(theIonisationTable && theIonisationTable != theDEDXTable && !baseParticle) {
01877 theIonisationTable->clearAndDestroy();
01878 delete theIonisationTable;
01879 }
01880 theIonisationTable = p;
01881 } else if(fIsSubIonisation == tType && theIonisationSubTable != p) {
01882 if(theIonisationSubTable && theIonisationSubTable != theDEDXSubTable && !baseParticle) {
01883 theIonisationSubTable->clearAndDestroy();
01884 delete theIonisationSubTable;
01885 }
01886 theIonisationSubTable = p;
01887 }
01888 }
01889
01890
01891
01892 void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p)
01893 {
01894 if(theCSDARangeTable != p) { theCSDARangeTable = p; }
01895
01896 if(p) {
01897 size_t n = p->length();
01898 G4PhysicsVector* pv;
01899 G4double emax = maxKinEnergyCSDA;
01900
01901 for (size_t i=0; i<n; ++i) {
01902 pv = (*p)[i];
01903 G4double rmax = 0.0;
01904 if(pv) { rmax = pv->Value(emax); }
01905 else {
01906 pv = (*p)[(*theDensityIdx)[i]];
01907 if(pv) { rmax = pv->Value(emax)/(*theDensityFactor)[i]; }
01908 }
01909 theRangeAtMaxEnergy[i] = rmax;
01910
01911
01912 }
01913 }
01914 }
01915
01916
01917
01918 void G4VEnergyLossProcess::SetRangeTableForLoss(G4PhysicsTable* p)
01919 {
01920 if(theRangeTableForLoss != p) {
01921 theRangeTableForLoss = p;
01922 if(1 < verboseLevel) {
01923 G4cout << "### Set Range table " << p
01924 << " for " << particle->GetParticleName()
01925 << " and process " << GetProcessName() << G4endl;
01926 }
01927 }
01928 }
01929
01930
01931
01932 void G4VEnergyLossProcess::SetSecondaryRangeTable(G4PhysicsTable* p)
01933 {
01934 if(theSecondaryRangeTable != p) {
01935 theSecondaryRangeTable = p;
01936 if(1 < verboseLevel) {
01937 G4cout << "### Set SecondaryRange table " << p
01938 << " for " << particle->GetParticleName()
01939 << " and process " << GetProcessName() << G4endl;
01940 }
01941 }
01942 }
01943
01944
01945
01946 void G4VEnergyLossProcess::SetInverseRangeTable(G4PhysicsTable* p)
01947 {
01948 if(theInverseRangeTable != p) {
01949 theInverseRangeTable = p;
01950 if(1 < verboseLevel) {
01951 G4cout << "### Set InverseRange table " << p
01952 << " for " << particle->GetParticleName()
01953 << " and process " << GetProcessName() << G4endl;
01954 }
01955 }
01956 }
01957
01958
01959
01960 void G4VEnergyLossProcess::SetLambdaTable(G4PhysicsTable* p)
01961 {
01962 if(1 < verboseLevel) {
01963 G4cout << "### Set Lambda table " << p
01964 << " for " << particle->GetParticleName()
01965 << " and process " << GetProcessName() << G4endl;
01966 }
01967 if(theLambdaTable != p) { theLambdaTable = p; }
01968 tablesAreBuilt = true;
01969
01970 if(theLambdaTable) {
01971 size_t n = theLambdaTable->length();
01972 G4PhysicsVector* pv = (*theLambdaTable)[0];
01973 G4double e, ss, smax, emax;
01974
01975 size_t i;
01976
01977
01978 for (i=0; i<n; ++i) {
01979 pv = (*theLambdaTable)[i];
01980 if(pv) {
01981 size_t nb = pv->GetVectorLength();
01982 emax = DBL_MAX;
01983 smax = 0.0;
01984 if(nb > 0) {
01985 for (size_t j=0; j<nb; ++j) {
01986 e = pv->Energy(j);
01987 ss = (*pv)(j);
01988 if(ss > smax) {
01989 smax = ss;
01990 emax = e;
01991 }
01992 }
01993 }
01994 theEnergyOfCrossSectionMax[i] = emax;
01995 theCrossSectionMax[i] = smax;
01996 if(1 < verboseLevel) {
01997 G4cout << "For " << particle->GetParticleName()
01998 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
01999 << " lambda= " << smax << G4endl;
02000 }
02001 }
02002 }
02003
02004 for (i=0; i<n; ++i) {
02005 pv = (*theLambdaTable)[i];
02006 if(!pv){
02007 G4int j = (*theDensityIdx)[i];
02008 theEnergyOfCrossSectionMax[i] = theEnergyOfCrossSectionMax[j];
02009 theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
02010 }
02011 }
02012 }
02013 }
02014
02015
02016
02017 void G4VEnergyLossProcess::SetSubLambdaTable(G4PhysicsTable* p)
02018 {
02019 if(theSubLambdaTable != p) {
02020 theSubLambdaTable = p;
02021 if(1 < verboseLevel) {
02022 G4cout << "### Set SebLambda table " << p
02023 << " for " << particle->GetParticleName()
02024 << " and process " << GetProcessName() << G4endl;
02025 }
02026 }
02027 }
02028
02029
02030
02031 const G4Element* G4VEnergyLossProcess::GetCurrentElement() const
02032 {
02033 const G4Element* elm = 0;
02034 if(currentModel) { elm = currentModel->GetCurrentElement(); }
02035 return elm;
02036 }
02037
02038
02039
02040 void G4VEnergyLossProcess::SetCrossSectionBiasingFactor(G4double f,
02041 G4bool flag)
02042 {
02043 if(f > 0.0) {
02044 biasFactor = f;
02045 weightFlag = flag;
02046 if(1 < verboseLevel) {
02047 G4cout << "### SetCrossSectionBiasingFactor: for "
02048 << " process " << GetProcessName()
02049 << " biasFactor= " << f << " weightFlag= " << flag
02050 << G4endl;
02051 }
02052 }
02053 }
02054
02055
02056
02057 void
02058 G4VEnergyLossProcess::ActivateForcedInteraction(G4double length,
02059 const G4String& region,
02060 G4bool flag)
02061 {
02062 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
02063 if(1 < verboseLevel) {
02064 G4cout << "### ActivateForcedInteraction: for "
02065 << " process " << GetProcessName()
02066 << " length(mm)= " << length/mm
02067 << " in G4Region <" << region
02068 << "> weightFlag= " << flag
02069 << G4endl;
02070 }
02071 weightFlag = flag;
02072 biasManager->ActivateForcedInteraction(length, region);
02073 }
02074
02075
02076
02077 void
02078 G4VEnergyLossProcess::ActivateSecondaryBiasing(const G4String& region,
02079 G4double factor,
02080 G4double energyLimit)
02081 {
02082 if (0.0 <= factor) {
02083
02084
02085 if(0.0 == factor && secondaryParticle != G4Electron::Electron())
02086 { return; }
02087
02088 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
02089 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
02090 if(1 < verboseLevel) {
02091 G4cout << "### ActivateSecondaryBiasing: for "
02092 << " process " << GetProcessName()
02093 << " factor= " << factor
02094 << " in G4Region <" << region
02095 << "> energyLimit(MeV)= " << energyLimit/MeV
02096 << G4endl;
02097 }
02098 }
02099 }
02100
02101
02102