G4EmBiasingManager Class Reference

#include <G4EmBiasingManager.hh>

Collaboration diagram for G4EmBiasingManager:

Public Member Functions
	G4EmBiasingManager ()
	~G4EmBiasingManager ()
void	Initialise (const G4ParticleDefinition &part, const G4String &procName, G4int verbose)
void	ActivateForcedInteraction (G4double length=0.0, const G4String &r="")
void	ActivateSecondaryBiasing (const G4String &region, G4double factor, G4double energyLimit)
G4double	GetStepLimit (G4int coupleIdx, G4double previousStep)
G4double	ApplySecondaryBiasing (std::vector< G4DynamicParticle *> &, const G4Track &track, G4VEmModel *currentModel, G4ParticleChangeForGamma *pParticleChange, G4double &eloss, G4int coupleIdx, G4double tcut, G4double safety=0.0)
G4double	ApplySecondaryBiasing (std::vector< G4DynamicParticle *> &, const G4Track &track, G4VEmModel *currentModel, G4ParticleChangeForLoss *pParticleChange, G4double &eloss, G4int coupleIdx, G4double tcut, G4double safety=0.0)
G4double	ApplySecondaryBiasing (std::vector< G4Track *> &, G4int coupleIdx)
G4bool	SecondaryBiasingRegion (G4int coupleIdx)
G4bool	ForcedInteractionRegion (G4int coupleIdx)
void	ResetForcedInteraction ()

Private Member Functions
void	ApplyRangeCut (std::vector< G4DynamicParticle *> &vd, const G4Track &track, G4double &eloss, G4double safety)
G4double	ApplySplitting (std::vector< G4DynamicParticle *> &vd, const G4Track &track, G4VEmModel *currentModel, G4int index, G4double tcut)
G4double	ApplyRussianRoulette (std::vector< G4DynamicParticle *> &vd, G4int index)
	G4EmBiasingManager (G4EmBiasingManager &)
G4EmBiasingManager &	operator= (const G4EmBiasingManager &right)

Private Attributes
G4int	nForcedRegions
G4int	nSecBiasedRegions
std::vector< const G4Region * >	forcedRegions
std::vector< G4double >	lengthForRegion
std::vector< const G4Region * >	secBiasedRegions
std::vector< G4double >	secBiasedWeight
std::vector< G4double >	secBiasedEnegryLimit
std::vector< G4int >	nBremSplitting
std::vector< G4int >	idxForcedCouple
std::vector< G4int >	idxSecBiasedCouple
std::vector< G4DynamicParticle * >	tmpSecondaries
G4VEnergyLossProcess *	eIonisation
const G4ParticleDefinition *	theElectron
G4double	fSafetyMin
G4double	currentStepLimit
G4bool	startTracking

Detailed Description

Definition at line 67 of file G4EmBiasingManager.hh.

Constructor & Destructor Documentation

G4EmBiasingManager() [1/2]

G4EmBiasingManager::G4EmBiasingManager

(

)

Definition at line 65 of file G4EmBiasingManager.cc.

  : nForcedRegions(0),nSecBiasedRegions(0),eIonisation(nullptr),
    currentStepLimit(0.0),startTracking(true)
{
  fSafetyMin = 1.e-6*mm;
  theElectron = G4Electron::Electron();
}

Here is the call graph for this function:

~G4EmBiasingManager()

G4EmBiasingManager::~G4EmBiasingManager

(

)

Definition at line 75 of file G4EmBiasingManager.cc.

{}

G4EmBiasingManager() [2/2]

G4EmBiasingManager::G4EmBiasingManager ( G4EmBiasingManager &  )

private

Member Function Documentation

ActivateForcedInteraction()

void G4EmBiasingManager::ActivateForcedInteraction	(	G4double	length = 0.0,
		const G4String &	r = ""
	)

Definition at line 146 of file G4EmBiasingManager.cc.

{
  G4RegionStore* regionStore = G4RegionStore::GetInstance();
  G4String name = rname;
  if(name == "" || name == "world" || name == "World") {
    name = "DefaultRegionForTheWorld";
  }
  const G4Region* reg = regionStore->GetRegion(name, false);
  if(!reg) { 
    G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
           << " G4Region <"
           << rname << "> is unknown" << G4endl;
    return; 
  }

  // the region is in the list
  if (0 < nForcedRegions) {
    for (G4int i=0; i<nForcedRegions; ++i) {
      if (reg == forcedRegions[i]) {
        lengthForRegion[i] = val; 
        return;
      }
    }
  }
  if(val < 0.0) { 
    G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
           << val << " < 0.0, so no activation for the G4Region <"
           << rname << ">" << G4endl;
    return; 
  }

  // new region 
  forcedRegions.push_back(reg);
  lengthForRegion.push_back(val);
  ++nForcedRegions;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ActivateSecondaryBiasing()

void G4EmBiasingManager::ActivateSecondaryBiasing	(	const G4String &	region,
		G4double	factor,
		G4double	energyLimit
	)

Definition at line 187 of file G4EmBiasingManager.cc.

{
  //G4cout << "G4EmBiasingManager::ActivateSecondaryBiasing: "
  //         << rname << " F= " << factor << " E(MeV)= " << energyLimit/MeV
  //         << G4endl; 
  G4RegionStore* regionStore = G4RegionStore::GetInstance();
  G4String name = rname;
  if(name == "" || name == "world" || name == "World") {
    name = "DefaultRegionForTheWorld";
  }
  const G4Region* reg = regionStore->GetRegion(name, false);
  if(!reg) { 
    G4cout << "### G4EmBiasingManager::ActivateBremsstrahlungSplitting "
           << "WARNING: G4Region <"
           << rname << "> is unknown" << G4endl;
    return; 
  }

  // Range cut
  G4int nsplit = 0;
  G4double w = factor;

  // splitting
  if(factor >= 1.0) {
    nsplit = G4lrint(factor);
    w = 1.0/G4double(nsplit);

    // Russian roulette 
  } else if(0.0 < factor) { 
    nsplit = 1;
    w = 1.0/factor;
  }

  // the region is in the list - overwrite parameters
  if (0 < nSecBiasedRegions) {
    for (G4int i=0; i<nSecBiasedRegions; ++i) {
      if (reg == secBiasedRegions[i]) {
        secBiasedWeight[i] = w;
        nBremSplitting[i]  = nsplit; 
        secBiasedEnegryLimit[i] = energyLimit;
        return;
      }
    }
  }
  /*
    G4cout << "### G4EmBiasingManager::ActivateSecondaryBiasing: "
           << " nsplit= " << nsplit << " for the G4Region <"
           << rname << ">" << G4endl; 
  */

  // new region 
  secBiasedRegions.push_back(reg);
  secBiasedWeight.push_back(w);
  nBremSplitting.push_back(nsplit);
  secBiasedEnegryLimit.push_back(energyLimit);
  ++nSecBiasedRegions;
  //G4cout << "nSecBiasedRegions= " << nSecBiasedRegions << G4endl;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyRangeCut()

void G4EmBiasingManager::ApplyRangeCut ( std::vector< G4DynamicParticle *> &  vd, const G4Track &  track, G4double &  eloss, G4double  safety  )

private

Definition at line 402 of file G4EmBiasingManager.cc.

{
  size_t n = vd.size();
  if(!eIonisation) { 
    eIonisation = 
      G4LossTableManager::Instance()->GetEnergyLossProcess(theElectron);
  }
  if(eIonisation) { 
    for(size_t k=0; k<n; ++k) {
      const G4DynamicParticle* dp = vd[k];
      if(dp->GetDefinition() == theElectron) {
        G4double e = dp->GetKineticEnergy();
        if(eIonisation->GetRangeForLoss(e, track.GetMaterialCutsCouple()) 
           < safety) {
          eloss += e;
          delete dp;
          vd[k] = 0;
        }
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyRussianRoulette()

G4double G4EmBiasingManager::ApplyRussianRoulette ( std::vector< G4DynamicParticle *> &  vd, G4int  index  )

inlineprivate

Definition at line 194 of file G4EmBiasingManager.hh.

{
  size_t n = vd.size();
  G4double weight = secBiasedWeight[index];
  for(size_t k=0; k<n; ++k) {
    if(G4UniformRand()*weight > 1.0) {
      const G4DynamicParticle* dp = vd[k];
      delete dp;
      vd[k] = 0;
    }
  }
  return weight;
}

Here is the caller graph for this function:

ApplySecondaryBiasing() [1/3]

G4double G4EmBiasingManager::ApplySecondaryBiasing	(	std::vector< G4DynamicParticle *> &	vd,
		const G4Track &	track,
		G4VEmModel *	currentModel,
		G4ParticleChangeForGamma *	pParticleChange,
		G4double &	eloss,
		G4int	coupleIdx,
		G4double	tcut,
		G4double	safety = 0.0
	)

Definition at line 320 of file G4EmBiasingManager.cc.

{
  G4int index = idxSecBiasedCouple[coupleIdx];
  G4double weight = 1.0;
  if(0 <= index) {
    size_t n = vd.size();

    // the check cannot be applied per secondary particle
    // because weight correction is common, so the first
    // secondary is checked
    if(0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {

      G4int nsplit = nBremSplitting[index];

      // Range cut
      if(0 == nsplit) { 
        if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }

        // Russian Roulette
      } else if(1 == nsplit) { 
        weight = ApplyRussianRoulette(vd, index);

        // Splitting
      } else {
        G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
        G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();

        weight = ApplySplitting(vd, track, currentModel, index, tcut);

        pPartChange->SetProposedKineticEnergy(tmpEnergy);
        pPartChange->ProposeMomentumDirection(tmpMomDir);
      }
    }
  }
  return weight;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ApplySecondaryBiasing() [2/3]

G4double G4EmBiasingManager::ApplySecondaryBiasing	(	std::vector< G4DynamicParticle *> &	vd,
		const G4Track &	track,
		G4VEmModel *	currentModel,
		G4ParticleChangeForLoss *	pParticleChange,
		G4double &	eloss,
		G4int	coupleIdx,
		G4double	tcut,
		G4double	safety = 0.0
	)

Definition at line 272 of file G4EmBiasingManager.cc.

{
  G4int index = idxSecBiasedCouple[coupleIdx];
  G4double weight = 1.0;
  if(0 <= index) {
    size_t n = vd.size();

    // the check cannot be applied per secondary particle
    // because weight correction is common, so the first
    // secondary is checked
    if(0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {

      G4int nsplit = nBremSplitting[index];

      // Range cut
      if(0 == nsplit) { 
        if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }

        // Russian Roulette
      } else if(1 == nsplit) { 
        weight = ApplyRussianRoulette(vd, index);

        // Splitting
      } else {
        G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
        G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();

        weight = ApplySplitting(vd, track, currentModel, index, tcut);

        pPartChange->SetProposedKineticEnergy(tmpEnergy);
        pPartChange->ProposeMomentumDirection(tmpMomDir);
      }
    }
  }
  return weight;
}

Here is the call graph for this function:

ApplySecondaryBiasing() [3/3]

G4double G4EmBiasingManager::ApplySecondaryBiasing	(	std::vector< G4Track *> &	track,
		G4int	coupleIdx
	)

Definition at line 368 of file G4EmBiasingManager.cc.

{
  G4int index = idxSecBiasedCouple[coupleIdx];
  G4double weight = 1.0;
  if(0 <= index) {
    size_t n = track.size();

    // the check cannot be applied per secondary particle
    // because weight correction is common, so the first
    // secondary is checked
    if(0 < n && track[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {

      G4int nsplit = nBremSplitting[index];

        // Russian Roulette only
      if(1 == nsplit) { 
        weight = secBiasedWeight[index];
        for(size_t k=0; k<n; ++k) {
          if(G4UniformRand()*weight > 1.0) {
            const G4Track* t = track[k];
            delete t;
            track[k] = 0;
          }
        }
      }
    }
  }
  return weight;
}

ApplySplitting()

G4double G4EmBiasingManager::ApplySplitting ( std::vector< G4DynamicParticle *> &  vd, const G4Track &  track, G4VEmModel *  currentModel, G4int  index, G4double  tcut  )

private

Definition at line 430 of file G4EmBiasingManager.cc.

{
  // method is applied only if 1 secondary created PostStep 
  // in the case of many secodndaries there is a contrudition
  G4double weight = 1.0;
  size_t n = vd.size();
  G4double w = secBiasedWeight[index];

  if(1 != n || 1.0 <= w) { return weight; }

  G4double trackWeight = track.GetWeight();
  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();

  G4int nsplit = nBremSplitting[index];

  // double splitting is supressed 
  if(1 < nsplit && trackWeight>w) {

    weight = w;
    // start from 1, because already one secondary created
    if(nsplit > (G4int)tmpSecondaries.size()) { 
      tmpSecondaries.reserve(nsplit);
    }
    const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
    for(G4int k=1; k<nsplit; ++k) {  
      tmpSecondaries.clear();
      currentModel->SampleSecondaries(&tmpSecondaries, couple, dynParticle, 
                                      tcut);
      for (size_t kk=0; kk<tmpSecondaries.size(); ++kk) {
        vd.push_back(tmpSecondaries[kk]);
      }
    }
  }
  return weight;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ForcedInteractionRegion()

G4bool G4EmBiasingManager::ForcedInteractionRegion ( G4int  coupleIdx )

inline

Definition at line 177 of file G4EmBiasingManager.hh.

{
  G4bool res = false;
  if(nForcedRegions > 0) {
    if(idxForcedCouple[coupleIdx] >= 0) { res = true; }
  }
  return res;
}

Here is the caller graph for this function:

GetStepLimit()

G4double G4EmBiasingManager::GetStepLimit	(	G4int	coupleIdx,
		G4double	previousStep
	)

Definition at line 250 of file G4EmBiasingManager.cc.

{
  if(startTracking) {
    startTracking = false;
    G4int i = idxForcedCouple[coupleIdx];
    if(i < 0) {
      currentStepLimit = DBL_MAX;
    } else {
      currentStepLimit = lengthForRegion[i];
      if(currentStepLimit > 0.0) { currentStepLimit *= G4UniformRand(); }
    }
  } else {
    currentStepLimit -= previousStep;
  }
  if(currentStepLimit < 0.0) { currentStepLimit = 0.0; }
  return currentStepLimit;
}

Here is the caller graph for this function:

Initialise()

void G4EmBiasingManager::Initialise	(	const G4ParticleDefinition &	part,
		const G4String &	procName,
		G4int	verbose
	)

Definition at line 80 of file G4EmBiasingManager.cc.

{
  //G4cout << "G4EmBiasingManager::Initialise for "
  //         << part.GetParticleName()
  //         << " and " << procName << G4endl;
  const G4ProductionCutsTable* theCoupleTable=
    G4ProductionCutsTable::GetProductionCutsTable();
  size_t numOfCouples = theCoupleTable->GetTableSize();

  if(0 < nForcedRegions) { idxForcedCouple.resize(numOfCouples, -1); }
  if(0 < nSecBiasedRegions) { idxSecBiasedCouple.resize(numOfCouples, -1); }

  // Deexcitation
  for (size_t j=0; j<numOfCouples; ++j) {
    const G4MaterialCutsCouple* couple =
      theCoupleTable->GetMaterialCutsCouple(j);
    const G4ProductionCuts* pcuts = couple->GetProductionCuts();
    if(0 <  nForcedRegions) {
      for(G4int i=0; i<nForcedRegions; ++i) {
        if(forcedRegions[i]) {
          if(pcuts == forcedRegions[i]->GetProductionCuts()) { 
            idxForcedCouple[j] = i;
            break; 
          }
        }
      }
    }
    if(0 < nSecBiasedRegions) { 
      for(G4int i=0; i<nSecBiasedRegions; ++i) {
        if(secBiasedRegions[i]) {
          if(pcuts == secBiasedRegions[i]->GetProductionCuts()) { 
            idxSecBiasedCouple[j] = i;
            break; 
          }
        }
      }
    }
  }
  if (nForcedRegions > 0 && 0 < verbose) {
    G4cout << " Forced Interaction is activated for "
           << part.GetParticleName() << " and " 
           << procName 
           << " inside G4Regions: " << G4endl;
    for (G4int i=0; i<nForcedRegions; ++i) {
      const G4Region* r = forcedRegions[i];
      if(r) { G4cout << "           " << r->GetName() << G4endl; }
    }
  }
  if (nSecBiasedRegions > 0 && 0 < verbose) {
    G4cout << " Secondary biasing is activated for " 
           << part.GetParticleName() << " and " 
           << procName 
           << " inside G4Regions: " << G4endl;
    for (G4int i=0; i<nSecBiasedRegions; ++i) {
      const G4Region* r = secBiasedRegions[i];
      if(r) { 
        G4cout << "           " << r->GetName() 
               << "  BiasingWeight= " << secBiasedWeight[i] << G4endl; 
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

operator=()

G4EmBiasingManager& G4EmBiasingManager::operator= ( const G4EmBiasingManager &  right )

private

ResetForcedInteraction()

void G4EmBiasingManager::ResetForcedInteraction ( )

inline

Definition at line 186 of file G4EmBiasingManager.hh.

{
  startTracking = true;
}

Here is the caller graph for this function:

SecondaryBiasingRegion()

G4bool G4EmBiasingManager::SecondaryBiasingRegion ( G4int  coupleIdx )

inline

Definition at line 168 of file G4EmBiasingManager.hh.

{
  G4bool res = false;
  if(nSecBiasedRegions > 0) {
    if(idxSecBiasedCouple[coupleIdx] >= 0) { res = true; }
  }
  return res;
}

Here is the caller graph for this function:

Member Data Documentation

currentStepLimit

G4double G4EmBiasingManager::currentStepLimit

private

Definition at line 163 of file G4EmBiasingManager.hh.

eIonisation

G4VEnergyLossProcess* G4EmBiasingManager::eIonisation

private

Definition at line 158 of file G4EmBiasingManager.hh.

forcedRegions

std::vector<const G4Region*> G4EmBiasingManager::forcedRegions

private

Definition at line 146 of file G4EmBiasingManager.hh.

fSafetyMin

G4double G4EmBiasingManager::fSafetyMin

private

Definition at line 162 of file G4EmBiasingManager.hh.

idxForcedCouple

std::vector<G4int> G4EmBiasingManager::idxForcedCouple

private

Definition at line 153 of file G4EmBiasingManager.hh.

idxSecBiasedCouple

std::vector<G4int> G4EmBiasingManager::idxSecBiasedCouple

private

Definition at line 154 of file G4EmBiasingManager.hh.

lengthForRegion

std::vector<G4double> G4EmBiasingManager::lengthForRegion

private

Definition at line 147 of file G4EmBiasingManager.hh.

nBremSplitting

std::vector<G4int> G4EmBiasingManager::nBremSplitting

private

Definition at line 151 of file G4EmBiasingManager.hh.

nForcedRegions

G4int G4EmBiasingManager::nForcedRegions

private

Definition at line 144 of file G4EmBiasingManager.hh.

nSecBiasedRegions

G4int G4EmBiasingManager::nSecBiasedRegions

private

Definition at line 145 of file G4EmBiasingManager.hh.

secBiasedEnegryLimit

std::vector<G4double> G4EmBiasingManager::secBiasedEnegryLimit

private

Definition at line 150 of file G4EmBiasingManager.hh.

secBiasedRegions

std::vector<const G4Region*> G4EmBiasingManager::secBiasedRegions

private

Definition at line 148 of file G4EmBiasingManager.hh.

secBiasedWeight

std::vector<G4double> G4EmBiasingManager::secBiasedWeight

private

Definition at line 149 of file G4EmBiasingManager.hh.

startTracking

G4bool G4EmBiasingManager::startTracking

private

Definition at line 164 of file G4EmBiasingManager.hh.

theElectron

const G4ParticleDefinition* G4EmBiasingManager::theElectron

private

Definition at line 160 of file G4EmBiasingManager.hh.

tmpSecondaries

std::vector<G4DynamicParticle*> G4EmBiasingManager::tmpSecondaries

private

Definition at line 156 of file G4EmBiasingManager.hh.

---

The documentation for this class was generated from the following files:

• G4EmBiasingManager.hh
• G4EmBiasingManager.cc