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)
Protected Member Functions
virtual G4double	MaxSecondaryEnergy (const G4ParticleDefinition *, G4double)

---

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(), and G4VEmModel::SetHighEnergyLimit().

  : 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.

                                                        {

  // Range vs energy table objects are deleted and the container is cleared
  RangeEnergyTable::iterator iterRange = r.begin();
  RangeEnergyTable::iterator iterRange_end = r.end();

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

  // Energy vs range table objects are deleted and the container is cleared
  EnergyRangeTable::iterator iterEnergy = E.begin();
  EnergyRangeTable::iterator iterEnergy_end = E.end();

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

  // 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();

  // The Bragg ion and Bethe Bloch objects are deleted
  delete betheBlochModel;
  delete braggIonModel;
}

---

Member Function Documentation

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

Definition at line 1246 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 *	,
		G4double	,
		G4double	,
		G4double	,
		G4double	,
		G4double
	)			[virtual]

Reimplemented from G4VEmModel.

Definition at line 370 of file G4IonParametrisedLossModel.cc.

References G4cout, G4endl, and MaxSecondaryEnergy().

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 *	,
		const G4ParticleDefinition *	,
		G4double	,
		G4double
	)			[virtual]

Reimplemented from G4VEmModel.

Definition at line 454 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 *	,
		const G4ParticleDefinition *	,
		G4double	,
		G4double
	)

Definition at line 1186 of file G4IonParametrisedLossModel.cc.

References G4cout, and G4endl.

Referenced by CorrectionsAlongStep().

                                          {

  G4double loss = 0.0;

  UpdateRangeCache(particle, matCutsCouple);

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

  if(energyRange != 0 && rangeEnergy != 0) {
     G4bool b;

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

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

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

        range =  energyRange -> GetValue(lowerEnEdge, b);
        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 -> GetValue(range - stepLength, b);
        loss = kineticEnergy - energy;      
     }
  }

  if(loss < 0.0) loss = 0.0;

  return loss;
}

void G4IonParametrisedLossModel::CorrectionsAlongStep	(	const G4MaterialCutsCouple *	,
		const G4DynamicParticle *	,
		G4double &	,
		G4double &	,
		G4double
	)			[virtual]

Reimplemented from G4VEmModel.

Definition at line 919 of file G4IonParametrisedLossModel.cc.

References ComputeDEDXPerVolume(), ComputeLossForStep(), DBL_MAX, G4cout, G4endl, G4VEmModel::GetModelOfFluctuations(), and G4VEmModel::LowEnergyLimit().

                                              {

  // ############## 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 *	,
		const G4ParticleDefinition *	,
		G4double	,
		G4double	,
		G4double
	)			[virtual]

Reimplemented from G4VEmModel.

Definition at line 435 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 1328 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]

Definition at line 58 of file G4IonParametrisedLossModel.icc.

References GetChargeSquareRatio(), G4Material::GetTotNbOfElectPerVolume(), and MaxSecondaryEnergy().

Referenced by ComputeDEDXPerVolume().

                                                          {

  // ############## Mean energy transferred to delta-rays ###################
  // Computes the mean energy transfered to delta-rays per unit length,
  // considering only delta-rays with energies above the energy threshold 
  // (energy cut)
  //
  // The mean energy transfer rate is derived by using the differential
  // cross section given in the references below.
  //
  // 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)


  //   *** Variables:
  //   kineticEnergy = kinetic energy of projectile
  //   totEnergy     = total energy of projectile, i.e. kinetic energy
  //                   plus rest energy (Mc^2)
  //   betaSquared   = beta of projectile squared, calculated as
  //                      beta^2 = 1 - 1 / (E/Mc^2)^2
  //                             = T * ( E + Mc^2 ) / E^2
  //                   where T = kineticEnergy, E = totEnergy
  //   cutEnergy     = energy threshold for secondary particle production
  //                   i.e. energy cut, below which energy transfered to 
  //                   electrons is treated as continuous loss of projectile
  //   maxKinEnergy  = maximum energy transferable to secondary electrons
  //   meanRate      = mean kinetic energy of delta ray (per unit length) 
  //                   (above cutEnergy)  

  G4double meanRate = 0.0;

  G4double maxKinEnergy = MaxSecondaryEnergy(particle, kineticEnergy);

  if (cutEnergy < maxKinEnergy) {

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

    G4double cutMaxEnergyRatio = cutEnergy / maxKinEnergy;

    meanRate = 
        (- std::log(cutMaxEnergyRatio) - (1.0 - cutMaxEnergyRatio) * betaSquared) * 
        CLHEP::twopi_mc2_rcl2 * 
        (material->GetTotNbOfElectPerVolume()) / betaSquared;

    meanRate *= GetChargeSquareRatio(particle, material, kineticEnergy);
  }
  
  return meanRate;
}

G4double G4IonParametrisedLossModel::GetChargeSquareRatio	(	const G4ParticleDefinition *	,
		const G4Material *	,
		G4double
	)			[virtual]

Reimplemented from G4VEmModel.

Definition at line 228 of file G4IonParametrisedLossModel.cc.

Referenced by ComputeDEDXPerVolume(), and DeltaRayMeanEnergyTransferRate().

                                                     {    // Kinetic energy

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

G4double G4IonParametrisedLossModel::GetParticleCharge	(	const G4ParticleDefinition *	,
		const G4Material *	,
		G4double
	)			[virtual]

Reimplemented from G4VEmModel.

Definition at line 246 of file G4IonParametrisedLossModel.cc.

                                                     {   // Kinetic energy 

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

void G4IonParametrisedLossModel::Initialise	(	const G4ParticleDefinition *	,
		const G4DataVector &
	)			[virtual]

Implements G4VEmModel.

Definition at line 256 of file G4IonParametrisedLossModel.cc.

References G4cout, G4endl, G4VEmModel::GetName(), G4VEmModel::GetParticleChangeForLoss(), G4ProductionCutsTable::GetProductionCutsTable(), 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]

Definition at line 130 of file G4IonParametrisedLossModel.icc.

Referenced by PrintDEDXTableHandlers().

                                                {          // Target material

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

  for(;iterTables != iterTables_end; iterTables++) {
      G4bool isApplicable = (*iterTables) -> 
                       IsApplicable(particle, material);
      if(isApplicable) {
         iter = iterTables;
         break;
      }
  }

  return iter;
}

G4double G4IonParametrisedLossModel::MaxSecondaryEnergy	(	const G4ParticleDefinition *	,
		G4double
	)			[protected, virtual]

Reimplemented from G4VEmModel.

Definition at line 200 of file G4IonParametrisedLossModel.cc.

Referenced by ComputeCrossSectionPerAtom(), and DeltaRayMeanEnergyTransferRate().

                                                     {

  // ############## 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 *
	)			[virtual]

Definition at line 190 of file G4IonParametrisedLossModel.cc.

                                                                           {

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

void G4IonParametrisedLossModel::PrintDEDXTable	(	const G4ParticleDefinition *	,
		const G4Material *	,
		G4double	,
		G4double	,
		G4int	,
		G4bool
	)

Definition at line 603 of file G4IonParametrisedLossModel.cc.

References ComputeDEDXPerVolume(), DBL_MAX, G4cout, G4endl, and G4VEmModel::GetName().

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 *	,
		const G4Material *	,
		G4double	,
		G4double	,
		G4int	,
		G4bool
	)

Definition at line 667 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 1289 of file G4IonParametrisedLossModel.cc.

References G4VEmModel::GetName().

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 * > *	,
		const G4MaterialCutsCouple *	,
		const G4DynamicParticle *	,
		G4double	,
		G4double
	)			[virtual]

Implements G4VEmModel.

Definition at line 691 of file G4IonParametrisedLossModel.cc.

References G4Electron::Definition(), G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetMomentumDirection(), and G4VEmModel::MaxSecondaryKinEnergy().

                                                           {


  // ############## 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]

Definition at line 152 of file G4IonParametrisedLossModel.icc.

                                                                         {

  if(ionEnergyLossLimit > 0 && ionEnergyLossLimit <=1) {

     energyLossLimit = ionEnergyLossLimit;
  }
}

---

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

• G4IonParametrisedLossModel.hh
• G4IonParametrisedLossModel.icc
• G4IonParametrisedLossModel.cc

---