G4hImpactIonisation Class Reference

#include <G4hImpactIonisation.hh>

Inheritance diagram for G4hImpactIonisation:

Public Member Functions
	G4hImpactIonisation (const G4String &processName="hImpactIoni")
	~G4hImpactIonisation ()
G4bool	IsApplicable (const G4ParticleDefinition &)
void	BuildPhysicsTable (const G4ParticleDefinition &aParticleType)
G4double	GetMeanFreePath (const G4Track &track, G4double previousStepSize, enum G4ForceCondition *condition)
void	PrintInfoDefinition () const
void	SetHighEnergyForProtonParametrisation (G4double energy)
void	SetLowEnergyForProtonParametrisation (G4double energy)
void	SetHighEnergyForAntiProtonParametrisation (G4double energy)
void	SetLowEnergyForAntiProtonParametrisation (G4double energy)
G4double	GetContinuousStepLimit (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety)
void	SetElectronicStoppingPowerModel (const G4ParticleDefinition *aParticle, const G4String &dedxTable)
void	SetNuclearStoppingPowerModel (const G4String &dedxTable)
void	SetNuclearStoppingOn ()
void	SetNuclearStoppingOff ()
void	SetBarkasOn ()
void	SetBarkasOff ()
void	SetPixe (const G4bool)
G4VParticleChange *	AlongStepDoIt (const G4Track &trackData, const G4Step &stepData)
G4VParticleChange *	PostStepDoIt (const G4Track &track, const G4Step &Step)
G4double	ComputeDEDX (const G4ParticleDefinition *aParticle, const G4MaterialCutsCouple *couple, G4double kineticEnergy)
void	SetCutForSecondaryPhotons (G4double cut)
void	SetCutForAugerElectrons (G4double cut)
void	ActivateAugerElectronProduction (G4bool val)
void	SetPixeCrossSectionK (const G4String &name)
void	SetPixeCrossSectionL (const G4String &name)
void	SetPixeCrossSectionM (const G4String &name)
void	SetPixeProjectileMinEnergy (G4double energy)
void	SetPixeProjectileMaxEnergy (G4double energy)
Public Member Functions inherited from G4hRDEnergyLoss
	G4hRDEnergyLoss (const G4String &)
	~G4hRDEnergyLoss ()
Public Member Functions inherited from G4VContinuousDiscreteProcess
	G4VContinuousDiscreteProcess (const G4String &, G4ProcessType aType=fNotDefined)
	G4VContinuousDiscreteProcess (G4VContinuousDiscreteProcess &)
virtual	~G4VContinuousDiscreteProcess ()
G4VContinuousDiscreteProcess &	operator= (const G4VContinuousDiscreteProcess &)=delete
virtual G4double	PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
virtual G4double	AlongStepGetPhysicalInteractionLength (const G4Track &, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety, G4GPILSelection *selection)
virtual G4double	AtRestGetPhysicalInteractionLength (const G4Track &, G4ForceCondition *)
virtual G4VParticleChange *	AtRestDoIt (const G4Track &, const G4Step &)
Public Member Functions inherited from G4VProcess
	G4VProcess (const G4String &aName="NoName", G4ProcessType aType=fNotDefined)
	G4VProcess (const G4VProcess &right)
virtual	~G4VProcess ()
G4VProcess &	operator= (const G4VProcess &)=delete
G4bool	operator== (const G4VProcess &right) const
G4bool	operator!= (const G4VProcess &right) const
G4double	GetCurrentInteractionLength () const
void	SetPILfactor (G4double value)
G4double	GetPILfactor () const
G4double	AlongStepGPIL (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)
G4double	AtRestGPIL (const G4Track &track, G4ForceCondition *condition)
G4double	PostStepGPIL (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
virtual void	PreparePhysicsTable (const G4ParticleDefinition &)
virtual G4bool	StorePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
virtual G4bool	RetrievePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
const G4String &	GetPhysicsTableFileName (const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
const G4String &	GetProcessName () const
G4ProcessType	GetProcessType () const
void	SetProcessType (G4ProcessType)
G4int	GetProcessSubType () const
void	SetProcessSubType (G4int)
virtual const G4VProcess *	GetCreatorProcess () const
virtual void	StartTracking (G4Track *)
virtual void	EndTracking ()
virtual void	SetProcessManager (const G4ProcessManager *)
virtual const G4ProcessManager *	GetProcessManager ()
virtual void	ResetNumberOfInteractionLengthLeft ()
G4double	GetNumberOfInteractionLengthLeft () const
G4double	GetTotalNumberOfInteractionLengthTraversed () const
G4bool	isAtRestDoItIsEnabled () const
G4bool	isAlongStepDoItIsEnabled () const
G4bool	isPostStepDoItIsEnabled () const
virtual void	DumpInfo () const
virtual void	ProcessDescription (std::ostream &outfile) const
void	SetVerboseLevel (G4int value)
G4int	GetVerboseLevel () const
virtual void	SetMasterProcess (G4VProcess *masterP)
const G4VProcess *	GetMasterProcess () const
virtual void	BuildWorkerPhysicsTable (const G4ParticleDefinition &part)
virtual void	PrepareWorkerPhysicsTable (const G4ParticleDefinition &)

Additional Inherited Members
Static Public Member Functions inherited from G4hRDEnergyLoss
static G4int	GetNumberOfProcesses ()
static void	SetNumberOfProcesses (G4int number)
static void	PlusNumberOfProcesses ()
static void	MinusNumberOfProcesses ()
static void	SetdRoverRange (G4double value)
static void	SetRndmStep (G4bool value)
static void	SetEnlossFluc (G4bool value)
static void	SetStepFunction (G4double c1, G4double c2)
Static Public Member Functions inherited from G4VProcess
static const G4String &	GetProcessTypeName (G4ProcessType)
Protected Member Functions inherited from G4hRDEnergyLoss
G4bool	CutsWhereModified ()
Protected Member Functions inherited from G4VContinuousDiscreteProcess
void	SetGPILSelection (G4GPILSelection selection)
G4GPILSelection	GetGPILSelection () const
Protected Member Functions inherited from G4VProcess
void	SubtractNumberOfInteractionLengthLeft (G4double prevStepSize)
void	ClearNumberOfInteractionLengthLeft ()
Static Protected Member Functions inherited from G4hRDEnergyLoss
static void	BuildDEDXTable (const G4ParticleDefinition &aParticleType)
Protected Attributes inherited from G4hRDEnergyLoss
const G4double	MaxExcitationNumber
const G4double	probLimFluct
const long	nmaxDirectFluct
const long	nmaxCont1
const long	nmaxCont2
G4PhysicsTable *	theLossTable
G4double	linLossLimit
G4double	MinKineticEnergy
Protected Attributes inherited from G4VProcess
const G4ProcessManager *	aProcessManager = nullptr
G4VParticleChange *	pParticleChange = nullptr
G4ParticleChange	aParticleChange
G4double	theNumberOfInteractionLengthLeft = -1.0
G4double	currentInteractionLength = -1.0
G4double	theInitialNumberOfInteractionLength = -1.0
G4String	theProcessName
G4String	thePhysicsTableFileName
G4ProcessType	theProcessType = fNotDefined
G4int	theProcessSubType = -1
G4double	thePILfactor = 1.0
G4int	verboseLevel = 0
G4bool	enableAtRestDoIt = true
G4bool	enableAlongStepDoIt = true
G4bool	enablePostStepDoIt = true
Static Protected Attributes inherited from G4hRDEnergyLoss
static G4ThreadLocal G4PhysicsTable *	theDEDXpTable = 0
static G4ThreadLocal G4PhysicsTable *	theDEDXpbarTable = 0
static G4ThreadLocal G4PhysicsTable *	theRangepTable = 0
static G4ThreadLocal G4PhysicsTable *	theRangepbarTable = 0
static G4ThreadLocal G4PhysicsTable *	theInverseRangepTable = 0
static G4ThreadLocal G4PhysicsTable *	theInverseRangepbarTable = 0
static G4ThreadLocal G4PhysicsTable *	theLabTimepTable = 0
static G4ThreadLocal G4PhysicsTable *	theLabTimepbarTable = 0
static G4ThreadLocal G4PhysicsTable *	theProperTimepTable = 0
static G4ThreadLocal G4PhysicsTable *	theProperTimepbarTable = 0
static G4ThreadLocal G4PhysicsTable **	RecorderOfpProcess = 0
static G4ThreadLocal G4PhysicsTable **	RecorderOfpbarProcess = 0
static G4ThreadLocal G4int	CounterOfpProcess = 0
static G4ThreadLocal G4int	CounterOfpbarProcess = 0
static G4ThreadLocal G4double	ParticleMass
static G4ThreadLocal G4double	ptableElectronCutInRange = 0.0
static G4ThreadLocal G4double	pbartableElectronCutInRange = 0.0
static G4ThreadLocal G4double	Charge
static G4ThreadLocal G4double	LowestKineticEnergy = 1e-05
static G4ThreadLocal G4double	HighestKineticEnergy = 1.e5
static G4ThreadLocal G4int	TotBin = 360
static G4ThreadLocal G4double	RTable
static G4ThreadLocal G4double	LOGRTable
static G4ThreadLocal G4double	dRoverRange = 0.20
static G4ThreadLocal G4double	finalRange = 0.2
static G4ThreadLocal G4double	c1lim = 0.20
static G4ThreadLocal G4double	c2lim = 0.32
static G4ThreadLocal G4double	c3lim = -0.032
static G4ThreadLocal G4bool	rndmStepFlag = false
static G4ThreadLocal G4bool	EnlossFlucFlag = true

Detailed Description

Definition at line 74 of file G4hImpactIonisation.hh.

Constructor & Destructor Documentation

G4hImpactIonisation()

G4hImpactIonisation::G4hImpactIonisation

(

const G4String &

processName = "hImpactIoni"

)

Definition at line 82 of file G4hImpactIonisation.cc.

  : G4hRDEnergyLoss(processName),
    betheBlochModel(0),
    protonModel(0),
    antiprotonModel(0),
    theIonEffChargeModel(0),
    theNuclearStoppingModel(0),
    theIonChuFluctuationModel(0),
    theIonYangFluctuationModel(0),
    protonTable("ICRU_R49p"),
    antiprotonTable("ICRU_R49p"),
    theNuclearTable("ICRU_R49"),
    nStopping(true),
    theBarkas(true),
    theMeanFreePathTable(0),
    paramStepLimit (0.005),
    pixeCrossSectionHandler(0)
{ 
  InitializeMe();
}

~G4hImpactIonisation()

G4hImpactIonisation::~G4hImpactIonisation

(

)

Definition at line 139 of file G4hImpactIonisation.cc.

{
  if (theMeanFreePathTable) 
    {
      theMeanFreePathTable->clearAndDestroy();
      delete theMeanFreePathTable;
    }
 
  if (betheBlochModel) delete betheBlochModel;
  if (protonModel) delete protonModel;
  if (antiprotonModel) delete antiprotonModel;
  if (theNuclearStoppingModel) delete theNuclearStoppingModel;
  if (theIonEffChargeModel) delete theIonEffChargeModel;
  if (theIonChuFluctuationModel) delete theIonChuFluctuationModel;
  if (theIonYangFluctuationModel) delete theIonYangFluctuationModel;
 
  delete pixeCrossSectionHandler;
 
  // ---- MGP ---- The following is to be checked
  //  if (shellVacancy) delete shellVacancy;
 
  cutForDelta.clear();
}

Member Function Documentation

ActivateAugerElectronProduction()

void G4hImpactIonisation::ActivateAugerElectronProduction

(

G4bool

val

)

Definition at line 1714 of file G4hImpactIonisation.cc.

{
  atomicDeexcitation.ActivateAugerElectronProduction(val);
}

AlongStepDoIt()

G4VParticleChange * G4hImpactIonisation::AlongStepDoIt ( const G4Track & trackData, const G4Step & stepData )

virtual

Reimplemented from G4VContinuousDiscreteProcess.

Definition at line 725 of file G4hImpactIonisation.cc.

{
  // compute the energy loss after a step
  G4Proton* proton = G4Proton::Proton();
  G4AntiProton* antiproton = G4AntiProton::AntiProton();
  G4double finalT = 0.;
 
  aParticleChange.Initialize(track) ;
 
  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
  const G4Material* material = couple->GetMaterial();
 
  // get the actual (true) Step length from step
  const G4double stepLength = step.GetStepLength() ;
 
  const G4DynamicParticle* particle = track.GetDynamicParticle() ;
 
  G4double kineticEnergy = particle->GetKineticEnergy() ;
  G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
  G4double tScaled = kineticEnergy * massRatio ;
  G4double eLoss = 0.0 ;
  G4double nLoss = 0.0 ;
 
 
  // very small particle energy
  if (kineticEnergy < MinKineticEnergy) 
    {
      eLoss = kineticEnergy ;
      // particle energy outside tabulated energy range
    } 
  
  else if( kineticEnergy > HighestKineticEnergy) 
    {
      eLoss = stepLength * fdEdx ;
      // big step
    } 
  else if (stepLength >= fRangeNow ) 
    {
      eLoss = kineticEnergy ;
      
      // tabulated range
    } 
  else 
    {
      // step longer than linear step limit
      if(stepLength > linLossLimit * fRangeNow) 
    {
      G4double rScaled = fRangeNow * massRatio * chargeSquare ;
      G4double sScaled = stepLength * massRatio * chargeSquare ;
      
      if(charge > 0.0) 
        {
          eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled, couple) -
        G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled-sScaled,couple) ;
        
        } 
      else 
        {
          // Antiproton
          eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled,couple) -
        G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled-sScaled,couple) ;
        }
      eLoss /= massRatio ;
      
      // Barkas correction at big step      
      eLoss += fBarkas * stepLength;
      
      // step shorter than linear step limit
    } 
      else 
    {
      eLoss = stepLength *fdEdx  ;
    }
      if (nStopping && tScaled < protonHighEnergy) 
    {
      nLoss = (theNuclearStoppingModel->TheValue(particle, material)) * stepLength;
    }
    }
  
  if (eLoss < 0.0) eLoss = 0.0;
 
  finalT = kineticEnergy - eLoss - nLoss;
 
  if ( EnlossFlucFlag && 0.0 < eLoss && finalT > MinKineticEnergy) 
    {
      
      //  now the electron loss with fluctuation
      eLoss = ElectronicLossFluctuation(particle, couple, eLoss, stepLength) ;
      if (eLoss < 0.0) eLoss = 0.0;
      finalT = kineticEnergy - eLoss - nLoss;
    }
 
  //  stop particle if the kinetic energy <= MinKineticEnergy
  if (finalT*massRatio <= MinKineticEnergy ) 
    {
      
      finalT = 0.0;
      if (!particle->GetDefinition()->GetProcessManager()->GetAtRestProcessVector()->size())
        aParticleChange.ProposeTrackStatus(fStopAndKill);
      else
        aParticleChange.ProposeTrackStatus(fStopButAlive);
    }
 
  aParticleChange.ProposeEnergy( finalT );
  eLoss = kineticEnergy-finalT;
 
  aParticleChange.ProposeLocalEnergyDeposit(eLoss);
  return &aParticleChange ;
}

BuildPhysicsTable()

void G4hImpactIonisation::BuildPhysicsTable ( const G4ParticleDefinition & aParticleType )

virtual

Reimplemented from G4VProcess.

Definition at line 198 of file G4hImpactIonisation.cc.

{
 
  // Verbose print-out
  if(verboseLevel > 0) 
    {
      G4cout << "G4hImpactIonisation::BuildPhysicsTable for "
         << particleDef.GetParticleName()
         << " mass(MeV)= " << particleDef.GetPDGMass()/MeV
         << " charge= " << particleDef.GetPDGCharge()/eplus
         << " type= " << particleDef.GetParticleType()
         << G4endl;
      
      if(verboseLevel > 1) 
    {
      G4ProcessVector* pv = particleDef.GetProcessManager()->GetProcessList();
      
      G4cout << " 0: " << (*pv)[0]->GetProcessName() << " " << (*pv)[0]
         << " 1: " << (*pv)[1]->GetProcessName() << " " << (*pv)[1]
        //        << " 2: " << (*pv)[2]->GetProcessName() << " " << (*pv)[2]
         << G4endl;
      G4cout << "ionModel= " << theIonEffChargeModel
         << " MFPtable= " << theMeanFreePathTable
         << " iniMass= " << initialMass
         << G4endl;
    }
    }
  // End of verbose print-out
 
  if (particleDef.GetParticleType() == "nucleus" &&
      particleDef.GetParticleName() != "GenericIon" &&
      particleDef.GetParticleSubType() == "generic")
    {
 
      G4EnergyLossTables::Register(&particleDef,
                   theDEDXpTable,
                   theRangepTable,
                   theInverseRangepTable,
                   theLabTimepTable,
                   theProperTimepTable,
                   LowestKineticEnergy, HighestKineticEnergy,
                   proton_mass_c2/particleDef.GetPDGMass(),
                   TotBin);
 
      return;
    }
 
  if( !CutsWhereModified() && theLossTable) return;
 
  InitializeParametrisation() ;
  G4Proton* proton = G4Proton::Proton();
  G4AntiProton* antiproton = G4AntiProton::AntiProton();
 
  charge = particleDef.GetPDGCharge() / eplus;
  chargeSquare = charge*charge ;
 
  const G4ProductionCutsTable* theCoupleTable=
    G4ProductionCutsTable::GetProductionCutsTable();
  G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
 
  cutForDelta.clear();
  cutForGamma.clear();
 
  for (G4int j=0; j<numOfCouples; ++j) {
 
    // get material parameters needed for the energy loss calculation
    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
    const G4Material* material= couple->GetMaterial();
 
    // the cut cannot be below lowest limit
    G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[j];
    if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
 
    G4double excEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
 
    tCut = std::max(tCut,excEnergy);
    cutForDelta.push_back(tCut);
 
    // the cut cannot be below lowest limit
    tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
    if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
    tCut = std::max(tCut,minGammaEnergy);
    cutForGamma.push_back(tCut);
  }
 
  if(verboseLevel > 0) {
    G4cout << "Cuts are defined " << G4endl;
  }
 
  if(0.0 < charge)
    {
      {
        BuildLossTable(*proton) ;
 
    //      The following vector has a fixed dimension (see src/G4hImpactLoss.cc for more details)        
    //      It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
    //        G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", particleDef=" << particleDef.GetParticleName() << ", proton=" << proton << ", theLossTable=" << theLossTable << ", CounterOfpProcess=" << CounterOfpProcess << G4endl;
        
        RecorderOfpProcess[CounterOfpProcess] = theLossTable ;
        CounterOfpProcess++;
      }
    } else {
    {
      BuildLossTable(*antiproton) ;
        
      //      The following vector has a fixed dimension (see src/G4hImpactLoss.cc for more details)        
      //      It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
      //        G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", particleDef=" << particleDef.GetParticleName() << ", antiproton=" << antiproton << ", theLossTable=" << theLossTable << ", CounterOfpbarProcess=" << CounterOfpbarProcess << G4endl;
        
      RecorderOfpbarProcess[CounterOfpbarProcess] = theLossTable ;
      CounterOfpbarProcess++;
    }
  }
 
  if(verboseLevel > 0) {
    G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
           << "Loss table is built "
      //       << theLossTable
       << G4endl;
  }
 
  BuildLambdaTable(particleDef) ;
  //  BuildDataForFluorescence(particleDef);
 
  if(verboseLevel > 1) {
    G4cout << (*theMeanFreePathTable) << G4endl;
  }
 
  if(verboseLevel > 0) {
    G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
           << "DEDX table will be built "
      //       << theDEDXpTable << " " << theDEDXpbarTable
      //       << " " << theRangepTable << " " << theRangepbarTable
       << G4endl;
  }
 
  BuildDEDXTable(particleDef) ;
 
  if(verboseLevel > 1) {
    G4cout << (*theDEDXpTable) << G4endl;
  }
 
  if((&particleDef == proton) ||  (&particleDef == antiproton)) PrintInfoDefinition() ;
 
  if(verboseLevel > 0) {
    G4cout << "G4hImpactIonisation::BuildPhysicsTable: end for "
           << particleDef.GetParticleName() << G4endl;
  }
}

ComputeDEDX()

G4double G4hImpactIonisation::ComputeDEDX	(	const G4ParticleDefinition *	aParticle,
		const G4MaterialCutsCouple *	couple,
		G4double	kineticEnergy )

Definition at line 1273 of file G4hImpactIonisation.cc.

{
  const G4Material* material = couple->GetMaterial();
  G4Proton* proton = G4Proton::Proton();
  G4AntiProton* antiproton = G4AntiProton::AntiProton();
  G4double dedx = 0.;
 
  G4double tScaled = kineticEnergy * proton_mass_c2 / (aParticle->GetPDGMass()) ;
  charge  = aParticle->GetPDGCharge() ;
 
  if( charge > 0.) 
    {
      if (tScaled > protonHighEnergy) 
    {
      dedx = G4EnergyLossTables::GetDEDX(proton,tScaled,couple) ;  
    }
      else 
    {
      dedx = ProtonParametrisedDEDX(couple,tScaled) ;
    }
    } 
  else 
    {
      if (tScaled > antiprotonHighEnergy) 
    {
      dedx = G4EnergyLossTables::GetDEDX(antiproton,tScaled,couple);
    } 
      else
    {
      dedx = AntiProtonParametrisedDEDX(couple,tScaled) ;
    }
    }
  dedx *= theIonEffChargeModel->TheValue(aParticle, material, kineticEnergy) ;
  
  return dedx ;
}

GetContinuousStepLimit()

G4double G4hImpactIonisation::GetContinuousStepLimit ( const G4Track & track, G4double previousStepSize, G4double currentMinimumStep, G4double & currentSafety )

inlinevirtual

Implements G4VContinuousDiscreteProcess.

Definition at line 293 of file G4hImpactIonisation.hh.

{
  G4double step = GetConstraints(track.GetDynamicParticle(),track.GetMaterialCutsCouple()) ;
 
  // ---- MGP ---- The following line, taken as is from G4hLowEnergyIonisation,
  // is meaningless: currentMinimumStep is passed by value,
  // therefore any local modification to it has no effect
 
  if ((step > 0.) && (step < currentMinimumStep)) currentMinimumStep = step ;
 
  return step ;
}

GetMeanFreePath()

G4double G4hImpactIonisation::GetMeanFreePath ( const G4Track & track, G4double previousStepSize, enum G4ForceCondition * condition )

virtual

Implements G4hRDEnergyLoss.

Definition at line 596 of file G4hImpactIonisation.cc.

{
  const G4DynamicParticle* dynamicParticle = track.GetDynamicParticle();
  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
  const G4Material* material = couple->GetMaterial();
 
  G4double meanFreePath = DBL_MAX;
  // ---- MGP ---- What is the meaning of the local variable isOutOfRange?
  G4bool isOutRange = false;
 
  *condition = NotForced ;
 
  G4double kineticEnergy = (dynamicParticle->GetKineticEnergy())*initialMass/(dynamicParticle->GetMass());
  charge = dynamicParticle->GetCharge()/eplus;
  chargeSquare = theIonEffChargeModel->TheValue(dynamicParticle, material);
 
  if (kineticEnergy < LowestKineticEnergy) 
    {
      meanFreePath = DBL_MAX;
    }
  else 
    {
      if (kineticEnergy > HighestKineticEnergy)  kineticEnergy = HighestKineticEnergy;
      meanFreePath = (((*theMeanFreePathTable)(couple->GetIndex()))->
              GetValue(kineticEnergy,isOutRange))/chargeSquare;
    }
 
  return meanFreePath ;
}

IsApplicable()

G4bool G4hImpactIonisation::IsApplicable ( const G4ParticleDefinition & particle )

inlinevirtual

Reimplemented from G4VProcess.

Definition at line 310 of file G4hImpactIonisation.hh.

{
  // ---- MGP ---- Better criterion for applicability to be defined;
  // now hard-coded particle mass > 0.1 * proton_mass
 
  return (particle.GetPDGCharge() != 0.0 && particle.GetPDGMass() > CLHEP::proton_mass_c2*0.1);
}

PostStepDoIt()

G4VParticleChange * G4hImpactIonisation::PostStepDoIt ( const G4Track & track, const G4Step & Step )

virtual

Implements G4hRDEnergyLoss.

Definition at line 959 of file G4hImpactIonisation.cc.

{
  // Units are expressed in GEANT4 internal units.
 
  //  std::cout << "----- Calling PostStepDoIt ----- " << std::endl;
 
  aParticleChange.Initialize(track) ;
  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();  
  const G4DynamicParticle* aParticle = track.GetDynamicParticle() ;
  
  // Some kinematics
 
  G4ParticleDefinition* definition = track.GetDefinition();
  G4double mass = definition->GetPDGMass();
  G4double kineticEnergy = aParticle->GetKineticEnergy();
  G4double totalEnergy = kineticEnergy + mass ;
  G4double pSquare = kineticEnergy *( totalEnergy + mass) ;
  G4double eSquare = totalEnergy * totalEnergy;
  G4double betaSquare = pSquare / eSquare;
  G4ThreeVector particleDirection = aParticle->GetMomentumDirection() ;
 
  G4double gamma = kineticEnergy / mass + 1.;
  G4double r = electron_mass_c2 / mass;
  G4double tMax = 2. * electron_mass_c2 *(gamma*gamma - 1.) / (1. + 2.*gamma*r + r*r);
 
  // Validity range for delta electron cross section
  G4double deltaCut = cutForDelta[couple->GetIndex()];
 
  // This should not be a case
  if (deltaCut >= tMax)
    return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
 
  G4double xc = deltaCut / tMax;
  G4double rate = tMax / totalEnergy;
  rate = rate*rate ;
  G4double spin = aParticle->GetDefinition()->GetPDGSpin() ;
 
  // Sampling follows ...
  G4double x = 0.;
  G4double gRej = 0.;
 
  do {
    x = xc / (1. - (1. - xc) * G4UniformRand());
    
    if (0.0 == spin) 
      {
    gRej = 1.0 - betaSquare * x ;
      }
    else if (0.5 == spin) 
      {
    gRej = (1. - betaSquare * x + 0.5 * x*x * rate) / (1. + 0.5 * rate) ;
      } 
    else 
      {
    gRej = (1. - betaSquare * x ) * (1. + x/(3.*xc)) +
      x*x * rate * (1. + 0.5*x/xc) / 3.0 /
      (1. + 1./(3.*xc) + rate *(1.+ 0.5/xc) / 3.);
      }
    
  } while( G4UniformRand() > gRej );
 
  G4double deltaKineticEnergy = x * tMax;
  G4double deltaTotalMomentum = std::sqrt(deltaKineticEnergy * 
                      (deltaKineticEnergy + 2. * electron_mass_c2 ));
  G4double totalMomentum = std::sqrt(pSquare) ;
  G4double cosTheta = deltaKineticEnergy * (totalEnergy + electron_mass_c2) / (deltaTotalMomentum*totalMomentum);
 
  //  protection against cosTheta > 1 or < -1 
  if ( cosTheta < -1. ) cosTheta = -1.;
  if ( cosTheta > 1. ) cosTheta = 1.;
 
  //  direction of the delta electron 
  G4double phi = twopi * G4UniformRand();
  G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
  G4double dirX = sinTheta * std::cos(phi);
  G4double dirY = sinTheta * std::sin(phi);
  G4double dirZ = cosTheta;
 
  G4ThreeVector deltaDirection(dirX,dirY,dirZ);
  deltaDirection.rotateUz(particleDirection);
 
  // create G4DynamicParticle object for delta ray
  G4DynamicParticle* deltaRay = new G4DynamicParticle;
  deltaRay->SetKineticEnergy( deltaKineticEnergy );
  deltaRay->SetMomentumDirection(deltaDirection.x(),
                 deltaDirection.y(),
                 deltaDirection.z());
  deltaRay->SetDefinition(G4Electron::Electron());
 
  // fill aParticleChange
  G4double finalKineticEnergy = kineticEnergy - deltaKineticEnergy;
  std::size_t totalNumber = 1;
 
  // Atomic relaxation
 
  // ---- MGP ---- Temporary limitation: currently PIXE only for incident protons
 
  std::size_t nSecondaries = 0;
  std::vector<G4DynamicParticle*>* secondaryVector = 0;
 
  if (definition == G4Proton::ProtonDefinition())
    {
      const G4Material* material = couple->GetMaterial();
 
      // Lazy initialization of pixeCrossSectionHandler
      if (pixeCrossSectionHandler == 0)
    {
      // Instantiate pixeCrossSectionHandler with selected shell cross section models
      // Ownership of interpolation is transferred to pixeCrossSectionHandler
      G4IInterpolator* interpolation = new G4LogLogInterpolator();
      pixeCrossSectionHandler = 
        new G4PixeCrossSectionHandler(interpolation,modelK,modelL,modelM,eMinPixe,eMaxPixe);
      G4String fileName("proton");
      pixeCrossSectionHandler->LoadShellData(fileName);
      //      pixeCrossSectionHandler->PrintData();
    }
      
      // Select an atom in the current material based on the total shell cross sections
      G4int Z = pixeCrossSectionHandler->SelectRandomAtom(material,kineticEnergy);
      //      std::cout << "G4hImpactIonisation::PostStepDoIt - Z = " << Z << std::endl;
 
      //      G4double microscopicCross = MicroscopicCrossSection(*definition,
      //                              kineticEnergy,
      //                      Z, deltaCut);
  //    G4double crossFromShells = pixeCrossSectionHandler->FindValue(Z,kineticEnergy);
 
      //std::cout << "G4hImpactIonisation: Z= "
      //        << Z
      //        << ", energy = "
      //        << kineticEnergy/MeV
      //        <<" MeV, microscopic = "
      //        << microscopicCross/barn 
      //        << " barn, from shells = "
      //        << crossFromShells/barn
      //        << " barn" 
      //        << std::endl;
 
      // Select a shell in the target atom based on the individual shell cross sections
      G4int shellIndex = pixeCrossSectionHandler->SelectRandomShell(Z,kineticEnergy);
    
      G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
      const G4AtomicShell* atomicShell = transitionManager->Shell(Z,shellIndex);
      G4double bindingEnergy = atomicShell->BindingEnergy();
      
      //     if (verboseLevel > 1) 
      //       {
      //     G4cout << "G4hImpactIonisation::PostStepDoIt - Z = " 
      //         << Z 
      //     << ", shell = " 
      //     << shellIndex
      //     << ", bindingE (keV) = " 
      //     << bindingEnergy/keV
      //     << G4endl;
      //       }
      
      // Generate PIXE if binding energy larger than cut for photons or electrons
      
      G4ParticleDefinition* type = 0;
      
      if (finalKineticEnergy >= bindingEnergy)
    //    && (bindingEnergy >= minGammaEnergy ||  bindingEnergy >= minElectronEnergy) ) 
    {
      // Vacancy in subshell shellIndex; shellId is the subshell identifier in EADL jargon 
      G4int shellId = atomicShell->ShellId();
      // Atomic relaxation: generate secondaries
      secondaryVector = atomicDeexcitation.GenerateParticles(Z, shellId);
 
      // ---- Debug ----
      //std::cout << "ShellId = "
      //    <<shellId << " ---- Atomic relaxation secondaries: ---- " 
      //    << secondaryVector->size()
      //    << std::endl;
 
      // ---- End debug ---
 
      if (secondaryVector != 0) 
        {
          nSecondaries = secondaryVector->size();
          for (std::size_t i = 0; i<nSecondaries; i++) 
        {
          G4DynamicParticle* aSecondary = (*secondaryVector)[i];
          if (aSecondary) 
            {
              G4double e = aSecondary->GetKineticEnergy();
              type = aSecondary->GetDefinition();
 
              // ---- Debug ----
              //if (type == G4Gamma::GammaDefinition())
              //    {           
              //      std::cout << "Z = " << Z 
              //            << ", shell: " << shellId
              //            << ", PIXE photon energy (keV) = " << e/keV 
              //            << std::endl;
              //    }
              // ---- End debug ---
 
              if (e < finalKineticEnergy &&
              ((type == G4Gamma::Gamma() && e > minGammaEnergy ) ||
               (type == G4Electron::Electron() && e > minElectronEnergy ))) 
            {
              // Subtract the energy of the emitted secondary from the primary
              finalKineticEnergy -= e;
              totalNumber++;
              // ---- Debug ----
              //if (type == G4Gamma::GammaDefinition())
              //    {           
              //      std::cout << "Z = " << Z 
              //            << ", shell: " << shellId
              //            << ", PIXE photon energy (keV) = " << e/keV
              //            << std::endl;
              //    }
              // ---- End debug ---
            } 
              else 
            {
              // The atomic relaxation product has energy below the cut
              // ---- Debug ----
              // if (type == G4Gamma::GammaDefinition())
              //    {           
              //      std::cout << "Z = " << Z 
              //
              //            << ", PIXE photon energy = " << e/keV 
              //            << " keV below threshold " << minGammaEnergy/keV << " keV"
              //            << std::endl;
              //    }   
              // ---- End debug ---
     
              delete aSecondary;
              (*secondaryVector)[i] = 0;
            }
            }
        }
        }
    }
    }
 
 
  // Save delta-electrons
 
  G4double eDeposit = 0.;
 
  if (finalKineticEnergy > MinKineticEnergy)
    {
      G4double finalPx = totalMomentum*particleDirection.x() - deltaTotalMomentum*deltaDirection.x();
      G4double finalPy = totalMomentum*particleDirection.y() - deltaTotalMomentum*deltaDirection.y();
      G4double finalPz = totalMomentum*particleDirection.z() - deltaTotalMomentum*deltaDirection.z();
      G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz) ;
      finalPx /= finalMomentum;
      finalPy /= finalMomentum;
      finalPz /= finalMomentum;
 
      aParticleChange.ProposeMomentumDirection( finalPx,finalPy,finalPz );
    }
  else
    {
      eDeposit = finalKineticEnergy;
      finalKineticEnergy = 0.;
      aParticleChange.ProposeMomentumDirection(particleDirection.x(),
                           particleDirection.y(),
                           particleDirection.z());
      if(!aParticle->GetDefinition()->GetProcessManager()->
     GetAtRestProcessVector()->size())
        aParticleChange.ProposeTrackStatus(fStopAndKill);
      else
        aParticleChange.ProposeTrackStatus(fStopButAlive);
    }
 
  aParticleChange.ProposeEnergy(finalKineticEnergy);
  aParticleChange.ProposeLocalEnergyDeposit (eDeposit);
  aParticleChange.SetNumberOfSecondaries((G4int)totalNumber);
  aParticleChange.AddSecondary(deltaRay);
 
  // ---- Debug ----
  //  std::cout << "RDHadronIonisation - finalKineticEnergy (MeV) = " 
  //        << finalKineticEnergy/MeV 
  //        << ", delta KineticEnergy (keV) = " 
  //        << deltaKineticEnergy/keV 
  //        << ", energy deposit (MeV) = "
  //        << eDeposit/MeV
  //        << std::endl;
  // ---- End debug ---
  
  // Save Fluorescence and Auger
 
  if (secondaryVector != 0) 
    { 
      for (std::size_t l = 0; l < nSecondaries; l++) 
    { 
      G4DynamicParticle* secondary = (*secondaryVector)[l];
      if (secondary) aParticleChange.AddSecondary(secondary);
 
      // ---- Debug ----
      //if (secondary != 0) 
      // {
      //   if (secondary->GetDefinition() == G4Gamma::GammaDefinition())
      //    {
      //      G4double eX = secondary->GetKineticEnergy();          
      //      std::cout << " PIXE photon of energy " << eX/keV 
      //            << " keV added to ParticleChange; total number of secondaries is " << totalNumber 
      //            << std::endl;
      //    }
      //}
      // ---- End debug ---   
 
    }
      delete secondaryVector;
    }
 
  return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
}

PrintInfoDefinition()

void G4hImpactIonisation::PrintInfoDefinition

(

)

const

Definition at line 1721 of file G4hImpactIonisation.cc.

{
  G4String comments = "  Knock-on electron cross sections . ";
  comments += "\n        Good description above the mean excitation energy.\n";
  comments += "        Delta ray energy sampled from  differential Xsection.";
 
  G4cout << G4endl << GetProcessName() << ":  " << comments
         << "\n        PhysicsTables from " << LowestKineticEnergy / eV << " eV "
         << " to " << HighestKineticEnergy / TeV << " TeV "
         << " in " << TotBin << " bins."
     << "\n        Electronic stopping power model is  "
     << protonTable
         << "\n        from " << protonLowEnergy / keV << " keV "
         << " to " << protonHighEnergy / MeV << " MeV " << "." << G4endl ;
  G4cout << "\n        Parametrisation model for antiprotons is  "
         << antiprotonTable
         << "\n        from " << antiprotonLowEnergy / keV << " keV "
         << " to " << antiprotonHighEnergy / MeV << " MeV " << "." << G4endl ;
  if(theBarkas){
    G4cout << "        Parametrization of the Barkas effect is switched on."
       << G4endl ;
  }
  if(nStopping) {
    G4cout << "        Nuclear stopping power model is " << theNuclearTable
       << G4endl ;
  }
 
  G4bool printHead = true;
 
  const G4ProductionCutsTable* theCoupleTable=
    G4ProductionCutsTable::GetProductionCutsTable();
  G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
 
  // loop for materials
 
  for (G4int j=0 ; j < numOfCouples; ++j) {
 
    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
    const G4Material* material= couple->GetMaterial();
    G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
    G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
 
    if(excitationEnergy > deltaCutNow) {
      if(printHead) {
        printHead = false ;
 
        G4cout << "           material       min.delta energy(keV) " << G4endl;
        G4cout << G4endl;
      }
 
      G4cout << std::setw(20) << material->GetName()
         << std::setw(15) << excitationEnergy/keV << G4endl;
    }
  }
}

Referenced by BuildPhysicsTable().

SetBarkasOff()

void G4hImpactIonisation::SetBarkasOff ( )

inline

Definition at line 147 of file G4hImpactIonisation.hh.

{theBarkas = false;};

SetBarkasOn()

void G4hImpactIonisation::SetBarkasOn ( )

inline

Definition at line 144 of file G4hImpactIonisation.hh.

{theBarkas = true;};

SetCutForAugerElectrons()

void G4hImpactIonisation::SetCutForAugerElectrons

(

G4double

cut

)

Definition at line 1707 of file G4hImpactIonisation.cc.

{
  minElectronEnergy = cut;
}

SetCutForSecondaryPhotons()

void G4hImpactIonisation::SetCutForSecondaryPhotons

(

G4double

cut

)

Definition at line 1700 of file G4hImpactIonisation.cc.

{
  minGammaEnergy = cut;
}

SetElectronicStoppingPowerModel()

void G4hImpactIonisation::SetElectronicStoppingPowerModel	(	const G4ParticleDefinition *	aParticle,
		const G4String &	dedxTable )

Definition at line 164 of file G4hImpactIonisation.cc.

{
  if (particle->GetPDGCharge() > 0 ) 
    {
      // Positive charge
      SetProtonElectronicStoppingPowerModel(dedxTable) ;
    } 
  else 
    {
      // Antiprotons
      SetAntiProtonElectronicStoppingPowerModel(dedxTable) ;
    }
}

SetHighEnergyForAntiProtonParametrisation()

void G4hImpactIonisation::SetHighEnergyForAntiProtonParametrisation ( G4double energy )

inline

Definition at line 109 of file G4hImpactIonisation.hh.

{antiprotonHighEnergy = energy;} ;

SetHighEnergyForProtonParametrisation()

void G4hImpactIonisation::SetHighEnergyForProtonParametrisation ( G4double energy )

inline

Definition at line 99 of file G4hImpactIonisation.hh.

{protonHighEnergy = energy;} ;

SetLowEnergyForAntiProtonParametrisation()

void G4hImpactIonisation::SetLowEnergyForAntiProtonParametrisation ( G4double energy )

inline

Definition at line 114 of file G4hImpactIonisation.hh.

{antiprotonLowEnergy = energy;} ;

SetLowEnergyForProtonParametrisation()

void G4hImpactIonisation::SetLowEnergyForProtonParametrisation ( G4double energy )

inline

Definition at line 104 of file G4hImpactIonisation.hh.

{protonLowEnergy = energy;} ;

SetNuclearStoppingOff()

void G4hImpactIonisation::SetNuclearStoppingOff ( )

inline

Definition at line 141 of file G4hImpactIonisation.hh.

{nStopping = false;};

SetNuclearStoppingOn()

void G4hImpactIonisation::SetNuclearStoppingOn ( )

inline

Definition at line 138 of file G4hImpactIonisation.hh.

{nStopping = true;};

Referenced by SetNuclearStoppingPowerModel().

SetNuclearStoppingPowerModel()

void G4hImpactIonisation::SetNuclearStoppingPowerModel ( const G4String & dedxTable )

inline

Definition at line 130 of file G4hImpactIonisation.hh.

  {theNuclearTable = dedxTable; SetNuclearStoppingOn();};

SetPixe()

void G4hImpactIonisation::SetPixe ( const G4bool )

inline

Definition at line 150 of file G4hImpactIonisation.hh.

{pixeIsActive = true;};

SetPixeCrossSectionK()

void G4hImpactIonisation::SetPixeCrossSectionK ( const G4String & name )

inline

Definition at line 176 of file G4hImpactIonisation.hh.

{ modelK = name; }

SetPixeCrossSectionL()

void G4hImpactIonisation::SetPixeCrossSectionL ( const G4String & name )

inline

Definition at line 177 of file G4hImpactIonisation.hh.

{ modelL = name; }

SetPixeCrossSectionM()

void G4hImpactIonisation::SetPixeCrossSectionM ( const G4String & name )

inline

Definition at line 178 of file G4hImpactIonisation.hh.

{ modelM = name; }

SetPixeProjectileMaxEnergy()

void G4hImpactIonisation::SetPixeProjectileMaxEnergy ( G4double energy )

inline

Definition at line 180 of file G4hImpactIonisation.hh.

{ eMaxPixe = energy; }

SetPixeProjectileMinEnergy()

void G4hImpactIonisation::SetPixeProjectileMinEnergy ( G4double energy )

inline

Definition at line 179 of file G4hImpactIonisation.hh.

{ eMinPixe = energy; }

---

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

• G4hImpactIonisation.hh
• G4hImpactIonisation.cc