93 secondaryParticle(nullptr),
94 buildLambdaTable(true),
96 theLambdaTable(nullptr),
97 theLambdaTablePrim(nullptr),
98 theDensityFactor(nullptr),
99 theDensityIdx(nullptr),
102 startFromNull(false),
104 currentModel(nullptr),
106 currentParticle(nullptr),
107 currentCouple(nullptr)
113 minKinEnergy = 0.1*
keV;
114 maxKinEnergy = 100.0*
TeV;
117 actBinning = actSpline = actMinKinEnergy = actMaxKinEnergy =
false;
123 biasFactor = fFactor = 1.0;
130 theCuts = theCutsGamma = theCutsElectron = theCutsPositron =
nullptr;
134 secParticles.reserve(5);
136 baseMaterial = currentMaterial =
nullptr;
138 preStepLambda = preStepKinEnergy = 0.0;
141 idxLambda = idxLambdaPrim = currentCoupleIndex
142 = basedCoupleIndex = 0;
145 biasManager =
nullptr;
150 secID = fluoID = augerID = biasID = -1;
151 mainSecondaries = 100;
170 delete theLambdaTable;
172 if(theLambdaTablePrim) {
174 delete theLambdaTablePrim;
185 void G4VEmProcess::Clear()
187 currentCouple =
nullptr;
190 idxLambda = idxLambdaPrim = 0;
207 modelManager->
AddEmModel(order, p, fm, region);
215 G4int n = emModels.size();
217 for(
G4int i=n; i<=index; ++i) {emModels.push_back(
nullptr);}
227 if(index >= 0 && index <
G4int(emModels.size())) { p = emModels[index]; }
264 return modelManager->
GetModel(idx, ver);
274 if(masterProcess && masterProcess !=
this) { isMaster =
false; }
282 if(pname !=
"deuteron" && pname !=
"triton" &&
283 pname !=
"alpha" && pname !=
"He3" &&
284 pname !=
"alpha+" && pname !=
"helium" &&
285 pname !=
"hydrogen") {
292 G4cout <<
"G4VEmProcess::PreparePhysicsTable() for "
299 if(particle != &part) {
return; }
312 theEnergyOfCrossSectionMax.resize(n, 0.0);
313 theCrossSectionMax.resize(n,
DBL_MAX);
316 if(!actMinKinEnergy) { minKinEnergy = theParameters->
MinKinEnergy(); }
317 if(!actMaxKinEnergy) { maxKinEnergy = theParameters->
MaxKinEnergy(); }
318 if(!actSpline) { splineFlag = theParameters->
Spline(); }
328 for(
G4int i=0; i<numberOfModels; ++i) {
330 if(0 == i) { currentModel = mod; }
339 theCuts = modelManager->
Initialise(particle,secondaryParticle,
346 if(buildLambdaTable && isMaster){
352 if(isMaster && minKinEnergyPrim < maxKinEnergy){
380 if(masterProc && masterProc !=
this) { isMaster =
false; }
384 G4cout <<
"### G4VEmProcess::BuildPhysicsTable() for "
386 <<
" and particle " << num
387 <<
" buildLambdaTable= " << buildLambdaTable
388 <<
" isMaster= " << isMaster
392 if(particle == &part) {
403 }
else if(theLambdaTablePrim) {
408 if(theLambdaTable) { FindLambdaMax(); }
413 for(
G4int i=0; i<numberOfModels; ++i) {
422 if(buildLambdaTable || minKinEnergyPrim < maxKinEnergy) {
431 num ==
"e+" || num ==
"mu+" ||
432 num ==
"mu-" || num ==
"proton"||
433 num ==
"pi+" || num ==
"pi-" ||
434 num ==
"kaon+" || num ==
"kaon-" ||
435 num ==
"alpha" || num ==
"anti_proton" ||
436 num ==
"GenericIon")))
438 PrintInfoProcess(part);
442 G4cout <<
"### G4VEmProcess::BuildPhysicsTable() done for "
444 <<
" and particle " << num
451 void G4VEmProcess::BuildLambdaTable()
454 G4cout <<
"G4EmProcess::BuildLambdaTable() for process "
474 scale =
G4Log(scale);
475 if(actBinning) { nbin =
std::max(nbin, nLambdaBins); }
478 for(
size_t i=0; i<numOfCouples; ++i) {
487 if(buildLambdaTable) {
488 delete (*theLambdaTable)[i];
501 if(emax <= emin) { emax = 2*emin; }
503 if(bin < 3) { bin = 3; }
511 if(minKinEnergyPrim < maxKinEnergy) {
512 delete (*theLambdaTablePrim)[i];
517 if(bin < 3) { bin = 3; }
520 bVectorPrim = aVectorPrim;
535 if(buildLambdaTable) { FindLambdaMax(); }
538 G4cout <<
"Lambda table is built for "
549 G4cout << std::setprecision(6);
552 if(integral) {
G4cout <<
", integral: 1 "; }
553 if(applyCuts) {
G4cout <<
", applyCuts: 1 "; }
555 if(biasFactor != 1.0) {
G4cout <<
" BiasingFactor= " << biasFactor; }
556 G4cout <<
" BuildTable= " << buildLambdaTable;
558 if(buildLambdaTable) {
559 if(particle == &part) {
560 size_t length = theLambdaTable->
length();
561 for(
size_t i=0; i<length; ++i) {
564 G4cout <<
" Lambda table from ";
568 if(emin > minKinEnergy) {
G4cout <<
"threshold "; }
572 <<
", " <<
G4lrint(nbin/std::log10(emax/emin))
573 <<
" bins per decade, spline: "
580 G4cout <<
" Used Lambda table of "
584 if(minKinEnergyPrim < maxKinEnergy) {
585 if(particle == &part) {
586 size_t length = theLambdaTablePrim->
length();
587 for(
size_t i=0; i<length; ++i) {
590 G4cout <<
" LambdaPrime table from "
601 G4cout <<
" Used LambdaPrime table of "
610 G4cout <<
" LambdaTable address= " << theLambdaTable <<
G4endl;
611 if(theLambdaTable && particle == &part) {
650 if(!currentModel->
IsActive(preStepKinEnergy)) {
661 return biasManager->
GetStepLimit(currentCoupleIndex, previousStepSize);
667 if(preStepKinEnergy < mfpKinEnergy) {
668 if (integral) { ComputeIntegralLambda(preStepKinEnergy); }
669 else { preStepLambda = GetCurrentLambda(preStepKinEnergy); }
672 if(preStepLambda <= 0.0) {
679 if(preStepLambda > 0.0) {
746 <<
" E(MeV)= " << finalT/
MeV
747 <<
" preLambda= " << preStepLambda <<
" < "
748 << lx <<
" (postLambda) "
764 weight /= biasFactor;
779 secParticles.clear();
783 (*theCuts)[currentCoupleIndex]);
785 G4int num0 = secParticles.size();
793 currentCoupleIndex, (*theCuts)[currentCoupleIndex],
803 G4int num = secParticles.size();
810 for (
G4int i=0; i<num; ++i) {
811 if (secParticles[i]) {
818 if (e < (*theCutsGamma)[currentCoupleIndex]) { good =
false; }
820 }
else if (p == theElectron) {
821 if (e < (*theCutsElectron)[currentCoupleIndex]) { good =
false; }
823 }
else if (p == thePositron) {
825 e < (*theCutsPositron)[currentCoupleIndex]) {
880 if(masterProc && masterProc !=
this) {
return yes; }
882 if ( theLambdaTable && part == particle) {
890 <<
" in the directory <" << directory
893 G4cout <<
"Fail to store Physics Table for "
896 <<
" in the directory <" << directory
900 if ( theLambdaTablePrim && part == particle) {
906 G4cout <<
"Physics table prim is stored for "
909 <<
" in the directory <" << directory
912 G4cout <<
"Fail to store Physics Table Prim for "
915 <<
" in the directory <" << directory
929 G4cout <<
"G4VEmProcess::RetrievePhysicsTable() for "
935 if((!buildLambdaTable && minKinEnergyPrim > maxKinEnergy)
936 || particle != part) {
return yes; }
941 if(buildLambdaTable) {
947 G4cout <<
"Lambda table for " << particleName
948 <<
" is Retrieved from <"
952 if(theParameters->
Spline()) {
953 size_t n = theLambdaTable->
length();
954 for(
size_t i=0; i<
n; ++i) {
955 if((* theLambdaTable)[i]) {
956 (* theLambdaTable)[i]->SetSpline(
true);
962 G4cout <<
"Lambda table for " << particleName <<
" in file <"
963 << filename <<
"> is not exist"
968 if(minKinEnergyPrim < maxKinEnergy) {
974 G4cout <<
"Lambda table prim for " << particleName
975 <<
" is Retrieved from <"
979 if(theParameters->
Spline()) {
980 size_t n = theLambdaTablePrim->
length();
981 for(
size_t i=0; i<
n; ++i) {
982 if((* theLambdaTablePrim)[i]) {
983 (* theLambdaTablePrim)[i]->SetSpline(
true);
989 G4cout <<
"Lambda table prim for " << particleName <<
" in file <"
990 << filename <<
"> is not exist"
1006 DefineMaterial(couple);
1008 if(buildLambdaTable && theLambdaTable) {
1009 cross = GetCurrentLambda(kineticEnergy);
1017 if(cross < 0.0) { cross = 0.0; }
1038 if(0.0 < preStepLambda) { x = 1.0/preStepLambda; }
1059 void G4VEmProcess::FindLambdaMax()
1062 G4cout <<
"### G4VEmProcess::FindLambdaMax: "
1066 size_t n = theLambdaTable->
length();
1073 for (i=0; i<
n; ++i) {
1074 pv = (*theLambdaTable)[i];
1080 for (
size_t j=0; j<nb; ++j) {
1089 theEnergyOfCrossSectionMax[i] =
emax;
1090 theCrossSectionMax[i] =
smax;
1093 <<
" Max CS at i= " << i <<
" emax(MeV)= " << emax/
MeV
1094 <<
" lambda= " << smax <<
G4endl;
1099 for (i=0; i<
n; ++i) {
1100 pv = (*theLambdaTable)[i];
1102 G4int j = (*theDensityIdx)[i];
1103 theEnergyOfCrossSectionMax[i] = theEnergyOfCrossSectionMax[j];
1104 theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
1137 G4cout <<
"### SetCrossSectionBiasingFactor: for "
1140 <<
" biasFactor= " << f <<
" weightFlag= " << flag
1154 G4cout <<
"### ActivateForcedInteraction: for "
1157 <<
" length(mm)= " << length/
mm
1158 <<
" in G4Region <" << r
1159 <<
"> weightFlag= " << flag
1173 if (0.0 <= factor) {
1182 G4cout <<
"### ActivateSecondaryBiasing: for "
1184 <<
" factor= " << factor
1185 <<
" in G4Region <" << region
1186 <<
"> energyLimit(MeV)= " << energyLimit/
MeV
1196 if(5 < n && n < 10000000) {
1201 PrintWarning(
"SetLambdaBinning", e);
1209 if(1.e-3*
eV < e && e < maxKinEnergy) {
1210 nLambdaBins =
G4lrint(nLambdaBins*
G4Log(maxKinEnergy/e)
1211 /
G4Log(maxKinEnergy/minKinEnergy));
1213 actMinKinEnergy =
true;
1214 }
else { PrintWarning(
"SetMinKinEnergy", e); }
1221 if(minKinEnergy < e && e < 1.e+6*
TeV) {
1222 nLambdaBins =
G4lrint(nLambdaBins*
G4Log(e/minKinEnergy)
1223 /
G4Log(maxKinEnergy/minKinEnergy));
1225 actMaxKinEnergy =
true;
1226 }
else { PrintWarning(
"SetMaxKinEnergy", e); }
1234 e <= theParameters->
MaxKinEnergy()) { minKinEnergyPrim = e; }
1235 else { PrintWarning(
"SetMinKinEnergyPrim", e); }
1242 G4String ss =
"G4VEmProcess::" + tit;
1244 ed <<
"Parameter is out of range: " << val
1245 <<
" it will have no effect!\n" <<
" Process "
void ActivateForcedInteraction(G4double length=0.0, const G4String &r="", G4bool flag=true)
G4int NumberOfRegionModels(size_t index_couple) const
G4double condition(const G4ErrorSymMatrix &m)
G4PhysicsTable * LambdaTable() const
G4int NumberOfBinsPerDecade() const
const G4VProcess * GetMasterProcess() const
virtual void PrintInfo()=0
G4int GetParentID() 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)
G4int NumberOfBins() const
virtual void StartTracking(G4Track *) override
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)
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
G4int WorkerVerbose() const
G4double MaxKinEnergy() const
static G4LossTableManager * Instance()
G4int GetNumberOfRegionModels(size_t couple_index) const
void ActivateForcedInteraction(G4double length=0.0, const G4String &r="")
G4double GetMaxEnergy() const
static constexpr double mm
G4double ApplySecondaryBiasing(std::vector< G4DynamicParticle * > &, const G4Track &track, G4VEmModel *currentModel, G4ParticleChangeForGamma *pParticleChange, G4double &eloss, G4int coupleIdx, G4double tcut, G4double safety=0.0)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
G4VEmModel * GetModel(G4int idx, G4bool ver=false)
G4double HighEnergyLimit() const
G4int GetNumberOfModels() const
G4VEmProcess(const G4String &name, G4ProcessType type=fElectromagnetic)
const G4DynamicParticle * GetDynamicParticle() const
G4bool SecondaryBiasingRegion(G4int coupleIdx)
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
G4double ComputeCrossSectionPerAtom(G4double kineticEnergy, G4double Z, G4double A=0., G4double cut=0.0)
G4VEmModel * EmModel(G4int index=1) const
void DeRegister(G4VEnergyLossProcess *p)
void UpdateEmModel(const G4String &model_name, G4double emin, G4double emax)
G4double MscThetaLimit() const
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &) override
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 G4bool StorePhysicsTable(const G4ParticleDefinition *, const G4String &directory, G4bool ascii=false) override
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
size_t GetVectorLength() const
const G4String & GetParticleSubType() const
void ClearNumberOfInteractionLengthLeft()
void FillLambdaVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4bool startFromNull=true, G4EmTableType t=fRestricted)
virtual G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
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)
virtual G4bool RetrievePhysicsTable(const G4ParticleDefinition *, const G4String &directory, G4bool ascii) override
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetCreatorModelIndex(G4int idx)
virtual void PreparePhysicsTable(const G4ParticleDefinition &) override
void Initialise(const G4ParticleDefinition &part, const G4String &procName, G4int verbose)
void SetHighEnergyLimit(G4double)
void SetMinKinEnergyPrim(G4double e)
void SetCrossSectionBiasingFactor(G4double f, G4bool flag=true)
static constexpr double TeV
virtual G4double MinPrimaryEnergy(const G4ParticleDefinition *, const G4Material *)
G4double MeanFreePath(const G4Track &track)
void SetLambdaBinning(G4int nbins)
void ProposeWeight(G4double finalWeight)
void SetSecondaryWeightByProcess(G4bool)
G4double LambdaFactor() const
void SetEmModel(G4VEmModel *, G4int index=1)
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4double GetKineticEnergy() const
G4VEmModel * GetRegionModel(G4int idx, size_t index_couple)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
size_t GetTableSize() const
G4VEmModel * GetRegionModel(G4int idx, size_t couple_index) const
void SetParticle(const G4ParticleDefinition *p)
G4double GetStepLimit(G4int coupleIdx, G4double previousStep)
G4double currentInteractionLength
const G4ParticleDefinition * GetParticleDefinition() const
static constexpr double eV
const G4String & GetParticleType() const
void Register(G4VEnergyLossProcess *p)
G4double MinKinEnergy() const
const G4ParticleDefinition * GetParticleDefinition() const
G4double GetGlobalTime() const
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)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
const G4DataVector * Initialise(const G4ParticleDefinition *part, const G4ParticleDefinition *secPart, G4double minSubRange, G4int verb)
static const G4double emax
void SetMasterThread(G4bool val)
G4double G4Log(G4double x)
void ActivateSecondaryBiasing(const G4String ®ion, G4double factor, G4double energyLimit)
void AddEmModel(G4int, G4VEmModel *, const G4Region *region=nullptr)
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
virtual void ProcessDescription(std::ostream &outFile) const
G4double GetProposedKineticEnergy() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4ProcessManager * GetProcessManager() const
void InitialiseBaseMaterials(G4PhysicsTable *table)
void SetNumberOfSecondaries(G4int totSecondaries)
static G4int Register(const G4String &)
G4double GetLocalEnergyDeposit() const
G4StepPoint * GetPostStepPoint() const
static G4EmParameters * Instance()
G4VParticleChange * pParticleChange
const G4Element * GetCurrentElement() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double GeV
G4double GetSafety() const
void PreparePhysicsTable(const G4ParticleDefinition *aParticle, G4VEnergyLossProcess *p, G4bool theMaster)
static G4Electron * Electron()
void SetFluoFlag(G4bool val)
static constexpr double MeV
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
void ProposeTrackStatus(G4TrackStatus status)
void DefineRegParamForEM(G4VEmProcess *) const
void InitializeForPostStep(const G4Track &)
static G4bool RetrievePhysicsTable(G4PhysicsTable *physTable, const G4String &fileName, G4bool ascii)
static constexpr double keV
G4double CrossSectionPerVolume(G4double kineticEnergy, const G4MaterialCutsCouple *couple)
G4int GetProcessSubType() const
void SetPolarAngleLimit(G4double)
void SetVerboseLevel(G4int value)
const G4Material * GetMaterial() const
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
const std::vector< G4int > * GetCoupleIndexes()
const G4Element * GetCurrentElement() const