93 buildLambdaTable(true),
96 theLambdaTablePrim(0),
101 startFromNull(false),
111 minKinEnergy = 0.1*
keV;
112 maxKinEnergy = 10.0*
TeV;
120 polarAngleLimit = 0.0;
130 secParticles.reserve(5);
135 idxLambda = idxLambdaPrim = 0;
143 secID = fluoID = augerID = biasID = -1;
144 mainSecondaries = 100;
156 <<
" " <<
this <<
" " << theLambdaTable <<
G4endl;
159 delete theLambdaTable;
160 delete theLambdaTablePrim;
170 void G4VEmProcess::Clear()
175 idxLambda = idxLambdaPrim = 0;
192 modelManager->
AddEmModel(order, p, fm, region);
200 G4int n = emModels.size();
201 if(index >= n) {
for(
G4int i=n; i<=
index; ++i) {emModels.push_back(0);} }
210 if(index >= 0 && index <
G4int(emModels.size())) { p = emModels[
index]; }
226 return modelManager->
GetModel(idx, ver);
241 if(pname !=
"deuteron" && pname !=
"triton" &&
242 pname !=
"alpha" && pname !=
"He3" &&
243 pname !=
"alpha+" && pname !=
"helium" &&
244 pname !=
"hydrogen") {
251 G4cout <<
"G4VEmProcess::PreparePhysicsTable() for "
258 if(particle != &part) {
return; }
271 theEnergyOfCrossSectionMax.resize(n, 0.0);
272 theCrossSectionMax.resize(n,
DBL_MAX);
276 for(
G4int i=0; i<numberOfModels; ++i) {
278 if(0 == i) { currentModel = mod; }
287 theCuts = modelManager->
Initialise(particle,secondaryParticle,
294 if(buildLambdaTable && isMaster){
300 if(isMaster && minKinEnergyPrim < maxKinEnergy){
334 G4cout <<
"G4VEmProcess::BuildPhysicsTable() for "
336 <<
" and particle " << num
337 <<
" buildLambdaTable= " << buildLambdaTable
341 if(particle == &part) {
352 }
else if(theLambdaTablePrim) {
357 if(theLambdaTable) { FindLambdaMax(); }
362 for(
G4int i=0; i<numberOfModels; ++i) {
371 if(buildLambdaTable || minKinEnergyPrim < maxKinEnergy) {
380 num ==
"e+" || num ==
"mu+" ||
381 num ==
"mu-" || num ==
"proton"||
382 num ==
"pi+" || num ==
"pi-" ||
383 num ==
"kaon+" || num ==
"kaon-" ||
384 num ==
"alpha" || num ==
"anti_proton" ||
385 num ==
"GenericIon")))
387 PrintInfoProcess(part);
391 G4cout <<
"G4VEmProcess::BuildPhysicsTable() done for "
393 <<
" and particle " << num
400 void G4VEmProcess::BuildLambdaTable()
403 G4cout <<
"G4EmProcess::BuildLambdaTable() for process "
422 if(minKinEnergyPrim < maxKinEnergy) { emax1 = minKinEnergyPrim; }
424 for(
size_t i=0; i<numOfCouples; ++i) {
433 if(buildLambdaTable) {
434 delete (*theLambdaTable)[i];
447 if(emax <= emin) { emax = 2*emin; }
449 if(bin < 3) { bin = 3; }
457 if(minKinEnergyPrim < maxKinEnergy) {
458 delete (*theLambdaTablePrim)[i];
463 G4lrint(nLambdaBins*
G4Log(maxKinEnergy/minKinEnergyPrim)/scale);
464 if(bin < 3) { bin = 3; }
467 bVectorPrim = aVectorPrim;
482 if(buildLambdaTable) { FindLambdaMax(); }
485 G4cout <<
"Lambda table is built for "
498 if(integral) {
G4cout <<
", integral: 1 "; }
499 if(applyCuts) {
G4cout <<
", applyCuts: 1 "; }
501 if(biasFactor != 1.0) {
G4cout <<
" BiasingFactor= " << biasFactor; }
503 if(buildLambdaTable) {
504 size_t length = theLambdaTable->
length();
505 for(
size_t i=0; i<length; ++i) {
508 G4cout <<
" Lambda table from "
520 if(minKinEnergyPrim < maxKinEnergy) {
521 size_t length = theLambdaTablePrim->
length();
522 for(
size_t i=0; i<length; ++i) {
525 G4cout <<
" LambdaPrime table from "
541 G4cout <<
" LambdaTable address= " << theLambdaTable <<
G4endl;
542 if(theLambdaTable) {
G4cout << (*theLambdaTable) <<
G4endl; }
581 if(!currentModel->
IsActive(preStepKinEnergy)) {
591 return biasManager->
GetStepLimit(currentCoupleIndex, previousStepSize);
597 if(preStepKinEnergy < mfpKinEnergy) {
598 if (integral) { ComputeIntegralLambda(preStepKinEnergy); }
599 else { preStepLambda = GetCurrentLambda(preStepKinEnergy); }
602 if(preStepLambda <= 0.0) {
609 if(preStepLambda > 0.0) {
635 G4cout <<
"G4VEmProcess::PostStepGetPhysicalInteractionLength ";
638 <<
" in Material " << currentMaterial->
GetName()
639 <<
" Ekin(MeV)= " << preStepKinEnergy/
MeV
641 G4cout <<
" MeanFreePath = " << currentInteractionLength/
cm <<
"[cm]"
642 <<
" InteractionLength= " << x/
cm <<
"[cm] " <<
G4endl;
679 <<
" E(MeV)= " << finalT/
MeV
680 <<
" preLambda= " << preStepLambda <<
" < "
681 << lx <<
" (postLambda) "
697 weight /= biasFactor;
712 secParticles.clear();
716 (*theCuts)[currentCoupleIndex]);
718 G4int num0 = secParticles.size();
726 currentCoupleIndex, (*theCuts)[currentCoupleIndex],
736 G4int num = secParticles.size();
743 for (
G4int i=0; i<num; ++i) {
744 if (secParticles[i]) {
751 if (e < (*theCutsGamma)[currentCoupleIndex]) { good =
false; }
753 }
else if (p == theElectron) {
754 if (e < (*theCutsElectron)[currentCoupleIndex]) { good =
false; }
756 }
else if (p == thePositron) {
758 e < (*theCutsPositron)[currentCoupleIndex]) {
813 if ( theLambdaTable && part == particle) {
821 <<
" in the directory <" << directory
824 G4cout <<
"Fail to store Physics Table for "
827 <<
" in the directory <" << directory
831 if ( theLambdaTablePrim && part == particle) {
837 G4cout <<
"Physics table prim is stored for "
840 <<
" in the directory <" << directory
843 G4cout <<
"Fail to store Physics Table Prim for "
846 <<
" in the directory <" << directory
860 G4cout <<
"G4VEmProcess::RetrievePhysicsTable() for "
866 if((!buildLambdaTable && minKinEnergyPrim > maxKinEnergy)
867 || particle != part) {
return yes; }
872 if(buildLambdaTable) {
878 G4cout <<
"Lambda table for " << particleName
879 <<
" is Retrieved from <"
884 size_t n = theLambdaTable->
length();
885 for(
size_t i=0; i<
n; ++i) {
886 if((* theLambdaTable)[i]) {
887 (* theLambdaTable)[i]->SetSpline(
true);
893 G4cout <<
"Lambda table for " << particleName <<
" in file <"
894 << filename <<
"> is not exist"
899 if(minKinEnergyPrim < maxKinEnergy) {
905 G4cout <<
"Lambda table prim for " << particleName
906 <<
" is Retrieved from <"
911 size_t n = theLambdaTablePrim->
length();
912 for(
size_t i=0; i<
n; ++i) {
913 if((* theLambdaTablePrim)[i]) {
914 (* theLambdaTablePrim)[i]->SetSpline(
true);
920 G4cout <<
"Lambda table prim for " << particleName <<
" in file <"
921 << filename <<
"> is not exist"
937 DefineMaterial(couple);
939 if(buildLambdaTable && theLambdaTable) {
940 cross = GetCurrentLambda(kineticEnergy);
948 if(cross < 0.0) { cross = 0.0; }
969 if(0.0 < preStepLambda) { x = 1.0/preStepLambda; }
990 void G4VEmProcess::FindLambdaMax()
993 G4cout <<
"### G4VEmProcess::FindLambdaMax: "
997 size_t n = theLambdaTable->
length();
1004 for (i=0; i<
n; ++i) {
1005 pv = (*theLambdaTable)[i];
1011 for (
size_t j=0; j<nb; ++j) {
1020 theEnergyOfCrossSectionMax[i] = emax;
1021 theCrossSectionMax[i] =
smax;
1024 <<
" Max CS at i= " << i <<
" emax(MeV)= " << emax/
MeV
1025 <<
" lambda= " << smax <<
G4endl;
1030 for (i=0; i<
n; ++i) {
1031 pv = (*theLambdaTable)[i];
1033 G4int j = (*theDensityIdx)[i];
1034 theEnergyOfCrossSectionMax[i] = theEnergyOfCrossSectionMax[j];
1035 theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
1068 G4cout <<
"### SetCrossSectionBiasingFactor: for "
1071 <<
" biasFactor= " << f <<
" weightFlag= " << flag
1085 G4cout <<
"### ActivateForcedInteraction: for "
1088 <<
" length(mm)= " << length/
mm
1089 <<
" in G4Region <" << r
1090 <<
"> weightFlag= " << flag
1104 if (0.0 <= factor) {
1113 G4cout <<
"### ActivateSecondaryBiasing: for "
1115 <<
" factor= " << factor
1116 <<
" in G4Region <" << region
1117 <<
"> energyLimit(MeV)= " << energyLimit/
MeV
1127 nLambdaBins =
G4lrint(nLambdaBins*std::log(maxKinEnergy/e)
1128 /std::log(maxKinEnergy/minKinEnergy));
1136 nLambdaBins =
G4lrint(nLambdaBins*std::log(e/minKinEnergy)
1137 /std::log(maxKinEnergy/minKinEnergy));
void ActivateForcedInteraction(G4double length=0.0, const G4String &r="", G4bool flag=true)
G4double condition(const G4ErrorSymMatrix &m)
G4PhysicsTable * LambdaTable() const
const G4VProcess * GetMasterProcess() const
virtual void PrintInfo()=0
G4int GetParentID() const
G4bool SplineFlag() const
G4double GetLambda(G4double &kinEnergy, const G4MaterialCutsCouple *couple)
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4PhysicsTable * LambdaTablePrim() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4VEmModel * SelectModel(G4double &kinEnergy, size_t index)
const std::vector< G4double > * GetDensityFactors()
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
static G4LossTableManager * Instance()
void ActivateForcedInteraction(G4double length=0.0, const G4String &r="")
G4double GetMaxEnergy() const
G4double ApplySecondaryBiasing(std::vector< G4DynamicParticle * > &, const G4Track &track, G4VEmModel *currentModel, G4ParticleChangeForGamma *pParticleChange, G4double &eloss, G4int coupleIdx, G4double tcut, G4double safety=0.0)
void UpdateEmModel(const G4String &, G4double, G4double)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
G4VEmProcess(const G4String &name, G4ProcessType type=fElectromagnetic)
const G4DynamicParticle * GetDynamicParticle() const
G4bool SecondaryBiasingRegion(G4int coupleIdx)
G4double ComputeCrossSectionPerAtom(G4double kineticEnergy, G4double Z, G4double A=0., G4double cut=0.0)
G4VEmModel * EmModel(G4int index=1) const
G4VEmModel * GetModel(G4int, G4bool ver=false)
void DeRegister(G4VEnergyLossProcess *p)
const G4String & GetName() const
const G4ThreeVector & GetPosition() const
G4ParticleChangeForGamma fParticleChange
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
G4TrackStatus GetTrackStatus() const
G4bool ForcedInteractionRegion(G4int coupleIdx)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
virtual void InitialiseProcess(const G4ParticleDefinition *)=0
G4double theNumberOfInteractionLengthLeft
G4bool GetFlag(size_t idx) const
void SetTouchableHandle(const G4TouchableHandle &apValue)
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
void ResetForcedInteraction()
G4double GetParentWeight() const
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void BuildPhysicsTable(const G4ParticleDefinition &)
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
size_t GetVectorLength() const
const G4String & GetParticleSubType() const
void ClearNumberOfInteractionLengthLeft()
void FillLambdaVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4bool startFromNull=true, G4EmTableType t=fRestricted)
G4ProcessManager * GetProcessManager() const
void FillSecondDerivatives()
void UpdateEmModel(const G4String &, G4double, G4double)
const G4String & GetPhysicsTableFileName(const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
const G4String & GetParticleName() const
void ActivateSecondaryBiasing(const G4String ®ion, G4double factor, G4double energyLimit)
G4LossTableBuilder * GetTableBuilder()
void SetWeight(G4double aValue)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetCreatorModelIndex(G4int idx)
void Initialise(const G4ParticleDefinition &part, const G4String &procName, G4int verbose)
void SetHighEnergyLimit(G4double)
void SetCrossSectionBiasingFactor(G4double f, G4bool flag=true)
virtual G4double MinPrimaryEnergy(const G4ParticleDefinition *, const G4Material *)
G4double MeanFreePath(const G4Track &track)
void ProposeWeight(G4double finalWeight)
void SetSecondaryWeightByProcess(G4bool)
void SetEmModel(G4VEmModel *, G4int index=1)
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4bool RetrievePhysicsTable(const G4ParticleDefinition *, const G4String &directory, G4bool ascii)
size_t GetTableSize() const
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=0)
void SetParticle(const G4ParticleDefinition *p)
G4double GetStepLimit(G4int coupleIdx, G4double previousStep)
G4double currentInteractionLength
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
const G4ParticleDefinition * GetParticleDefinition() const
const G4String & GetParticleType() const
void Register(G4VEnergyLossProcess *p)
const G4ParticleDefinition * GetParticleDefinition() const
G4double GetGlobalTime() const
void PreparePhysicsTable(const G4ParticleDefinition &)
G4double Energy(size_t index) const
const G4String & GetProcessName() const
void AddSecondary(G4Track *aSecondary)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4TouchableHandle & GetTouchableHandle() const
G4bool IsActive(G4double kinEnergy)
void DumpModelList(G4int verb)
const G4DataVector * Initialise(const G4ParticleDefinition *part, const G4ParticleDefinition *secPart, G4double minSubRange, G4int verb)
void SetMasterThread(G4bool val)
G4double G4Log(G4double x)
void ActivateSecondaryBiasing(const G4String ®ion, G4double factor, G4double energyLimit)
static G4ProductionCutsTable * GetProductionCutsTable()
G4PhysicsVector * LambdaPhysicsVector(const G4MaterialCutsCouple *)
static G4Positron * Positron()
static G4GenericIon * GenericIon()
G4int NumberOfModels() const
void SetMaxKinEnergy(G4double e)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
void AddEmModel(G4int, G4VEmModel *, const G4Region *region=0)
G4double GetProposedKineticEnergy() const
void InitialiseBaseMaterials(G4PhysicsTable *table)
void SetNumberOfSecondaries(G4int totSecondaries)
G4double GetLocalEnergyDeposit() const
G4StepPoint * GetPostStepPoint() const
G4VParticleChange * pParticleChange
const G4Element * GetCurrentElement() const
G4double GetSafety() const
void PreparePhysicsTable(const G4ParticleDefinition *aParticle, G4VEnergyLossProcess *p, G4bool theMaster)
static G4Electron * Electron()
void SetFluoFlag(G4bool val)
G4bool StorePhysicsTable(const G4String &filename, G4bool ascii=false)
G4TrackStatus GetTrackStatus() const
G4VEmModel * GetModelByIndex(G4int idx=0, G4bool ver=false) const
G4VAtomDeexcitation * AtomDeexcitation()
void SetMinKinEnergy(G4double e)
G4double theInitialNumberOfInteractionLength
static G4int Register(G4String &)
void ProposeTrackStatus(G4TrackStatus status)
void InitializeForPostStep(const G4Track &)
void StartTracking(G4Track *)
static G4bool RetrievePhysicsTable(G4PhysicsTable *physTable, const G4String &fileName, G4bool ascii)
G4double CrossSectionPerVolume(G4double kineticEnergy, const G4MaterialCutsCouple *couple)
G4bool StorePhysicsTable(const G4ParticleDefinition *, const G4String &directory, G4bool ascii=false)
G4int GetProcessSubType() const
void SetPolarAngleLimit(G4double)
G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void SetVerboseLevel(G4int value)
const G4Material * GetMaterial() const
const std::vector< G4int > * GetCoupleIndexes()
const G4Element * GetCurrentElement() const