#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 ()

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

G4int	operator== (const G4VProcess &right) const

G4int	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 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

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 previousStepSize)

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

G4VParticleChange *	pParticleChange

G4ParticleChange	aParticleChange

G4double	theNumberOfInteractionLengthLeft

G4double	currentInteractionLength

G4double	theInitialNumberOfInteractionLength

G4String	theProcessName

G4String	thePhysicsTableFileName

G4ProcessType	theProcessType

G4int	theProcessSubType

G4double	thePILfactor

G4bool	enableAtRestDoIt

G4bool	enableAlongStepDoIt

G4bool	enablePostStepDoIt

G4int	verboseLevel

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 75 of file G4hImpactIonisation.hh.

Constructor & Destructor Documentation

G4hImpactIonisation::G4hImpactIonisation ( const G4String & processName = "hImpactIoni" )

Definition at line 83 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 ( )

Definition at line 132 of file G4hImpactIonisation.cc.

References G4PhysicsTable::clearAndDestroy().

 {
   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

void G4hImpactIonisation::ActivateAugerElectronProduction ( G4bool val )

Definition at line 1707 of file G4hImpactIonisation.cc.

References G4AtomicDeexcitation::ActivateAugerElectronProduction().

 {
   atomicDeexcitation.ActivateAugerElectronProduction(val);
 }

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

virtual

Reimplemented from G4VContinuousDiscreteProcess.

Definition at line 718 of file G4hImpactIonisation.cc.

References G4AntiProton::AntiProton(), G4VProcess::aParticleChange, G4hRDEnergyLoss::EnlossFlucFlag, fStopAndKill, fStopButAlive, G4ProcessManager::GetAtRestProcessVector(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4MaterialCutsCouple::GetMaterial(), G4Track::GetMaterialCutsCouple(), G4EnergyLossTables::GetPreciseEnergyFromRange(), G4ParticleDefinition::GetProcessManager(), G4Step::GetStepLength(), G4hRDEnergyLoss::HighestKineticEnergy, G4ParticleChange::Initialize(), G4hRDEnergyLoss::linLossLimit, G4hRDEnergyLoss::MinKineticEnergy, G4ParticleChange::ProposeEnergy(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4VParticleChange::ProposeTrackStatus(), G4Proton::Proton(), python.hepunit::proton_mass_c2, G4ProcessVector::size(), and G4VLowEnergyModel::TheValue().

 {
   // 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 ;
 }

void G4hImpactIonisation::BuildPhysicsTable ( const G4ParticleDefinition & aParticleType )

virtual

Reimplemented from G4VProcess.

Definition at line 191 of file G4hImpactIonisation.cc.

References G4AntiProton::AntiProton(), G4hRDEnergyLoss::BuildDEDXTable(), G4hRDEnergyLoss::CounterOfpbarProcess, G4hRDEnergyLoss::CounterOfpProcess, G4hRDEnergyLoss::CutsWhereModified(), python.hepunit::eplus, G4cout, G4endl, G4ProductionCutsTable::GetEnergyCutsVector(), G4Material::GetIonisation(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4IonisParamMat::GetMeanExcitationEnergy(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetParticleType(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4ProcessManager::GetProcessList(), G4ParticleDefinition::GetProcessManager(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), G4hRDEnergyLoss::HighestKineticEnergy, G4hRDEnergyLoss::LowestKineticEnergy, eplot::material, G4INCL::Math::max(), python.hepunit::MeV, PrintInfoDefinition(), G4InuclParticleNames::proton, G4Proton::Proton(), python.hepunit::proton_mass_c2, G4hRDEnergyLoss::RecorderOfpbarProcess, G4hRDEnergyLoss::RecorderOfpProcess, G4EnergyLossTables::Register(), G4hRDEnergyLoss::theDEDXpTable, G4hRDEnergyLoss::theInverseRangepTable, G4hRDEnergyLoss::theLabTimepTable, G4hRDEnergyLoss::theLossTable, G4hRDEnergyLoss::theProperTimepTable, G4hRDEnergyLoss::theRangepTable, G4hRDEnergyLoss::TotBin, and G4VProcess::verboseLevel.

 {
 
   // 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();
   size_t numOfCouples = theCoupleTable->GetTableSize();
 
   cutForDelta.clear();
   cutForGamma.clear();
 
   for (size_t 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;
   }
 }

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

Definition at line 1266 of file G4hImpactIonisation.cc.

References G4AntiProton::AntiProton(), G4EnergyLossTables::GetDEDX(), G4MaterialCutsCouple::GetMaterial(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4Proton::Proton(), python.hepunit::proton_mass_c2, and G4VLowEnergyModel::TheValue().

 {
   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 ;
 }

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

inlinevirtual

Implements G4VContinuousDiscreteProcess.

Definition at line 294 of file G4hImpactIonisation.hh.

References G4Track::GetDynamicParticle(), and G4Track::GetMaterialCutsCouple().

 {
   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 ;
 }

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

virtual

Implements G4hRDEnergyLoss.

Definition at line 589 of file G4hImpactIonisation.cc.

References DBL_MAX, python.hepunit::eplus, G4DynamicParticle::GetCharge(), G4Track::GetDynamicParticle(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4MaterialCutsCouple::GetMaterial(), G4Track::GetMaterialCutsCouple(), G4hRDEnergyLoss::HighestKineticEnergy, G4hRDEnergyLoss::LowestKineticEnergy, NotForced, and G4VLowEnergyModel::TheValue().

 {
   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 ;
 }

G4bool G4hImpactIonisation::IsApplicable ( const G4ParticleDefinition & particle )

inlinevirtual

Reimplemented from G4VProcess.

Definition at line 311 of file G4hImpactIonisation.hh.

References G4ParticleDefinition::GetPDGCharge(), and G4ParticleDefinition::GetPDGMass().

 {
   // ---- 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);
 }

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

virtual

Implements G4hRDEnergyLoss.

Definition at line 952 of file G4hImpactIonisation.cc.

References G4ParticleChange::AddSecondary(), G4VProcess::aParticleChange, G4InuclSpecialFunctions::bindingEnergy(), G4AtomicShell::BindingEnergy(), G4Electron::Electron(), python.hepunit::electron_mass_c2, fStopAndKill, fStopButAlive, G4UniformRand, G4Gamma::Gamma(), G4AtomicDeexcitation::GenerateParticles(), G4Track::GetDefinition(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4MaterialCutsCouple::GetMaterial(), G4Track::GetMaterialCutsCouple(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGMass(), G4ParticleDefinition::GetPDGSpin(), G4ParticleDefinition::GetProcessManager(), G4ParticleChange::Initialize(), G4AtomicTransitionManager::Instance(), G4PixeCrossSectionHandler::LoadShellData(), G4hRDEnergyLoss::MinKineticEnergy, G4VContinuousDiscreteProcess::PostStepDoIt(), G4ParticleChange::ProposeEnergy(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4Proton::ProtonDefinition(), CLHEP::Hep3Vector::rotateUz(), G4PixeCrossSectionHandler::SelectRandomAtom(), G4PixeCrossSectionHandler::SelectRandomShell(), G4DynamicParticle::SetDefinition(), G4DynamicParticle::SetKineticEnergy(), G4DynamicParticle::SetMomentumDirection(), G4VParticleChange::SetNumberOfSecondaries(), G4AtomicTransitionManager::Shell(), G4AtomicShell::ShellId(), python.hepunit::twopi, test::x, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

 {
   // 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;
   size_t totalNumber = 1;
 
   // Atomic relaxation
 
   // ---- MGP ---- Temporary limitation: currently PIXE only for incident protons
 
   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 (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(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 (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);
 }

void G4hImpactIonisation::PrintInfoDefinition ( ) const

Definition at line 1714 of file G4hImpactIonisation.cc.

References python.hepunit::eV, G4cout, G4endl, G4MaterialCutsCouple::GetIndex(), G4Material::GetIonisation(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4IonisParamMat::GetMeanExcitationEnergy(), G4Material::GetName(), G4VProcess::GetProcessName(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), G4hRDEnergyLoss::HighestKineticEnergy, python.hepunit::keV, G4hRDEnergyLoss::LowestKineticEnergy, python.hepunit::MeV, and python.hepunit::TeV.

Referenced by BuildPhysicsTable().

 {
   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();
   size_t numOfCouples = theCoupleTable->GetTableSize();
 
   // loop for materials
 
   for (size_t 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;
     }
   }
 }

void G4hImpactIonisation::SetBarkasOff ( )

inline

Definition at line 148 of file G4hImpactIonisation.hh.

148 {theBarkas = false;};

void G4hImpactIonisation::SetBarkasOn ( )

inline

Definition at line 145 of file G4hImpactIonisation.hh.

145 {theBarkas = true;};

void G4hImpactIonisation::SetCutForAugerElectrons ( G4double cut )

Definition at line 1700 of file G4hImpactIonisation.cc.

Referenced by eRositaPhysicsList::ConstructEM().

 {
   minElectronEnergy = cut;
 }

void G4hImpactIonisation::SetCutForSecondaryPhotons ( G4double cut )

Definition at line 1693 of file G4hImpactIonisation.cc.

Referenced by eRositaPhysicsList::ConstructEM().

 {
   minGammaEnergy = cut;
 }

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

Definition at line 157 of file G4hImpactIonisation.cc.

References G4ParticleDefinition::GetPDGCharge().

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