G4IonParametrisedLossModel Class Reference

#include <G4IonParametrisedLossModel.hh>

Inheritance diagram for G4IonParametrisedLossModel:

Public Member Functions
	G4IonParametrisedLossModel (const G4ParticleDefinition *particle=0, const G4String &name="ParamICRU73")
virtual	~G4IonParametrisedLossModel ()
virtual void	Initialise (const G4ParticleDefinition *, const G4DataVector &)
virtual G4double	MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
virtual G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
virtual G4double	CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double, G4double, G4double)
virtual G4double	ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double, G4double)
G4double	ComputeLossForStep (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double, G4double)
G4double	DeltaRayMeanEnergyTransferRate (const G4Material *, const G4ParticleDefinition *, G4double, G4double)
virtual void	SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double)
virtual G4double	GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double)
virtual G4double	GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double)
virtual void	CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &, G4double &, G4double)
G4bool	AddDEDXTable (const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=0)
G4bool	RemoveDEDXTable (const G4String &name)
void	DeactivateICRU73Scaling ()
LossTableList::iterator	IsApplicable (const G4ParticleDefinition *, const G4Material *)
void	PrintDEDXTable (const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
void	PrintDEDXTableHandlers (const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
void	SetEnergyLossLimit (G4double ionEnergyLossLimit)
Public Member Functions inherited from G4VEmModel
	G4VEmModel (const G4String &nam)
virtual	~G4VEmModel ()
virtual void	InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
virtual void	InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
virtual void	InitialiseForElement (const G4ParticleDefinition *, G4int Z)
virtual G4double	ChargeSquareRatio (const G4Track &)
virtual void	StartTracking (G4Track *)
virtual G4double	Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
virtual G4double	MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
virtual void	SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual void	DefineForRegion (const G4Region *)
void	InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
std::vector < G4EmElementSelector * > *	GetElementSelectors ()
void	SetElementSelectors (std::vector< G4EmElementSelector * > *)
G4double	ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
G4double	CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double	ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4int	SelectIsotopeNumber (const G4Element *)
const G4Element *	SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element *	SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4int	SelectRandomAtomNumber (const G4Material *)
void	SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
void	SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
G4ElementData *	GetElementData ()
G4PhysicsTable *	GetCrossSectionTable ()
G4VEmFluctuationModel *	GetModelOfFluctuations ()
G4VEmAngularDistribution *	GetAngularDistribution ()
void	SetAngularDistribution (G4VEmAngularDistribution *)
G4double	HighEnergyLimit () const
G4double	LowEnergyLimit () const
G4double	HighEnergyActivationLimit () const
G4double	LowEnergyActivationLimit () const
G4double	PolarAngleLimit () const
G4double	SecondaryThreshold () const
G4bool	LPMFlag () const
G4bool	DeexcitationFlag () const
G4bool	ForceBuildTableFlag () const
G4bool	UseAngularGeneratorFlag () const
void	SetAngularGeneratorFlag (G4bool)
void	SetHighEnergyLimit (G4double)
void	SetLowEnergyLimit (G4double)
void	SetActivationHighEnergyLimit (G4double)
void	SetActivationLowEnergyLimit (G4double)
G4bool	IsActive (G4double kinEnergy)
void	SetPolarAngleLimit (G4double)
void	SetSecondaryThreshold (G4double)
void	SetLPMFlag (G4bool val)
void	SetDeexcitationFlag (G4bool val)
void	SetForceBuildTable (G4bool val)
void	SetMasterThread (G4bool val)
G4bool	IsMaster () const
G4double	MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
const G4String &	GetName () const
void	SetCurrentCouple (const G4MaterialCutsCouple *)
const G4Element *	GetCurrentElement () const

Protected Member Functions
virtual G4double	MaxSecondaryEnergy (const G4ParticleDefinition *, G4double)
Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss *	GetParticleChangeForLoss ()
G4ParticleChangeForGamma *	GetParticleChangeForGamma ()
const G4MaterialCutsCouple *	CurrentCouple () const
void	SetCurrentElement (const G4Element *)

Additional Inherited Members
Protected Attributes inherited from G4VEmModel
G4ElementData *	fElementData
G4VParticleChange *	pParticleChange
G4PhysicsTable *	xSectionTable
const std::vector< G4double > *	theDensityFactor
const std::vector< G4int > *	theDensityIdx
size_t	idxTable

Detailed Description

Definition at line 93 of file G4IonParametrisedLossModel.hh.

Constructor & Destructor Documentation

G4IonParametrisedLossModel::G4IonParametrisedLossModel	(	const G4ParticleDefinition *	particle = 0,
		const G4String &	name = "ParamICRU73"
	)

Definition at line 103 of file G4IonParametrisedLossModel.cc.

References AddDEDXTable(), G4GenericIon::Definition(), G4VEmModel::HighEnergyLimit(), G4LossTableManager::Instance(), python.hepunit::MeV, G4VEmModel::SetHighEnergyLimit(), and python.hepunit::TeV.

  : G4VEmModel(nam),
    braggIonModel(0),
    betheBlochModel(0),
    nmbBins(90),
    nmbSubBins(100),
    particleChangeLoss(0),
    corrFactor(1.0),
    energyLossLimit(0.01),
    cutEnergies(0) 
{
  genericIon = G4GenericIon::Definition();
  genericIonPDGMass = genericIon -> GetPDGMass();
  corrections = G4LossTableManager::Instance() -> EmCorrections();
 
  // The upper limit of the current model is set to 100 TeV
  SetHighEnergyLimit(100.0 * TeV);

  // The Bragg ion and Bethe Bloch models are instantiated
  braggIonModel = new G4BraggIonModel();
  betheBlochModel = new G4BetheBlochModel();

  // By default ICRU 73 stopping power tables are loaded:
  AddDEDXTable("ICRU73",
           new G4IonStoppingData("ion_stopping_data/icru73"),
           new G4IonDEDXScalingICRU73());

  // The boundaries for the range tables are set
  lowerEnergyEdgeIntegr = 0.025 * MeV;
  upperEnergyEdgeIntegr = betheBlochModel -> HighEnergyLimit();

  // Cache parameters are set
  cacheParticle = 0;
  cacheMass = 0;
  cacheElecMassRatio = 0;
  cacheChargeSquare = 0;

  // Cache parameters are set
  rangeCacheParticle = 0; 
  rangeCacheMatCutsCouple = 0; 
  rangeCacheEnergyRange = 0; 
  rangeCacheRangeEnergy = 0;

  // Cache parameters are set
  dedxCacheParticle = 0;
  dedxCacheMaterial = 0;
  dedxCacheEnergyCut = 0;
  dedxCacheIter = lossTableList.end();
  dedxCacheTransitionEnergy = 0.0;  
  dedxCacheTransitionFactor = 0.0;
  dedxCacheGenIonMassRatio = 0.0;
}

G4IonParametrisedLossModel::~G4IonParametrisedLossModel ( )

virtual

Definition at line 160 of file G4IonParametrisedLossModel.cc.

                                                        {

  // dE/dx table objects are deleted and the container is cleared
  LossTableList::iterator iterTables = lossTableList.begin();
  LossTableList::iterator iterTables_end = lossTableList.end();

  for(;iterTables != iterTables_end; iterTables++) delete *iterTables;
  lossTableList.clear();
}

Member Function Documentation

G4bool G4IonParametrisedLossModel::AddDEDXTable	(	const G4String &	name,
		G4VIonDEDXTable *	table,
		G4VIonDEDXScalingAlgorithm *	algorithm = 0
	)

Definition at line 1226 of file G4IonParametrisedLossModel.cc.

References G4cerr, G4endl, and G4VEmModel::GetName().

Referenced by DeactivateICRU73Scaling(), and G4IonParametrisedLossModel().

                                                                       {

  if(table == 0) {
     G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
            << " add table: Invalid pointer."
            << G4endl;      

     return false;
  }

  // Checking uniqueness of name
  LossTableList::iterator iter = lossTableList.begin();
  LossTableList::iterator iter_end = lossTableList.end();
  
  for(;iter != iter_end; iter++) {
     G4String tableName = (*iter) -> GetName();

     if(tableName == nam) { 
        G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
               << " add table: Name already exists."      
               << G4endl;

        return false;
     }
  }

  G4VIonDEDXScalingAlgorithm* scalingAlgorithm = algorithm;
  if(scalingAlgorithm == 0) 
     scalingAlgorithm = new G4VIonDEDXScalingAlgorithm; 
 
  G4IonDEDXHandler* handler = 
                      new G4IonDEDXHandler(table, scalingAlgorithm, nam); 

  lossTableList.push_front(handler);

  return true;
}

G4double G4IonParametrisedLossModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  particle, G4double  kineticEnergy, G4double  atomicNumber, G4double  , G4double  cutEnergy, G4double  maxKinEnergy  )

virtual

Reimplemented from G4VEmModel.

Definition at line 352 of file G4IonParametrisedLossModel.cc.

References energy(), G4cout, G4endl, MaxSecondaryEnergy(), python.hepunit::MeV, G4INCL::Math::min(), right, and python.hepunit::twopi_mc2_rcl2.

Referenced by CrossSectionPerVolume().

                                                          {

  // ############## Cross section per atom  ################################
  // Function computes ionization cross section per atom
  //
  // See Geant4 physics reference manual (version 9.1), section 9.1.3
  // 
  // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
  //       B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
  //
  // (Implementation adapted from G4BraggIonModel)

  G4double crosssection     = 0.0;
  G4double tmax      = MaxSecondaryEnergy(particle, kineticEnergy);
  G4double maxEnergy = std::min(tmax, maxKinEnergy);

  if(cutEnergy < tmax) {

     G4double energy  = kineticEnergy + cacheMass;
     G4double betaSquared  = kineticEnergy * 
                                     (energy + cacheMass) / (energy * energy);

     crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy - 
                         betaSquared * std::log(maxEnergy / cutEnergy) / tmax;

     crosssection *= twopi_mc2_rcl2 * cacheChargeSquare / betaSquared;
  }

#ifdef PRINT_DEBUG_CS
  G4cout << "########################################################"
         << G4endl
         << "# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
         << G4endl
         << "# particle =" << particle -> GetParticleName() 
         <<  G4endl
         << "# cut(MeV) = " << cutEnergy/MeV  
         << G4endl;

  G4cout << "#"
         << std::setw(13) << std::right << "E(MeV)"
         << std::setw(14) << "CS(um)"
         << std::setw(14) << "E_max_sec(MeV)"
         << G4endl
         << "# ------------------------------------------------------"
         << G4endl;

  G4cout << std::setw(14) << std::right << kineticEnergy / MeV
         << std::setw(14) << crosssection / (um * um)
         << std::setw(14) << tmax / MeV
         << G4endl;
#endif

  crosssection *= atomicNumber;

  return crosssection;
}

G4double G4IonParametrisedLossModel::ComputeDEDXPerVolume ( const G4Material *  material, const G4ParticleDefinition *  particle, G4double  kineticEnergy, G4double  cutEnergy  )

virtual

Reimplemented from G4VEmModel.

Definition at line 436 of file G4IonParametrisedLossModel.cc.

References DeltaRayMeanEnergyTransferRate(), GetChargeSquareRatio(), and G4VEmModel::LowEnergyLimit().

Referenced by CorrectionsAlongStep(), and PrintDEDXTable().

                                          {

  // ############## dE/dx ##################################################
  // Function computes dE/dx values, where following rules are adopted:
  //   A. If the ion-material pair is covered by any native ion data
  //      parameterisation, then: 
  //      * This parameterization is used for energies below a given energy 
  //        limit,
  //      * whereas above the limit the Bethe-Bloch model is applied, in 
  //        combination with an effective charge estimate and high order 
  //        correction terms.
  //      A smoothing procedure is applied to dE/dx values computed with
  //      the second approach. The smoothing factor is based on the dE/dx
  //      values of both approaches at the transition energy (high order
  //      correction terms are included in the calculation of the transition
  //      factor).
  //   B. If the particle is a generic ion, the BraggIon and Bethe-Bloch
  //      models are used and a smoothing procedure is applied to values
  //      obtained with the second approach.
  //   C. If the ion-material is not covered by any ion data parameterization
  //      then:
  //      * The BraggIon model is used for energies below a given energy 
  //        limit,
  //      * whereas above the limit the Bethe-Bloch model is applied, in 
  //        combination with an effective charge estimate and high order 
  //        correction terms.
  //      Also in this case, a smoothing procedure is applied to dE/dx values 
  //      computed with the second model. 

  G4double dEdx = 0.0;
  UpdateDEDXCache(particle, material, cutEnergy);

  LossTableList::iterator iter = dedxCacheIter;

  if(iter != lossTableList.end()) {

     G4double transitionEnergy = dedxCacheTransitionEnergy;

     if(transitionEnergy > kineticEnergy) {

        dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);

        G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material, 
                                           particle, 
                                           kineticEnergy, 
                                           cutEnergy);
        dEdx -= dEdxDeltaRays;    
     }
     else {
        G4double massRatio = dedxCacheGenIonMassRatio;
                       
        G4double chargeSquare = 
                       GetChargeSquareRatio(particle, material, kineticEnergy);

        G4double scaledKineticEnergy = kineticEnergy * massRatio;
        G4double scaledTransitionEnergy = transitionEnergy * massRatio;

        G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();

        if(scaledTransitionEnergy >= lowEnergyLimit) {

           dEdx = betheBlochModel -> ComputeDEDXPerVolume(
                                      material, genericIon,
                          scaledKineticEnergy, cutEnergy);
        
           dEdx *= chargeSquare;

           dEdx += corrections -> ComputeIonCorrections(particle, 
                                                 material, kineticEnergy);

           G4double factor = 1.0 + dedxCacheTransitionFactor / 
                                                           kineticEnergy;

       dEdx *= factor;
    }
     }
  }
  else {
     G4double massRatio = 1.0;
     G4double chargeSquare = 1.0;

     if(particle != genericIon) {

        chargeSquare = GetChargeSquareRatio(particle, material, kineticEnergy);
        massRatio = genericIonPDGMass / particle -> GetPDGMass(); 
     }

     G4double scaledKineticEnergy = kineticEnergy * massRatio;

     G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
     if(scaledKineticEnergy < lowEnergyLimit) {
        dEdx = braggIonModel -> ComputeDEDXPerVolume(
                                      material, genericIon,
                          scaledKineticEnergy, cutEnergy);

        dEdx *= chargeSquare;
     }
     else {
        G4double dEdxLimitParam = braggIonModel -> ComputeDEDXPerVolume(
                                      material, genericIon,
                          lowEnergyLimit, cutEnergy);

        G4double dEdxLimitBetheBloch = betheBlochModel -> ComputeDEDXPerVolume(
                                      material, genericIon,
                          lowEnergyLimit, cutEnergy);

        if(particle != genericIon) {
           G4double chargeSquareLowEnergyLimit = 
               GetChargeSquareRatio(particle, material, 
                                    lowEnergyLimit / massRatio);

           dEdxLimitParam *= chargeSquareLowEnergyLimit;
           dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;

           dEdxLimitBetheBloch += 
                    corrections -> ComputeIonCorrections(particle, 
                                      material, lowEnergyLimit / massRatio);
        }

        G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
                               * lowEnergyLimit / scaledKineticEnergy);

        dEdx = betheBlochModel -> ComputeDEDXPerVolume(
                                      material, genericIon,
                          scaledKineticEnergy, cutEnergy);

        dEdx *= chargeSquare;

        if(particle != genericIon) {
           dEdx += corrections -> ComputeIonCorrections(particle, 
                                      material, kineticEnergy);
        }

        dEdx *= factor;
     }

  }

  if (dEdx < 0.0) dEdx = 0.0;

  return dEdx;
}

G4double G4IonParametrisedLossModel::ComputeLossForStep	(	const G4MaterialCutsCouple *	matCutsCouple,
		const G4ParticleDefinition *	particle,
		G4double	kineticEnergy,
		G4double	stepLength
	)

Definition at line 1167 of file G4IonParametrisedLossModel.cc.

References energy(), G4cout, G4endl, python.hepunit::mm, and G4VEmModel::Value().

Referenced by CorrectionsAlongStep().

                                          {

  G4double loss = 0.0;

  UpdateRangeCache(particle, matCutsCouple);

  G4PhysicsVector* energyRange = rangeCacheEnergyRange;
  G4PhysicsVector* rangeEnergy = rangeCacheRangeEnergy;

  if(energyRange != 0 && rangeEnergy != 0) {

     G4double lowerEnEdge = energyRange -> Energy( 0 );
     G4double lowerRangeEdge = rangeEnergy -> Energy( 0 );

     // Computing range for pre-step kinetic energy:
     G4double range = energyRange -> Value(kineticEnergy);

     // Energy below vector boundary:
     if(kineticEnergy < lowerEnEdge) {

        range =  energyRange -> Value(lowerEnEdge);
        range *= std::sqrt(kineticEnergy / lowerEnEdge);
     }

#ifdef PRINT_DEBUG
     G4cout << "G4IonParametrisedLossModel::ComputeLossForStep() range = " 
            << range / mm << " mm, step = " << stepLength / mm << " mm" 
            << G4endl;
#endif

     // Remaining range:
     G4double remRange = range - stepLength;

     // If range is smaller than step length, the loss is set to kinetic  
     // energy
     if(remRange < 0.0) loss = kineticEnergy;
     else if(remRange < lowerRangeEdge) {
 
        G4double ratio = remRange / lowerRangeEdge;
        loss = kineticEnergy - ratio * ratio * lowerEnEdge;
     }
     else {

        G4double energy = rangeEnergy -> Value(range - stepLength);
        loss = kineticEnergy - energy;      
     }
  }

  if(loss < 0.0) loss = 0.0;

  return loss;
}

void G4IonParametrisedLossModel::CorrectionsAlongStep ( const G4MaterialCutsCouple *  couple, const G4DynamicParticle *  dynamicParticle, G4double &  eloss, G4double &  , G4double  length  )

virtual

Reimplemented from G4VEmModel.

Definition at line 901 of file G4IonParametrisedLossModel.cc.

References ComputeDEDXPerVolume(), ComputeLossForStep(), DBL_MAX, energy(), G4cout, G4endl, G4VEmModel::GetModelOfFluctuations(), G4VEmModel::LowEnergyLimit(), python.hepunit::MeV, and right.

                                              {

  // ############## Corrections for along step energy loss calculation ######
  // The computed energy loss (due to electronic stopping) is overwritten 
  // by this function if an ion data parameterization is available for the 
  // current ion-material pair.
  // No action on the energy loss (due to electronic stopping) is performed 
  // if no parameterization is available. In this case the original 
  // generic ion tables (in combination with the effective charge) are used
  // in the along step DoIt function.
  //
  //
  // (Implementation partly adapted from G4BraggIonModel/G4BetheBlochModel)

  const G4ParticleDefinition* particle = dynamicParticle -> GetDefinition();
  const G4Material* material = couple -> GetMaterial();

  G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();

  if(kineticEnergy == eloss) { return; }

  G4double cutEnergy = DBL_MAX;
  size_t cutIndex = couple -> GetIndex();
  cutEnergy = cutEnergies[cutIndex];

  UpdateDEDXCache(particle, material, cutEnergy);

  LossTableList::iterator iter = dedxCacheIter;

  // If parameterization for ions is available the electronic energy loss
  // is overwritten
  if(iter != lossTableList.end()) {

     // The energy loss is calculated using the ComputeDEDXPerVolume function
     // and the step length (it is assumed that dE/dx does not change 
     // considerably along the step)
     eloss = 
        length * ComputeDEDXPerVolume(material, particle, 
                                      kineticEnergy, cutEnergy);

#ifdef PRINT_DEBUG
  G4cout.precision(6);      
  G4cout << "########################################################"
         << G4endl
         << "# G4IonParametrisedLossModel::CorrectionsAlongStep"
         << G4endl
         << "# cut(MeV) = " << cutEnergy/MeV 
         << G4endl;

  G4cout << "#" 
         << std::setw(13) << std::right << "E(MeV)"
         << std::setw(14) << "l(um)"
         << std::setw(14) << "l*dE/dx(MeV)"
         << std::setw(14) << "(l*dE/dx)/E"
         << G4endl
         << "# ------------------------------------------------------"
         << G4endl;

  G4cout << std::setw(14) << std::right << kineticEnergy / MeV
         << std::setw(14) << length / um
         << std::setw(14) << eloss / MeV
         << std::setw(14) << eloss / kineticEnergy * 100.0
         << G4endl;
#endif

     // If the energy loss exceeds a certain fraction of the kinetic energy
     // (the fraction is indicated by the parameter "energyLossLimit") then
     // the range tables are used to derive a more accurate value of the
     // energy loss
     if(eloss > energyLossLimit * kineticEnergy) {

        eloss = ComputeLossForStep(couple, particle, 
                                   kineticEnergy,length);

#ifdef PRINT_DEBUG
  G4cout << "# Correction applied:"
         << G4endl;

  G4cout << std::setw(14) << std::right << kineticEnergy / MeV
         << std::setw(14) << length / um
         << std::setw(14) << eloss / MeV
         << std::setw(14) << eloss / kineticEnergy * 100.0
         << G4endl;
#endif

     }
  }

  // For all corrections below a kinetic energy between the Pre- and 
  // Post-step energy values is used
  G4double energy = kineticEnergy - eloss * 0.5;
  if(energy < 0.0) energy = kineticEnergy * 0.5;

  G4double chargeSquareRatio = corrections ->
                                     EffectiveChargeSquareRatio(particle,
                                                                material,
                                                                energy);
  GetModelOfFluctuations() -> SetParticleAndCharge(particle, 
                                                   chargeSquareRatio);

  // A correction is applied considering the change of the effective charge 
  // along the step (the parameter "corrFactor" refers to the effective 
  // charge at the beginning of the step). Note: the correction is not
  // applied for energy loss values deriving directly from parameterized 
  // ion stopping power tables
  G4double transitionEnergy = dedxCacheTransitionEnergy;

  if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
     chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
                                                                material,
                                                                energy);

     G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
     eloss *= chargeSquareRatioCorr;
  }
  else if (iter == lossTableList.end()) {

     chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
                                                                material,
                                                                energy);

     G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
     eloss *= chargeSquareRatioCorr;
  }

  // Ion high order corrections are applied if the current model does not
  // overwrite the energy loss (i.e. when the effective charge approach is
  // used)
  if(iter == lossTableList.end()) {

     G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
     G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();

     // Corrections are only applied in the Bethe-Bloch energy region
     if(scaledKineticEnergy > lowEnergyLimit) 
       eloss += length * 
            corrections -> IonHighOrderCorrections(particle, couple, energy);
  }
}

G4double G4IonParametrisedLossModel::CrossSectionPerVolume ( const G4Material *  material, const G4ParticleDefinition *  particle, G4double  kineticEnergy, G4double  cutEnergy, G4double  maxEnergy  )

virtual

Reimplemented from G4VEmModel.

Definition at line 417 of file G4IonParametrisedLossModel.cc.

References ComputeCrossSectionPerAtom().

                                                       {

  G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
  G4double cross = ComputeCrossSectionPerAtom(particle, 
                                              kineticEnergy, 
                                              nbElecPerVolume, 0,
                                              cutEnergy,
                                              maxEnergy);

  return cross;
}

void G4IonParametrisedLossModel::DeactivateICRU73Scaling

(

)

Definition at line 1308 of file G4IonParametrisedLossModel.cc.

References AddDEDXTable(), and RemoveDEDXTable().

                                                         {

  RemoveDEDXTable("ICRU73");
  AddDEDXTable("ICRU73", new G4IonStoppingData("ion_stopping_data/icru73"));
}

G4double G4IonParametrisedLossModel::DeltaRayMeanEnergyTransferRate ( const G4Material *  , const G4ParticleDefinition *  , G4double  , G4double    )

inline

Referenced by ComputeDEDXPerVolume().

G4double G4IonParametrisedLossModel::GetChargeSquareRatio ( const G4ParticleDefinition *  particle, const G4Material *  material, G4double  kineticEnergy  )

virtual

Reimplemented from G4VEmModel.

Definition at line 210 of file G4IonParametrisedLossModel.cc.

Referenced by ComputeDEDXPerVolume().

                                                     {    // Kinetic energy

  G4double chargeSquareRatio = corrections ->
                                     EffectiveChargeSquareRatio(particle,
                                                                material,
                                                                kineticEnergy);
  corrFactor = chargeSquareRatio * 
                       corrections -> EffectiveChargeCorrection(particle,
                                                                material,
                                                                kineticEnergy);
  return corrFactor;
}

G4double G4IonParametrisedLossModel::GetParticleCharge ( const G4ParticleDefinition *  particle, const G4Material *  material, G4double  kineticEnergy  )

virtual

Reimplemented from G4VEmModel.

Definition at line 228 of file G4IonParametrisedLossModel.cc.

                                                     {   // Kinetic energy 

  return corrections -> GetParticleCharge(particle, material, kineticEnergy);
}

void G4IonParametrisedLossModel::Initialise ( const G4ParticleDefinition *  particle, const G4DataVector &  cuts  )

virtual

Implements G4VEmModel.

Definition at line 238 of file G4IonParametrisedLossModel.cc.

References G4cout, G4endl, G4VEmModel::GetName(), G4VEmModel::GetParticleChangeForLoss(), G4ProductionCutsTable::GetProductionCutsTable(), eplot::material, python.hepunit::second, and G4VEmModel::SetParticleChange().

                                                                 {

  // Cached parameters are reset
  cacheParticle = 0;
  cacheMass = 0;
  cacheElecMassRatio = 0;
  cacheChargeSquare = 0;
  
  // Cached parameters are reset
  rangeCacheParticle = 0; 
  rangeCacheMatCutsCouple = 0; 
  rangeCacheEnergyRange = 0; 
  rangeCacheRangeEnergy = 0;

  // Cached parameters are reset
  dedxCacheParticle = 0;
  dedxCacheMaterial = 0;
  dedxCacheEnergyCut = 0;
  dedxCacheIter = lossTableList.end();
  dedxCacheTransitionEnergy = 0.0;  
  dedxCacheTransitionFactor = 0.0;
  dedxCacheGenIonMassRatio = 0.0;

  // The cache of loss tables is cleared
  LossTableList::iterator iterTables = lossTableList.begin();
  LossTableList::iterator iterTables_end = lossTableList.end();

  for(;iterTables != iterTables_end; iterTables++) 
                                       (*iterTables) -> ClearCache();

  // Range vs energy and energy vs range vectors from previous runs are
  // cleared
  RangeEnergyTable::iterator iterRange = r.begin();
  RangeEnergyTable::iterator iterRange_end = r.end();

  for(;iterRange != iterRange_end; iterRange++) delete iterRange -> second;
  r.clear();

  EnergyRangeTable::iterator iterEnergy = E.begin();
  EnergyRangeTable::iterator iterEnergy_end = E.end();

  for(;iterEnergy != iterEnergy_end; iterEnergy++) delete iterEnergy -> second;
  E.clear();

  // The cut energies are (re)loaded
  size_t size = cuts.size();
  cutEnergies.clear();
  for(size_t i = 0; i < size; i++) cutEnergies.push_back(cuts[i]);

  // All dE/dx vectors are built
  const G4ProductionCutsTable* coupleTable=
                     G4ProductionCutsTable::GetProductionCutsTable();
  size_t nmbCouples = coupleTable -> GetTableSize();

#ifdef PRINT_TABLE_BUILT
    G4cout << "G4IonParametrisedLossModel::Initialise():"
           << " Building dE/dx vectors:"
           << G4endl;         
#endif

  for (size_t i = 0; i < nmbCouples; i++) {

    const G4MaterialCutsCouple* couple = 
                                     coupleTable -> GetMaterialCutsCouple(i);

    const G4Material* material = couple -> GetMaterial();
    //    G4ProductionCuts* productionCuts = couple -> GetProductionCuts();

    for(G4int atomicNumberIon = 3; atomicNumberIon < 102; atomicNumberIon++) {

       LossTableList::iterator iter = lossTableList.begin();
       LossTableList::iterator iter_end = lossTableList.end();

       for(;iter != iter_end; iter++) { 

          if(*iter == 0) {
              G4cout << "G4IonParametrisedLossModel::Initialise():"
                     << " Skipping illegal table."
                     << G4endl;         
          }

          G4bool isApplicable = 
                    (*iter) -> BuildDEDXTable(atomicNumberIon, material);
          if(isApplicable) {

#ifdef PRINT_TABLE_BUILT
             G4cout << "  Atomic Number Ion = " << atomicNumberIon
                    << ", Material = " << material -> GetName()
                    << ", Table = " << (*iter) -> GetName()
                    << G4endl;      
#endif
             break; 
      }
       }
    }
  }

  // The particle change object 
  if(! particleChangeLoss) {
    particleChangeLoss = GetParticleChangeForLoss();
    braggIonModel -> SetParticleChange(particleChangeLoss, 0);
    betheBlochModel -> SetParticleChange(particleChangeLoss, 0);
  }
 
  // The G4BraggIonModel and G4BetheBlochModel instances are initialised with
  // the same settings as the current model:
  braggIonModel -> Initialise(particle, cuts);
  betheBlochModel -> Initialise(particle, cuts);
}

LossTableList::iterator G4IonParametrisedLossModel::IsApplicable ( const G4ParticleDefinition *  , const G4Material *    )

inline

Referenced by PrintDEDXTableHandlers().

G4double G4IonParametrisedLossModel::MaxSecondaryEnergy ( const G4ParticleDefinition *  particle, G4double  kineticEnergy  )

protectedvirtual

Reimplemented from G4VEmModel.

Definition at line 182 of file G4IonParametrisedLossModel.cc.

References python.hepunit::electron_mass_c2.

Referenced by ComputeCrossSectionPerAtom().

                                                     {

  // ############## Maximum energy of secondaries ##########################
  // Function computes maximum energy of secondary electrons which are 
  // released by an ion
  //
  // See Geant4 physics reference manual (version 9.1), section 9.1.1
  // 
  // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
  //       C.Caso et al. (Part. Data Group), Europ. Phys. Jour. C 3 1 (1998).
  //       B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
  //
  // (Implementation adapted from G4BraggIonModel)

  if(particle != cacheParticle) UpdateCache(particle);

  G4double tau  = kineticEnergy/cacheMass;
  G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
                  (1. + 2.0 * (tau + 1.) * cacheElecMassRatio + 
                  cacheElecMassRatio * cacheElecMassRatio);

  return tmax;
}

G4double G4IonParametrisedLossModel::MinEnergyCut ( const G4ParticleDefinition *  , const G4MaterialCutsCouple *  couple  )

virtual

Reimplemented from G4VEmModel.

Definition at line 172 of file G4IonParametrisedLossModel.cc.

                                                                           {

  return couple -> GetMaterial() -> GetIonisation() -> 
                                                  GetMeanExcitationEnergy();
}

void G4IonParametrisedLossModel::PrintDEDXTable	(	const G4ParticleDefinition *	particle,
		const G4Material *	material,
		G4double	lowerBoundary,
		G4double	upperBoundary,
		G4int	numBins,
		G4bool	logScaleEnergy
	)

Definition at line 585 of file G4IonParametrisedLossModel.cc.

References python.hepunit::cm, python.hepunit::cm2, python.hepunit::cm3, ComputeDEDXPerVolume(), DBL_MAX, energy(), g(), G4cout, G4endl, G4VEmModel::GetName(), python.hepunit::MeV, and right.

Referenced by PrintDEDXTableHandlers().

                                          {     // Logarithmic scaling of energy

  G4double atomicMassNumber = particle -> GetAtomicMass();
  G4double materialDensity = material -> GetDensity();

  G4cout << "# dE/dx table for " << particle -> GetParticleName() 
         << " in material " << material -> GetName()
         << " of density " << materialDensity / g * cm3
         << " g/cm3"
         << G4endl
         << "# Projectile mass number A1 = " << atomicMassNumber
         << G4endl
         << "# ------------------------------------------------------"
         << G4endl;
  G4cout << "#"
         << std::setw(13) << std::right << "E"
         << std::setw(14) << "E/A1"
         << std::setw(14) << "dE/dx"
         << std::setw(14) << "1/rho*dE/dx"
         << G4endl;
  G4cout << "#"
         << std::setw(13) << std::right << "(MeV)"
         << std::setw(14) << "(MeV)"
         << std::setw(14) << "(MeV/cm)"
         << std::setw(14) << "(MeV*cm2/mg)"
         << G4endl
         << "# ------------------------------------------------------"
         << G4endl;

  G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
  G4double energyUpperBoundary = upperBoundary * atomicMassNumber; 

  if(logScaleEnergy) {

     energyLowerBoundary = std::log(energyLowerBoundary);
     energyUpperBoundary = std::log(energyUpperBoundary); 
  }

  G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) / 
                                                           G4double(nmbBins);

  for(int i = 0; i < numBins + 1; i++) {

      G4double energy = energyLowerBoundary + i * deltaEnergy;
      if(logScaleEnergy) energy = std::exp(energy);

      G4double dedx = ComputeDEDXPerVolume(material, particle, energy, DBL_MAX);
      G4cout.precision(6);
      G4cout << std::setw(14) << std::right << energy / MeV
             << std::setw(14) << energy / atomicMassNumber / MeV
             << std::setw(14) << dedx / MeV * cm
             << std::setw(14) << dedx / materialDensity / (MeV*cm2/(0.001*g)) 
             << G4endl;
  }
}

void G4IonParametrisedLossModel::PrintDEDXTableHandlers	(	const G4ParticleDefinition *	particle,
		const G4Material *	material,
		G4double	lowerBoundary,
		G4double	upperBoundary,
		G4int	numBins,
		G4bool	logScaleEnergy
	)

Definition at line 649 of file G4IonParametrisedLossModel.cc.

References IsApplicable(), and PrintDEDXTable().

                                          {     // Logarithmic scaling of energy

  LossTableList::iterator iter = lossTableList.begin();
  LossTableList::iterator iter_end = lossTableList.end();

  for(;iter != iter_end; iter++) { 
      G4bool isApplicable =  (*iter) -> IsApplicable(particle, material);
      if(isApplicable) {  
    (*iter) -> PrintDEDXTable(particle, material,
                                  lowerBoundary, upperBoundary, 
                                  numBins,logScaleEnergy); 
        break;
      } 
  }
}

G4bool G4IonParametrisedLossModel::RemoveDEDXTable

(

const G4String &

name

)

Definition at line 1269 of file G4IonParametrisedLossModel.cc.

References G4VEmModel::GetName(), and python.hepunit::second.

Referenced by DeactivateICRU73Scaling().

                                      {

  LossTableList::iterator iter = lossTableList.begin();
  LossTableList::iterator iter_end = lossTableList.end();
  
  for(;iter != iter_end; iter++) {
     G4String tableName = (*iter) -> GetName();

     if(tableName == nam) { 
        delete (*iter);

        // Remove from table list
        lossTableList.erase(iter);

        // Range vs energy and energy vs range vectors are cleared
        RangeEnergyTable::iterator iterRange = r.begin();
        RangeEnergyTable::iterator iterRange_end = r.end();

        for(;iterRange != iterRange_end; iterRange++) 
                                            delete iterRange -> second;
        r.clear();

        EnergyRangeTable::iterator iterEnergy = E.begin();
        EnergyRangeTable::iterator iterEnergy_end = E.end();

        for(;iterEnergy != iterEnergy_end; iterEnergy++) 
                                            delete iterEnergy -> second;
        E.clear();

        return true;
     }
  }

  return false;
}

void G4IonParametrisedLossModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  secondaries, const G4MaterialCutsCouple *  , const G4DynamicParticle *  particle, G4double  cutKinEnergySec, G4double  userMaxKinEnergySec  )

virtual

Implements G4VEmModel.

Definition at line 673 of file G4IonParametrisedLossModel.cc.

References G4Electron::Definition(), python.hepunit::electron_mass_c2, energy(), G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetMomentumDirection(), G4VEmModel::MaxSecondaryKinEnergy(), G4INCL::Math::min(), CLHEP::Hep3Vector::rotateUz(), python.hepunit::twopi, and CLHEP::Hep3Vector::unit().

                                               {


  // ############## Sampling of secondaries #################################
  // The probability density function (pdf) of the kinetic energy T of a 
  // secondary electron may be written as:
  //    pdf(T) = f(T) * g(T)
  // where 
  //    f(T) = (Tmax - Tcut) / (Tmax * Tcut) * (1 / T^2)
  //    g(T) = 1 - beta^2 * T / Tmax
  // where Tmax is the maximum kinetic energy of the secondary, Tcut
  // is the lower energy cut and beta is the kinetic energy of the 
  // projectile.
  //
  // Sampling of the kinetic energy of a secondary electron:
  //  1) T0 is sampled from f(T) using the cumulated distribution function
  //     F(T) = int_Tcut^T f(T')dT'
  //  2) T is accepted or rejected by evaluating the rejection function g(T)
  //     at the sampled energy T0 against a randomly sampled value
  //
  //
  // See Geant4 physics reference manual (version 9.1), section 9.1.4
  //
  //
  // Reference pdf: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
  //
  // (Implementation adapted from G4BraggIonModel)

  G4double rossiMaxKinEnergySec = MaxSecondaryKinEnergy(particle);
  G4double maxKinEnergySec = 
                        std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);

  if(cutKinEnergySec >= maxKinEnergySec) return;

  G4double kineticEnergy = particle -> GetKineticEnergy();
  G4ThreeVector direction = particle ->GetMomentumDirection();

  G4double energy  = kineticEnergy + cacheMass;
  G4double betaSquared  = kineticEnergy * 
                                    (energy + cacheMass) / (energy * energy);

  G4double kinEnergySec;
  G4double grej;

  do {

    // Sampling kinetic energy from f(T) (using F(T)):
    G4double xi = G4UniformRand();
    kinEnergySec = cutKinEnergySec * maxKinEnergySec / 
                        (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);

    // Deriving the value of the rejection function at the obtained kinetic 
    // energy:
    grej = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;

    if(grej > 1.0) {
        G4cout << "G4IonParametrisedLossModel::SampleSecondary Warning: "
               << "Majorant 1.0 < "
               << grej << " for e= " << kinEnergySec
               << G4endl;
    }

  } while( G4UniformRand() >= grej );

  G4double momentumSec =
           std::sqrt(kinEnergySec * (kinEnergySec + 2.0 * electron_mass_c2));

  G4double totMomentum = energy*std::sqrt(betaSquared);
  G4double cost = kinEnergySec * (energy + electron_mass_c2) /
                                   (momentumSec * totMomentum);
  if(cost > 1.0) cost = 1.0;
  G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));

  G4double phi = twopi * G4UniformRand() ;

  G4ThreeVector directionSec(sint*std::cos(phi),sint*std::sin(phi), cost) ;
  directionSec.rotateUz(direction);

  // create G4DynamicParticle object for delta ray
  G4DynamicParticle* delta = new G4DynamicParticle(G4Electron::Definition(),
                                                   directionSec,
                           kinEnergySec);

  secondaries -> push_back(delta);

  // Change kinematics of primary particle
  kineticEnergy       -= kinEnergySec;
  G4ThreeVector finalP = direction*totMomentum - directionSec*momentumSec;
  finalP               = finalP.unit();

  particleChangeLoss -> SetProposedKineticEnergy(kineticEnergy);
  particleChangeLoss -> SetProposedMomentumDirection(finalP);
}

void G4IonParametrisedLossModel::SetEnergyLossLimit ( G4double  ionEnergyLossLimit )

inline

---

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

• G4IonParametrisedLossModel.hh
• G4IonParametrisedLossModel.cc