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 "G4VEmProcess.hh"
00066 #include "G4PhysicalConstants.hh"
00067 #include "G4SystemOfUnits.hh"
00068 #include "G4ProcessManager.hh"
00069 #include "G4LossTableManager.hh"
00070 #include "G4LossTableBuilder.hh"
00071 #include "G4Step.hh"
00072 #include "G4ParticleDefinition.hh"
00073 #include "G4VEmModel.hh"
00074 #include "G4DataVector.hh"
00075 #include "G4PhysicsTable.hh"
00076 #include "G4PhysicsLogVector.hh"
00077 #include "G4VParticleChange.hh"
00078 #include "G4ProductionCutsTable.hh"
00079 #include "G4Region.hh"
00080 #include "G4Gamma.hh"
00081 #include "G4Electron.hh"
00082 #include "G4Positron.hh"
00083 #include "G4PhysicsTableHelper.hh"
00084 #include "G4EmBiasingManager.hh"
00085 #include "G4GenericIon.hh"
00086
00087
00088
00089 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
00090 G4VDiscreteProcess(name, type),
00091 secondaryParticle(0),
00092 buildLambdaTable(true),
00093 numberOfModels(0),
00094 theLambdaTable(0),
00095 theLambdaTablePrim(0),
00096 theDensityFactor(0),
00097 theDensityIdx(0),
00098 integral(false),
00099 applyCuts(false),
00100 startFromNull(false),
00101 splineFlag(true),
00102 currentModel(0),
00103 particle(0),
00104 currentParticle(0),
00105 currentCouple(0)
00106 {
00107 SetVerboseLevel(1);
00108
00109
00110 minKinEnergy = 0.1*keV;
00111 maxKinEnergy = 10.0*TeV;
00112 nLambdaBins = 77;
00113 minKinEnergyPrim = DBL_MAX;
00114
00115
00116 lambdaFactor = 0.8;
00117
00118
00119 polarAngleLimit = 0.0;
00120 biasFactor = 1.0;
00121
00122
00123 theGamma = G4Gamma::Gamma();
00124 theElectron = G4Electron::Electron();
00125 thePositron = G4Positron::Positron();
00126
00127 pParticleChange = &fParticleChange;
00128 secParticles.reserve(5);
00129
00130 preStepLambda = 0.0;
00131 mfpKinEnergy = DBL_MAX;
00132
00133 modelManager = new G4EmModelManager();
00134 biasManager = 0;
00135 biasFlag = false;
00136 weightFlag = false;
00137 (G4LossTableManager::Instance())->Register(this);
00138 warn = 0;
00139 }
00140
00141
00142
00143 G4VEmProcess::~G4VEmProcess()
00144 {
00145 if(1 < verboseLevel) {
00146 G4cout << "G4VEmProcess destruct " << GetProcessName()
00147 << " " << this << " " << theLambdaTable <<G4endl;
00148 }
00149 Clear();
00150 if(theLambdaTable) {
00151 theLambdaTable->clearAndDestroy();
00152 delete theLambdaTable;
00153 }
00154 if(theLambdaTablePrim) {
00155 theLambdaTablePrim->clearAndDestroy();
00156 delete theLambdaTablePrim;
00157 }
00158 delete modelManager;
00159 delete biasManager;
00160 (G4LossTableManager::Instance())->DeRegister(this);
00161 }
00162
00163
00164
00165 void G4VEmProcess::Clear()
00166 {
00167 currentCouple = 0;
00168 preStepLambda = 0.0;
00169 mfpKinEnergy = DBL_MAX;
00170 }
00171
00172
00173
00174 G4double G4VEmProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
00175 const G4Material*)
00176 {
00177 return 0.0;
00178 }
00179
00180
00181
00182 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p,
00183 const G4Region* region)
00184 {
00185 G4VEmFluctuationModel* fm = 0;
00186 modelManager->AddEmModel(order, p, fm, region);
00187 if(p) { p->SetParticleChange(pParticleChange); }
00188 }
00189
00190
00191
00192 void G4VEmProcess::SetModel(G4VEmModel* p, G4int index)
00193 {
00194 ++warn;
00195 if(warn < 10) {
00196 G4cout << "### G4VEmProcess::SetModel is obsolete method and will be "
00197 << "removed for the next release." << G4endl;
00198 G4cout << " Please, use SetEmModel" << G4endl;
00199 }
00200 G4int n = emModels.size();
00201 if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
00202 emModels[index] = p;
00203 }
00204
00205
00206
00207 G4VEmModel* G4VEmProcess::Model(G4int index)
00208 {
00209 if(warn < 10) {
00210 G4cout << "### G4VEmProcess::Model is obsolete method and will be "
00211 << "removed for the next release." << G4endl;
00212 G4cout << " Please, use EmModel" << G4endl;
00213 }
00214 G4VEmModel* p = 0;
00215 if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; }
00216 return p;
00217 }
00218
00219
00220
00221 void G4VEmProcess::SetEmModel(G4VEmModel* p, G4int index)
00222 {
00223 G4int n = emModels.size();
00224 if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
00225 emModels[index] = p;
00226 }
00227
00228
00229
00230 G4VEmModel* G4VEmProcess::EmModel(G4int index)
00231 {
00232 G4VEmModel* p = 0;
00233 if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; }
00234 return p;
00235 }
00236
00237
00238
00239 void G4VEmProcess::UpdateEmModel(const G4String& nam,
00240 G4double emin, G4double emax)
00241 {
00242 modelManager->UpdateEmModel(nam, emin, emax);
00243 }
00244
00245
00246
00247 G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver)
00248 {
00249 return modelManager->GetModel(idx, ver);
00250 }
00251
00252
00253
00254 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
00255 {
00256 if(!particle) { SetParticle(&part); }
00257
00258 if(part.GetParticleType() == "nucleus" &&
00259 part.GetParticleSubType() == "generic") {
00260
00261 G4String pname = part.GetParticleName();
00262 if(pname != "deuteron" && pname != "triton" &&
00263 pname != "alpha" && pname != "He3" &&
00264 pname != "alpha+" && pname != "helium" &&
00265 pname != "hydrogen") {
00266
00267 particle = G4GenericIon::GenericIon();
00268 }
00269 }
00270
00271 if(1 < verboseLevel) {
00272 G4cout << "G4VEmProcess::PreparePhysicsTable() for "
00273 << GetProcessName()
00274 << " and particle " << part.GetParticleName()
00275 << " local particle " << particle->GetParticleName()
00276 << G4endl;
00277 }
00278
00279 G4LossTableManager* man = G4LossTableManager::Instance();
00280 G4LossTableBuilder* bld = man->GetTableBuilder();
00281
00282 man->PreparePhysicsTable(&part, this);
00283
00284 if(particle == &part) {
00285 Clear();
00286 InitialiseProcess(particle);
00287
00288 const G4ProductionCutsTable* theCoupleTable=
00289 G4ProductionCutsTable::GetProductionCutsTable();
00290 size_t n = theCoupleTable->GetTableSize();
00291
00292 theEnergyOfCrossSectionMax.resize(n, 0.0);
00293 theCrossSectionMax.resize(n, DBL_MAX);
00294
00295
00296 numberOfModels = modelManager->NumberOfModels();
00297 for(G4int i=0; i<numberOfModels; ++i) {
00298 G4VEmModel* mod = modelManager->GetModel(i);
00299 if(0 == i) { currentModel = mod; }
00300 mod->SetPolarAngleLimit(polarAngleLimit);
00301 if(mod->HighEnergyLimit() > maxKinEnergy) {
00302 mod->SetHighEnergyLimit(maxKinEnergy);
00303 }
00304 }
00305
00306 if(man->AtomDeexcitation()) { modelManager->SetFluoFlag(true); }
00307 theCuts = modelManager->Initialise(particle,secondaryParticle,
00308 2.,verboseLevel);
00309 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
00310 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
00311 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
00312
00313
00314 if(buildLambdaTable){
00315 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
00316 bld->InitialiseBaseMaterials(theLambdaTable);
00317 }
00318
00319 if(minKinEnergyPrim < maxKinEnergy){
00320 theLambdaTablePrim =
00321 G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTablePrim);
00322 bld->InitialiseBaseMaterials(theLambdaTablePrim);
00323 }
00324
00325 if(biasManager) {
00326 biasManager->Initialise(part,GetProcessName(),verboseLevel);
00327 biasFlag = false;
00328 }
00329 }
00330 theDensityFactor = bld->GetDensityFactors();
00331 theDensityIdx = bld->GetCoupleIndexes();
00332 }
00333
00334
00335
00336 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
00337 {
00338 G4String num = part.GetParticleName();
00339 if(1 < verboseLevel) {
00340 G4cout << "G4VEmProcess::BuildPhysicsTable() for "
00341 << GetProcessName()
00342 << " and particle " << num
00343 << " buildLambdaTable= " << buildLambdaTable
00344 << G4endl;
00345 }
00346
00347 (G4LossTableManager::Instance())->BuildPhysicsTable(particle);
00348
00349 if(buildLambdaTable || minKinEnergyPrim < maxKinEnergy) {
00350 BuildLambdaTable();
00351 }
00352
00353
00354 if(1 < verboseLevel ||
00355 (0 < verboseLevel && (num == "gamma" || num == "e-" ||
00356 num == "e+" || num == "mu+" ||
00357 num == "mu-" || num == "proton"||
00358 num == "pi+" || num == "pi-" ||
00359 num == "kaon+" || num == "kaon-" ||
00360 num == "alpha" || num == "anti_proton" ||
00361 num == "GenericIon")))
00362 {
00363 particle = ∂
00364 PrintInfoDefinition();
00365 }
00366
00367 if(1 < verboseLevel) {
00368 G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
00369 << GetProcessName()
00370 << " and particle " << num
00371 << G4endl;
00372 }
00373 }
00374
00375
00376
00377 void G4VEmProcess::BuildLambdaTable()
00378 {
00379 if(1 < verboseLevel) {
00380 G4cout << "G4EmProcess::BuildLambdaTable() for process "
00381 << GetProcessName() << " and particle "
00382 << particle->GetParticleName() << " " << this
00383 << G4endl;
00384 }
00385
00386
00387 const G4ProductionCutsTable* theCoupleTable=
00388 G4ProductionCutsTable::GetProductionCutsTable();
00389 size_t numOfCouples = theCoupleTable->GetTableSize();
00390
00391 G4LossTableBuilder* bld = (G4LossTableManager::Instance())->GetTableBuilder();
00392
00393 G4PhysicsLogVector* aVector = 0;
00394 G4PhysicsLogVector* bVector = 0;
00395 G4PhysicsLogVector* aVectorPrim = 0;
00396 G4PhysicsLogVector* bVectorPrim = 0;
00397
00398 G4double scale = 1.0;
00399 G4double emax1 = maxKinEnergy;
00400 if(startFromNull || minKinEnergyPrim < maxKinEnergy ) {
00401 scale = std::log(maxKinEnergy/minKinEnergy);
00402 if(minKinEnergyPrim < maxKinEnergy) { emax1 = minKinEnergyPrim; }
00403 }
00404
00405 for(size_t i=0; i<numOfCouples; ++i) {
00406
00407 if (bld->GetFlag(i)) {
00408
00409
00410 const G4MaterialCutsCouple* couple =
00411 theCoupleTable->GetMaterialCutsCouple(i);
00412
00413
00414 if(buildLambdaTable) {
00415 delete (*theLambdaTable)[i];
00416
00417 G4bool startNull = startFromNull;
00418
00419 if(startFromNull || minKinEnergyPrim < maxKinEnergy) {
00420 G4double emin = MinPrimaryEnergy(particle,couple->GetMaterial());
00421 if(emin < minKinEnergy) {
00422 emin = minKinEnergy;
00423 startNull = false;
00424 }
00425 G4double emax = emax1;
00426 if(emax <= emin) { emax = 2*emin; }
00427 G4int bin =
00428 G4lrint(nLambdaBins*std::log(emax/emin)/scale);
00429 if(bin < 3) { bin = 3; }
00430 aVector = new G4PhysicsLogVector(emin, emax, bin);
00431
00432
00433 } else if(!bVector) {
00434 aVector =
00435 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
00436 bVector = aVector;
00437 } else {
00438 aVector = new G4PhysicsLogVector(*bVector);
00439 }
00440 aVector->SetSpline(splineFlag);
00441 modelManager->FillLambdaVector(aVector, couple, startNull);
00442 if(splineFlag) { aVector->FillSecondDerivatives(); }
00443 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
00444 }
00445
00446 if(minKinEnergyPrim < maxKinEnergy) {
00447 delete (*theLambdaTablePrim)[i];
00448
00449
00450 if(!bVectorPrim) {
00451 G4int bin =
00452 G4lrint(nLambdaBins*std::log(maxKinEnergy/minKinEnergyPrim)/scale);
00453 if(bin < 3) { bin = 3; }
00454 aVectorPrim =
00455 new G4PhysicsLogVector(minKinEnergyPrim, maxKinEnergy, bin);
00456 bVectorPrim = aVectorPrim;
00457 } else {
00458 aVectorPrim = new G4PhysicsLogVector(*bVectorPrim);
00459 }
00460
00461 aVectorPrim->SetSpline(true);
00462 modelManager->FillLambdaVector(aVectorPrim, couple, false,
00463 fIsCrossSectionPrim);
00464 aVectorPrim->FillSecondDerivatives();
00465 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTablePrim, i, aVectorPrim);
00466 }
00467 }
00468 }
00469
00470 if(buildLambdaTable) { FindLambdaMax(); }
00471
00472 if(1 < verboseLevel) {
00473 G4cout << "Lambda table is built for "
00474 << particle->GetParticleName()
00475 << G4endl;
00476 }
00477 }
00478
00479
00480
00481 void G4VEmProcess::PrintInfoDefinition()
00482 {
00483 if(verboseLevel > 0) {
00484 G4cout << G4endl << GetProcessName() << ": for "
00485 << particle->GetParticleName();
00486 if(integral) { G4cout << ", integral: 1 "; }
00487 if(applyCuts) { G4cout << ", applyCuts: 1 "; }
00488 G4cout << " SubType= " << GetProcessSubType();;
00489 if(biasFactor != 1.0) { G4cout << " BiasingFactor= " << biasFactor; }
00490 G4cout << G4endl;
00491 if(buildLambdaTable) {
00492 size_t length = theLambdaTable->length();
00493 for(size_t i=0; i<length; ++i) {
00494 G4PhysicsVector* v = (*theLambdaTable)[i];
00495 if(v) {
00496 G4cout << " Lambda table from "
00497 << G4BestUnit(minKinEnergy,"Energy")
00498 << " to "
00499 << G4BestUnit(v->GetMaxEnergy(),"Energy")
00500 << " in " << v->GetVectorLength()-1
00501 << " bins, spline: "
00502 << splineFlag
00503 << G4endl;
00504 break;
00505 }
00506 }
00507 }
00508 if(minKinEnergyPrim < maxKinEnergy) {
00509 size_t length = theLambdaTablePrim->length();
00510 for(size_t i=0; i<length; ++i) {
00511 G4PhysicsVector* v = (*theLambdaTablePrim)[i];
00512 if(v) {
00513 G4cout << " LambdaPrime table from "
00514 << G4BestUnit(minKinEnergyPrim,"Energy")
00515 << " to "
00516 << G4BestUnit(maxKinEnergy,"Energy")
00517 << " in " << v->GetVectorLength()-1
00518 << " bins "
00519 << G4endl;
00520 break;
00521 }
00522 }
00523 }
00524 PrintInfo();
00525 modelManager->DumpModelList(verboseLevel);
00526 }
00527
00528 if(verboseLevel > 2 && buildLambdaTable) {
00529 G4cout << " LambdaTable address= " << theLambdaTable << G4endl;
00530 if(theLambdaTable) { G4cout << (*theLambdaTable) << G4endl; }
00531 }
00532 }
00533
00534
00535
00536 void G4VEmProcess::StartTracking(G4Track* track)
00537 {
00538
00539 currentParticle = track->GetParticleDefinition();
00540 theNumberOfInteractionLengthLeft = -1.0;
00541
00542
00543 mfpKinEnergy = DBL_MAX;
00544
00545
00546 if(biasManager) {
00547 if(0 == track->GetParentID()) {
00548
00549 biasFlag = true;
00550 biasManager->ResetForcedInteraction();
00551 }
00552 }
00553 }
00554
00555
00556
00557 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
00558 const G4Track& track,
00559 G4double previousStepSize,
00560 G4ForceCondition* condition)
00561 {
00562 *condition = NotForced;
00563 G4double x = DBL_MAX;
00564
00565 preStepKinEnergy = track.GetKineticEnergy();
00566 DefineMaterial(track.GetMaterialCutsCouple());
00567 SelectModel(preStepKinEnergy, currentCoupleIndex);
00568
00569 if(!currentModel->IsActive(preStepKinEnergy)) { return x; }
00570
00571
00572 if(biasManager) {
00573 if(0 == track.GetParentID()) {
00574 if(biasFlag && biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
00575 return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize);
00576 }
00577 }
00578 }
00579
00580
00581 if(preStepKinEnergy < mfpKinEnergy) {
00582 if (integral) { ComputeIntegralLambda(preStepKinEnergy); }
00583 else { preStepLambda = GetCurrentLambda(preStepKinEnergy); }
00584
00585
00586 if(preStepLambda <= 0.0) {
00587 theNumberOfInteractionLengthLeft = -1.0;
00588 currentInteractionLength = DBL_MAX;
00589 }
00590 }
00591
00592
00593 if(preStepLambda > 0.0) {
00594
00595 if (theNumberOfInteractionLengthLeft < 0.0) {
00596
00597
00598 ResetNumberOfInteractionLengthLeft();
00599
00600 } else if(currentInteractionLength < DBL_MAX) {
00601
00602
00603 theNumberOfInteractionLengthLeft -= previousStepSize/currentInteractionLength;
00604
00605 if(theNumberOfInteractionLengthLeft < 0.) {
00606 theNumberOfInteractionLengthLeft = 0.0;
00607
00608 }
00609 }
00610
00611
00612 currentInteractionLength = 1.0/preStepLambda;
00613 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
00614 #ifdef G4VERBOSE
00615 if (verboseLevel>2){
00616 G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
00617 G4cout << "[ " << GetProcessName() << "]" << G4endl;
00618 G4cout << " for " << currentParticle->GetParticleName()
00619 << " in Material " << currentMaterial->GetName()
00620 << " Ekin(MeV)= " << preStepKinEnergy/MeV
00621 <<G4endl;
00622 G4cout << " MeanFreePath = " << currentInteractionLength/cm << "[cm]"
00623 << " InteractionLength= " << x/cm <<"[cm] " <<G4endl;
00624 }
00625 #endif
00626 }
00627 return x;
00628 }
00629
00630
00631
00632 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
00633 const G4Step& step)
00634 {
00635
00636 theNumberOfInteractionLengthLeft = -1.0;
00637 mfpKinEnergy = DBL_MAX;
00638
00639 fParticleChange.InitializeForPostStep(track);
00640
00641
00642
00643 if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
00644
00645 G4double finalT = track.GetKineticEnergy();
00646
00647
00648 if(biasFlag) {
00649 if(biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
00650 biasFlag = false;
00651 }
00652 }
00653
00654
00655 if (integral) {
00656 G4double lx = GetLambda(finalT, currentCouple);
00657 if(preStepLambda<lx && 1 < verboseLevel) {
00658 G4cout << "WARNING: for " << currentParticle->GetParticleName()
00659 << " and " << GetProcessName()
00660 << " E(MeV)= " << finalT/MeV
00661 << " preLambda= " << preStepLambda << " < "
00662 << lx << " (postLambda) "
00663 << G4endl;
00664 }
00665
00666 if(preStepLambda*G4UniformRand() > lx) {
00667 ClearNumberOfInteractionLengthLeft();
00668 return &fParticleChange;
00669 }
00670 }
00671
00672 SelectModel(finalT, currentCoupleIndex);
00673 if(!currentModel->IsActive(finalT)) { return &fParticleChange; }
00674
00675
00676 G4double weight = fParticleChange.GetParentWeight();
00677 if(weightFlag) {
00678 weight /= biasFactor;
00679 fParticleChange.ProposeWeight(weight);
00680 }
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693 secParticles.clear();
00694 currentModel->SampleSecondaries(&secParticles,
00695 currentCouple,
00696 track.GetDynamicParticle(),
00697 (*theCuts)[currentCoupleIndex]);
00698
00699
00700 if(biasManager) {
00701 if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) {
00702 G4double eloss = 0.0;
00703 weight *= biasManager->ApplySecondaryBiasing(secParticles,
00704 track, currentModel,
00705 &fParticleChange,
00706 eloss, currentCoupleIndex,
00707 (*theCuts)[currentCoupleIndex],
00708 step.GetPostStepPoint()->GetSafety());
00709 if(eloss > 0.0) {
00710 eloss += fParticleChange.GetLocalEnergyDeposit();
00711 fParticleChange.ProposeLocalEnergyDeposit(eloss);
00712 }
00713 }
00714 }
00715
00716
00717 G4int num = secParticles.size();
00718 if(num > 0) {
00719
00720 fParticleChange.SetNumberOfSecondaries(num);
00721 G4double edep = fParticleChange.GetLocalEnergyDeposit();
00722
00723 for (G4int i=0; i<num; ++i) {
00724 if (secParticles[i]) {
00725 G4DynamicParticle* dp = secParticles[i];
00726 const G4ParticleDefinition* p = dp->GetParticleDefinition();
00727 G4double e = dp->GetKineticEnergy();
00728 G4bool good = true;
00729 if(applyCuts) {
00730 if (p == theGamma) {
00731 if (e < (*theCutsGamma)[currentCoupleIndex]) { good = false; }
00732
00733 } else if (p == theElectron) {
00734 if (e < (*theCutsElectron)[currentCoupleIndex]) { good = false; }
00735
00736 } else if (p == thePositron) {
00737 if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
00738 e < (*theCutsPositron)[currentCoupleIndex]) {
00739 good = false;
00740 e += 2.0*electron_mass_c2;
00741 }
00742 }
00743
00744 }
00745 if (good) {
00746 G4Track* t = new G4Track(dp, track.GetGlobalTime(), track.GetPosition());
00747 t->SetTouchableHandle(track.GetTouchableHandle());
00748 t->SetWeight(weight);
00749 pParticleChange->AddSecondary(t);
00750
00751
00752 } else {
00753 delete dp;
00754 edep += e;
00755 }
00756 }
00757 }
00758 fParticleChange.ProposeLocalEnergyDeposit(edep);
00759 }
00760
00761 if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
00762 fAlive == fParticleChange.GetTrackStatus()) {
00763 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
00764 { fParticleChange.ProposeTrackStatus(fStopButAlive); }
00765 else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
00766 }
00767
00768
00769 return &fParticleChange;
00770 }
00771
00772
00773
00774 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
00775 const G4String& directory,
00776 G4bool ascii)
00777 {
00778 G4bool yes = true;
00779
00780 if ( theLambdaTable && part == particle) {
00781 const G4String name =
00782 GetPhysicsTableFileName(part,directory,"Lambda",ascii);
00783 yes = theLambdaTable->StorePhysicsTable(name,ascii);
00784
00785 if ( yes ) {
00786 G4cout << "Physics table is stored for " << particle->GetParticleName()
00787 << " and process " << GetProcessName()
00788 << " in the directory <" << directory
00789 << "> " << G4endl;
00790 } else {
00791 G4cout << "Fail to store Physics Table for "
00792 << particle->GetParticleName()
00793 << " and process " << GetProcessName()
00794 << " in the directory <" << directory
00795 << "> " << G4endl;
00796 }
00797 }
00798 if ( theLambdaTablePrim && part == particle) {
00799 const G4String name =
00800 GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
00801 yes = theLambdaTablePrim->StorePhysicsTable(name,ascii);
00802
00803 if ( yes ) {
00804 G4cout << "Physics table prim is stored for " << particle->GetParticleName()
00805 << " and process " << GetProcessName()
00806 << " in the directory <" << directory
00807 << "> " << G4endl;
00808 } else {
00809 G4cout << "Fail to store Physics Table Prim for "
00810 << particle->GetParticleName()
00811 << " and process " << GetProcessName()
00812 << " in the directory <" << directory
00813 << "> " << G4endl;
00814 }
00815 }
00816 return yes;
00817 }
00818
00819
00820
00821 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
00822 const G4String& directory,
00823 G4bool ascii)
00824 {
00825 if(1 < verboseLevel) {
00826 G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
00827 << part->GetParticleName() << " and process "
00828 << GetProcessName() << G4endl;
00829 }
00830 G4bool yes = true;
00831
00832 if((!buildLambdaTable && minKinEnergyPrim > maxKinEnergy)
00833 || particle != part) { return yes; }
00834
00835 const G4String particleName = part->GetParticleName();
00836 G4String filename;
00837
00838 if(buildLambdaTable) {
00839 filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
00840 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,
00841 filename,ascii);
00842 if ( yes ) {
00843 if (0 < verboseLevel) {
00844 G4cout << "Lambda table for " << particleName
00845 << " is Retrieved from <"
00846 << filename << ">"
00847 << G4endl;
00848 }
00849 if((G4LossTableManager::Instance())->SplineFlag()) {
00850 size_t n = theLambdaTable->length();
00851 for(size_t i=0; i<n; ++i) {
00852 if((* theLambdaTable)[i]) {
00853 (* theLambdaTable)[i]->SetSpline(true);
00854 }
00855 }
00856 }
00857 } else {
00858 if (1 < verboseLevel) {
00859 G4cout << "Lambda table for " << particleName << " in file <"
00860 << filename << "> is not exist"
00861 << G4endl;
00862 }
00863 }
00864 }
00865 if(minKinEnergyPrim < maxKinEnergy) {
00866 filename = GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
00867 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTablePrim,
00868 filename,ascii);
00869 if ( yes ) {
00870 if (0 < verboseLevel) {
00871 G4cout << "Lambda table prim for " << particleName
00872 << " is Retrieved from <"
00873 << filename << ">"
00874 << G4endl;
00875 }
00876 if((G4LossTableManager::Instance())->SplineFlag()) {
00877 size_t n = theLambdaTablePrim->length();
00878 for(size_t i=0; i<n; ++i) {
00879 if((* theLambdaTablePrim)[i]) {
00880 (* theLambdaTablePrim)[i]->SetSpline(true);
00881 }
00882 }
00883 }
00884 } else {
00885 if (1 < verboseLevel) {
00886 G4cout << "Lambda table prim for " << particleName << " in file <"
00887 << filename << "> is not exist"
00888 << G4endl;
00889 }
00890 }
00891 }
00892
00893 return yes;
00894 }
00895
00896
00897
00898 G4double
00899 G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
00900 const G4MaterialCutsCouple* couple)
00901 {
00902
00903 DefineMaterial(couple);
00904 G4double cross = 0.0;
00905 if(theLambdaTable) {
00906 cross = (*theDensityFactor)[currentCoupleIndex]*
00907 (((*theLambdaTable)[basedCoupleIndex])->Value(kineticEnergy));
00908 } else {
00909 SelectModel(kineticEnergy, currentCoupleIndex);
00910 cross = currentModel->CrossSectionPerVolume(currentMaterial,
00911 currentParticle,kineticEnergy);
00912 }
00913
00914 if(cross < 0.0) { cross = 0.0; }
00915 return cross;
00916 }
00917
00918
00919
00920 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
00921 G4double,
00922 G4ForceCondition* condition)
00923 {
00924 *condition = NotForced;
00925 return G4VEmProcess::MeanFreePath(track);
00926 }
00927
00928
00929
00930 G4double G4VEmProcess::MeanFreePath(const G4Track& track)
00931 {
00932 DefineMaterial(track.GetMaterialCutsCouple());
00933 preStepLambda = GetCurrentLambda(track.GetKineticEnergy());
00934 G4double x = DBL_MAX;
00935 if(0.0 < preStepLambda) { x = 1.0/preStepLambda; }
00936 return x;
00937 }
00938
00939
00940
00941 G4double
00942 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kineticEnergy,
00943 G4double Z, G4double A, G4double cut)
00944 {
00945 SelectModel(kineticEnergy, currentCoupleIndex);
00946 G4double x = 0.0;
00947 if(currentModel) {
00948 x = currentModel->ComputeCrossSectionPerAtom(currentParticle,kineticEnergy,
00949 Z,A,cut);
00950 }
00951 return x;
00952 }
00953
00954
00955
00956 void G4VEmProcess::FindLambdaMax()
00957 {
00958 if(1 < verboseLevel) {
00959 G4cout << "### G4VEmProcess::FindLambdaMax: "
00960 << particle->GetParticleName()
00961 << " and process " << GetProcessName() << G4endl;
00962 }
00963 size_t n = theLambdaTable->length();
00964 G4PhysicsVector* pv;
00965 G4double e, ss, emax, smax;
00966
00967 size_t i;
00968
00969
00970 for (i=0; i<n; ++i) {
00971 pv = (*theLambdaTable)[i];
00972 if(pv) {
00973 size_t nb = pv->GetVectorLength();
00974 emax = DBL_MAX;
00975 smax = 0.0;
00976 if(nb > 0) {
00977 for (size_t j=0; j<nb; ++j) {
00978 e = pv->Energy(j);
00979 ss = (*pv)(j);
00980 if(ss > smax) {
00981 smax = ss;
00982 emax = e;
00983 }
00984 }
00985 }
00986 theEnergyOfCrossSectionMax[i] = emax;
00987 theCrossSectionMax[i] = smax;
00988 if(1 < verboseLevel) {
00989 G4cout << "For " << particle->GetParticleName()
00990 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
00991 << " lambda= " << smax << G4endl;
00992 }
00993 }
00994 }
00995
00996 for (i=0; i<n; ++i) {
00997 pv = (*theLambdaTable)[i];
00998 if(!pv){
00999 G4int j = (*theDensityIdx)[i];
01000 theEnergyOfCrossSectionMax[i] = theEnergyOfCrossSectionMax[j];
01001 theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
01002 }
01003 }
01004 }
01005
01006
01007
01008 G4PhysicsVector*
01009 G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
01010 {
01011 G4PhysicsVector* v =
01012 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
01013 v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
01014 return v;
01015 }
01016
01017
01018
01019 const G4Element* G4VEmProcess::GetCurrentElement() const
01020 {
01021 const G4Element* elm = 0;
01022 if(currentModel) {elm = currentModel->GetCurrentElement(); }
01023 return elm;
01024 }
01025
01026
01027
01028 void G4VEmProcess::SetCrossSectionBiasingFactor(G4double f, G4bool flag)
01029 {
01030 if(f > 0.0) {
01031 biasFactor = f;
01032 weightFlag = flag;
01033 if(1 < verboseLevel) {
01034 G4cout << "### SetCrossSectionBiasingFactor: for "
01035 << particle->GetParticleName()
01036 << " and process " << GetProcessName()
01037 << " biasFactor= " << f << " weightFlag= " << flag
01038 << G4endl;
01039 }
01040 }
01041 }
01042
01043
01044
01045 void
01046 G4VEmProcess::ActivateForcedInteraction(G4double length, const G4String& r,
01047 G4bool flag)
01048 {
01049 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
01050 if(1 < verboseLevel) {
01051 G4cout << "### ActivateForcedInteraction: for "
01052 << particle->GetParticleName()
01053 << " and process " << GetProcessName()
01054 << " length(mm)= " << length/mm
01055 << " in G4Region <" << r
01056 << "> weightFlag= " << flag
01057 << G4endl;
01058 }
01059 weightFlag = flag;
01060 biasManager->ActivateForcedInteraction(length, r);
01061 }
01062
01063
01064
01065 void
01066 G4VEmProcess::ActivateSecondaryBiasing(const G4String& region,
01067 G4double factor,
01068 G4double energyLimit)
01069 {
01070 if (0.0 <= factor) {
01071
01072
01073 if(0.0 == factor && secondaryParticle != G4Electron::Electron())
01074 { return; }
01075
01076 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
01077 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
01078 if(1 < verboseLevel) {
01079 G4cout << "### ActivateSecondaryBiasing: for "
01080 << " process " << GetProcessName()
01081 << " factor= " << factor
01082 << " in G4Region <" << region
01083 << "> energyLimit(MeV)= " << energyLimit/MeV
01084 << G4endl;
01085 }
01086 }
01087 }
01088
01089
01090
01091 void G4VEmProcess::SetMinKinEnergy(G4double e)
01092 {
01093 nLambdaBins = G4lrint(nLambdaBins*std::log(maxKinEnergy/e)
01094 /std::log(maxKinEnergy/minKinEnergy));
01095 minKinEnergy = e;
01096 }
01097
01098
01099
01100 void G4VEmProcess::SetMaxKinEnergy(G4double e)
01101 {
01102 nLambdaBins = G4lrint(nLambdaBins*std::log(e/minKinEnergy)
01103 /std::log(maxKinEnergy/minKinEnergy));
01104 maxKinEnergy = e;
01105 }
01106
01107