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)

---

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(), G4InuclParticleNames::proton, 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(), 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, PrintInfoDefinition(), G4Proton::Proton(), G4InuclParticleNames::proton, 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(), G4InuclParticleNames::proton, 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
	)			[inline, virtual]

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, G4DynamicParticle::GetCharge(), G4Track::GetDynamicParticle(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4MaterialCutsCouple::GetMaterial(), G4Track::GetMaterialCutsCouple(), G4hRDEnergyLoss::HighestKineticEnergy, 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 &

)

[inline, virtual]

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, G4AtomicShell::BindingEnergy(), G4InuclSpecialFunctions::bindingEnergy(), G4Electron::Electron(), fStopAndKill, fStopButAlive, G4UniformRand, G4Gamma::Gamma(), G4AtomicDeexcitation::GenerateParticles(), G4DynamicParticle::GetDefinition(), G4Track::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(), G4PixeCrossSectionHandler::SelectRandomAtom(), G4PixeCrossSectionHandler::SelectRandomShell(), G4DynamicParticle::SetDefinition(), G4DynamicParticle::SetKineticEnergy(), G4DynamicParticle::SetMomentumDirection(), G4VParticleChange::SetNumberOfSecondaries(), G4AtomicTransitionManager::Shell(), and G4AtomicShell::ShellId().

{
  // 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 G4cout, G4endl, G4MaterialCutsCouple::GetIndex(), G4Material::GetIonisation(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4IonisParamMat::GetMeanExcitationEnergy(), G4Material::GetName(), G4VProcess::GetProcessName(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), and G4hRDEnergyLoss::HighestKineticEnergy.

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.

{theBarkas = false;};

void G4hImpactIonisation::SetBarkasOn

(

)

[inline]

Definition at line 145 of file G4hImpactIonisation.hh.

{theBarkas = true;};

void G4hImpactIonisation::SetCutForAugerElectrons

(

G4double

cut

)

Definition at line 1700 of file G4hImpactIonisation.cc.

{
  minElectronEnergy = cut;
}

void G4hImpactIonisation::SetCutForSecondaryPhotons

(

G4double

cut

)

Definition at line 1693 of file G4hImpactIonisation.cc.

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

void G4hImpactIonisation::SetHighEnergyForAntiProtonParametrisation

(

G4double

energy

)

[inline]

Definition at line 110 of file G4hImpactIonisation.hh.

{antiprotonHighEnergy = energy;} ;

void G4hImpactIonisation::SetHighEnergyForProtonParametrisation

(

G4double

energy

)

[inline]

Definition at line 100 of file G4hImpactIonisation.hh.

{protonHighEnergy = energy;} ;

void G4hImpactIonisation::SetLowEnergyForAntiProtonParametrisation

(

G4double

energy

)

[inline]

Definition at line 115 of file G4hImpactIonisation.hh.

{antiprotonLowEnergy = energy;} ;

void G4hImpactIonisation::SetLowEnergyForProtonParametrisation

(

G4double

energy

)

[inline]

Definition at line 105 of file G4hImpactIonisation.hh.

{protonLowEnergy = energy;} ;

void G4hImpactIonisation::SetNuclearStoppingOff

(

)

[inline]

Definition at line 142 of file G4hImpactIonisation.hh.

{nStopping = false;};

void G4hImpactIonisation::SetNuclearStoppingOn

(

)

[inline]

Definition at line 139 of file G4hImpactIonisation.hh.

Referenced by SetNuclearStoppingPowerModel().

{nStopping = true;};

void G4hImpactIonisation::SetNuclearStoppingPowerModel

(

const G4String &

dedxTable

)

[inline]

Definition at line 131 of file G4hImpactIonisation.hh.

References SetNuclearStoppingOn().

  {theNuclearTable = dedxTable; SetNuclearStoppingOn();};

void G4hImpactIonisation::SetPixe

(

const

G4bool

)

[inline]

Definition at line 151 of file G4hImpactIonisation.hh.

{pixeIsActive = true;};

void G4hImpactIonisation::SetPixeCrossSectionK

(

const G4String &

name

)

[inline]

Definition at line 177 of file G4hImpactIonisation.hh.

{ modelK = name; }

void G4hImpactIonisation::SetPixeCrossSectionL

(

const G4String &

name

)

[inline]

Definition at line 178 of file G4hImpactIonisation.hh.

{ modelL = name; }

void G4hImpactIonisation::SetPixeCrossSectionM

(

const G4String &

name

)

[inline]

Definition at line 179 of file G4hImpactIonisation.hh.

{ modelM = name; }

void G4hImpactIonisation::SetPixeProjectileMaxEnergy

(

G4double

energy

)

[inline]

Definition at line 181 of file G4hImpactIonisation.hh.

{ eMaxPixe = energy; }

void G4hImpactIonisation::SetPixeProjectileMinEnergy

(

G4double

energy

)

[inline]

Definition at line 180 of file G4hImpactIonisation.hh.

{ eMinPixe = energy; }

---

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

• G4hImpactIonisation.hh
• G4hImpactIonisation.cc

---