G4AdjointhIonisationModel Class Reference

#include <G4AdjointhIonisationModel.hh>

Inheritance diagram for G4AdjointhIonisationModel:

Public Member Functions
	G4AdjointhIonisationModel (G4ParticleDefinition *projectileDefinition)
virtual	~G4AdjointhIonisationModel ()
virtual void	SampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
void	RapidSampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
virtual G4double	DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
virtual G4double	AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
virtual G4double	GetSecondAdjEnergyMaxForScatProjToProjCase (G4double PrimAdjEnergy)
virtual G4double	GetSecondAdjEnergyMinForScatProjToProjCase (G4double PrimAdjEnergy, G4double Tcut=0)
virtual G4double	GetSecondAdjEnergyMaxForProdToProjCase (G4double PrimAdjEnergy)
virtual G4double	GetSecondAdjEnergyMinForProdToProjCase (G4double PrimAdjEnergy)
Public Member Functions inherited from G4VEmAdjointModel
	G4VEmAdjointModel (const G4String &nam)
virtual	~G4VEmAdjointModel ()
virtual void	SampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)=0
virtual G4double	AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
virtual G4double	GetAdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
virtual G4double	DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
virtual G4double	DiffCrossSectionPerAtomPrimToScatPrim (G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
virtual G4double	DiffCrossSectionPerVolumePrimToSecond (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
virtual G4double	DiffCrossSectionPerVolumePrimToScatPrim (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
virtual G4double	GetSecondAdjEnergyMaxForScatProjToProjCase (G4double PrimAdjEnergy)
virtual G4double	GetSecondAdjEnergyMinForScatProjToProjCase (G4double PrimAdjEnergy, G4double Tcut=0)
virtual G4double	GetSecondAdjEnergyMaxForProdToProjCase (G4double PrimAdjEnergy)
virtual G4double	GetSecondAdjEnergyMinForProdToProjCase (G4double PrimAdjEnergy)
void	DefineCurrentMaterial (const G4MaterialCutsCouple *couple)
std::vector< std::vector< double > * >	ComputeAdjointCrossSectionVectorPerAtomForSecond (G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
std::vector< std::vector< double > * >	ComputeAdjointCrossSectionVectorPerAtomForScatProj (G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
std::vector< std::vector< double > * >	ComputeAdjointCrossSectionVectorPerVolumeForSecond (G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
std::vector< std::vector< double > * >	ComputeAdjointCrossSectionVectorPerVolumeForScatProj (G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
void	SetCSMatrices (std::vector< G4AdjointCSMatrix * > *Vec1CSMatrix, std::vector< G4AdjointCSMatrix * > *Vec2CSMatrix)
G4ParticleDefinition *	GetAdjointEquivalentOfDirectPrimaryParticleDefinition ()
G4ParticleDefinition *	GetAdjointEquivalentOfDirectSecondaryParticleDefinition ()
G4double	GetHighEnergyLimit ()
G4double	GetLowEnergyLimit ()
void	SetHighEnergyLimit (G4double aVal)
void	SetLowEnergyLimit (G4double aVal)
void	DefineDirectEMModel (G4VEmModel *aModel)
void	SetAdjointEquivalentOfDirectPrimaryParticleDefinition (G4ParticleDefinition *aPart)
void	SetAdjointEquivalentOfDirectSecondaryParticleDefinition (G4ParticleDefinition *aPart)
void	SetSecondPartOfSameType (G4bool aBool)
G4bool	GetSecondPartOfSameType ()
void	SetUseMatrix (G4bool aBool)
void	SetUseMatrixPerElement (G4bool aBool)
void	SetUseOnlyOneMatrixForAllElements (G4bool aBool)
void	SetApplyCutInRange (G4bool aBool)
G4bool	GetUseMatrix ()
G4bool	GetUseMatrixPerElement ()
G4bool	GetUseOnlyOneMatrixForAllElements ()
G4bool	GetApplyCutInRange ()
G4String	GetName ()
virtual void	SetCSBiasingFactor (G4double aVal)

Additional Inherited Members
Protected Member Functions inherited from G4VEmAdjointModel
G4double	DiffCrossSectionFunction1 (G4double kinEnergyProj)
G4double	DiffCrossSectionFunction2 (G4double kinEnergyProj)
G4double	DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj (G4double EkinProd)
G4double	SampleAdjSecEnergyFromCSMatrix (size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
G4double	SampleAdjSecEnergyFromCSMatrix (G4double prim_energy, G4bool IsScatProjToProjCase)
void	SelectCSMatrix (G4bool IsScatProjToProjCase)
virtual G4double	SampleAdjSecEnergyFromDiffCrossSectionPerAtom (G4double prim_energy, G4bool IsScatProjToProjCase)
virtual void	CorrectPostStepWeight (G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
Protected Attributes inherited from G4VEmAdjointModel
G4VEmModel *	theDirectEMModel
G4VParticleChange *	pParticleChange
const G4String	name
G4int	ASelectedNucleus
G4int	ZSelectedNucleus
G4Material *	SelectedMaterial
G4double	kinEnergyProdForIntegration
G4double	kinEnergyScatProjForIntegration
G4double	kinEnergyProjForIntegration
std::vector< G4AdjointCSMatrix * > *	pOnCSMatrixForProdToProjBackwardScattering
std::vector< G4AdjointCSMatrix * > *	pOnCSMatrixForScatProjToProjBackwardScattering
std::vector< G4double >	CS_Vs_ElementForScatProjToProjCase
std::vector< G4double >	CS_Vs_ElementForProdToProjCase
G4double	lastCS
G4double	lastAdjointCSForScatProjToProjCase
G4double	lastAdjointCSForProdToProjCase
G4ParticleDefinition *	theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition *	theAdjEquivOfDirectSecondPartDef
G4ParticleDefinition *	theDirectPrimaryPartDef
G4bool	second_part_of_same_type
G4double	preStepEnergy
G4Material *	currentMaterial
G4MaterialCutsCouple *	currentCouple
size_t	currentMaterialIndex
size_t	currentCoupleIndex
G4double	currentTcutForDirectPrim
G4double	currentTcutForDirectSecond
G4bool	ApplyCutInRange
G4double	mass_ratio_product
G4double	mass_ratio_projectile
G4double	HighEnergyLimit
G4double	LowEnergyLimit
G4double	CS_biasing_factor
G4bool	UseMatrix
G4bool	UseMatrixPerElement
G4bool	UseOnlyOneMatrixForAllElements
size_t	indexOfUsedCrossSectionMatrix
size_t	model_index

Detailed Description

Definition at line 72 of file G4AdjointhIonisationModel.hh.

Constructor & Destructor Documentation

G4AdjointhIonisationModel()

G4AdjointhIonisationModel::G4AdjointhIonisationModel

(

G4ParticleDefinition *

projectileDefinition

)

Definition at line 46 of file G4AdjointhIonisationModel.cc.

                                                                                              :
  G4VEmAdjointModel("Adjoint_hIonisation")
{ 
 
 
 
  UseMatrix =true;
  UseMatrixPerElement = true;
  ApplyCutInRange = true;
  UseOnlyOneMatrixForAllElements = true; 
  CS_biasing_factor =1.;
  second_part_of_same_type =false;
  
  //The direct EM Modfel is taken has BetheBloch it is only used for the computation 
  // of the differential cross section.
  //The Bragg model could be used  as an alternative as  it offers the same differential cross section
  
  theDirectEMModel = new G4BetheBlochModel(projectileDefinition);
  theBraggDirectEMModel = new G4BraggModel(projectileDefinition); 
  theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
  
  theDirectPrimaryPartDef = projectileDefinition;
  theAdjEquivOfDirectPrimPartDef = 0;
  if (projectileDefinition == G4Proton::Proton()) {
    theAdjEquivOfDirectPrimPartDef = G4AdjointProton::AdjointProton();
    
  }
  
  DefineProjectileProperty();
}

~G4AdjointhIonisationModel()

G4AdjointhIonisationModel::~G4AdjointhIonisationModel ( )

virtual

Definition at line 78 of file G4AdjointhIonisationModel.cc.

{;}

Member Function Documentation

AdjointCrossSection()

G4double G4AdjointhIonisationModel::AdjointCrossSection ( const G4MaterialCutsCouple *  aCouple, G4double  primEnergy, G4bool  IsScatProjToProjCase  )

virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 438 of file G4AdjointhIonisationModel.cc.

{ 
  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
  DefineCurrentMaterial(aCouple);
    
  
  G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
  
  if (!IsScatProjToProjCase ){
    G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
    G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
    if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
        Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
    } 
    else Cross=0.;
        
    
 
    
    
    
  }
  else {
    G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
    G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
    G4double diff1=Emin_proj-primEnergy;
    G4double diff2=Emax_proj-primEnergy;
    G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
    //G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
    G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy;
    Cross*=(t1+t2);
 
    
  }
  lastCS =Cross;
  return Cross; 
}

DiffCrossSectionPerAtomPrimToSecond()

G4double G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond ( G4double  kinEnergyProj, G4double  kinEnergyProd, G4double  Z, G4double  A = 0.  )

virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 306 of file G4AdjointhIonisationModel.cc.

{//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
  
 
 
 G4double dSigmadEprod=0;
 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
 
 
 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile 
    G4double Tmax=kinEnergyProj;
    
        G4double E1=kinEnergyProd;
    G4double E2=kinEnergyProd*1.000001;
    G4double dE=(E2-E1);
    G4double sigma1,sigma2;
    if (kinEnergyProj >2.*MeV){
        sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
        sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
        }
    else {
        sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
        sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
    }
    
    
    dSigmadEprod=(sigma1-sigma2)/dE;
    if (dSigmadEprod>1.) {
        G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
        G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
        G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
        
    }
    
     
    
     //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
     //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
     //to test the rejection of a secondary
     //-------------------------
     
     //Source code taken from    G4BetheBlochModel::SampleSecondaries
     
     G4double deltaKinEnergy = kinEnergyProd;
     
     //Part of the taken code
     //----------------------
     
     
     
     // projectile formfactor - suppresion of high energy
         // delta-electron production at high energy
     G4double x = formfact*deltaKinEnergy;
          if(x > 1.e-6) {
        
        
        G4double totEnergy     = kinEnergyProj + mass;
            G4double etot2         = totEnergy*totEnergy;
        G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
        G4double f;
        G4double f1 = 0.0;
        f = 1.0 - beta2*deltaKinEnergy/Tmax;
            if( 0.5 == spin ) {
                f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
                f += f1;
            }
        G4double x1 = 1.0 + x;
            G4double gg  = 1.0/(x1*x1);
            if( 0.5 == spin ) {
                G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
                gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
            }
            if(gg > 1.0) {
                G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
                    << G4endl;
            gg=1.;
        }
        //G4cout<<"gg"<<gg<<G4endl;
        dSigmadEprod*=gg;       
         }
     
  }
 
 return dSigmadEprod;   
}

Referenced by RapidSampleSecondaries().

GetSecondAdjEnergyMaxForProdToProjCase()

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase ( G4double  PrimAdjEnergy )

virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 491 of file G4AdjointhIonisationModel.cc.

{ return HighEnergyLimit;
}

Referenced by AdjointCrossSection(), DiffCrossSectionPerAtomPrimToSecond(), and RapidSampleSecondaries().

GetSecondAdjEnergyMaxForScatProjToProjCase()

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase ( G4double  PrimAdjEnergy )

virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 479 of file G4AdjointhIonisationModel.cc.

{ 
  G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
  return Tmax;
}

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

GetSecondAdjEnergyMinForProdToProjCase()

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase ( G4double  PrimAdjEnergy )

virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 496 of file G4AdjointhIonisationModel.cc.

{  G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
   return Tmin;
}

Referenced by AdjointCrossSection(), DiffCrossSectionPerAtomPrimToSecond(), and RapidSampleSecondaries().

GetSecondAdjEnergyMinForScatProjToProjCase()

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase ( G4double  PrimAdjEnergy, G4double  Tcut = 0  )

virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 486 of file G4AdjointhIonisationModel.cc.

{ return PrimAdjEnergy+Tcut;
}

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

RapidSampleSecondaries()

void G4AdjointhIonisationModel::RapidSampleSecondaries	(	const G4Track &	aTrack,
		G4bool	IsScatProjToProjCase,
		G4ParticleChange *	fParticleChange
	)

Definition at line 160 of file G4AdjointhIonisationModel.cc.

{ 
 
 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
 
 
 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
 
 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
    return;
 }
  
 G4double projectileKinEnergy =0.;
 G4double eEnergy=0.;
 G4double newCS=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
 if (!IsScatProjToProjCase){//1/E^2 distribution
    
    eEnergy=adjointPrimKinEnergy;
    G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
        G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
    if (Emin>=Emax) return;
    G4double a=1./Emax;
    G4double b=1./Emin;
    newCS=newCS*(b-a)/eEnergy;
    
    projectileKinEnergy =1./(b- (b-a)*G4UniformRand()); 
    
    
 }
 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
    G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
    if (Emin>=Emax) return;
    G4double diff1=Emin-adjointPrimKinEnergy;
    G4double diff2=Emax-adjointPrimKinEnergy;
    
    G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
    G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
    /*G4double f31=diff1/Emin;
    G4double f32=diff2/Emax/f31;*/
    G4double t3=2.*std::log(Emax/Emin);
    G4double sum_t=t1+t2+t3;
    newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
    G4double t=G4UniformRand()*sum_t;
    if (t <=t1 ){
        G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ;
        projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
        
    }
    else if (t <=t2 )  {
        G4double q= G4UniformRand()*t2/adjointPrimKinEnergy;
        projectileKinEnergy =1./(1./Emin-q);
    }
    else {  
        projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
        
    }
    eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
    
    
 }
 
 
 
 G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy; 
  
  
                            
 //Weight correction
 //-----------------------
 //First w_corr is set to the ratio between adjoint total CS and fwd total CS
 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
 
 //G4cout<<w_corr<<G4endl;
 w_corr*=newCS/lastCS;
 //G4cout<<w_corr<<G4endl;
 //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the one consistent with the direct model
 //Here we consider the true  diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS
 
 G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1);
 w_corr*=diffCS/diffCS_perAtom_Used;
 //G4cout<<w_corr<<G4endl;
       
 G4double new_weight = aTrack.GetWeight()*w_corr;
 fParticleChange->SetParentWeightByProcess(false);
 fParticleChange->SetSecondaryWeightByProcess(false);
 fParticleChange->ProposeParentWeight(new_weight);
 
 
 
 
 //Kinematic: 
 //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest  and gives
 // him part of its  energy
 //----------------------------------------------------------------------------------------
 
 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;   
 
 
 
 //Companion
 //-----------
 G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
 if (IsScatProjToProjCase) {
    companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
 } 
 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;    
 
 
 //Projectile momentum
 //--------------------
 G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
 G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
 G4double phi =G4UniformRand()*2.*3.1415926;
 G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
 projectileMomentum.rotateUz(dir_parallel);
 
 
 
 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
    fParticleChange->ProposeTrackStatus(fStopAndKill);
    fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
    //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
 }
 else {
    fParticleChange->ProposeEnergy(projectileKinEnergy);
    fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
 }  
        
 
 
 
 
 
 
} 

Referenced by SampleSecondaries().

SampleSecondaries()

void G4AdjointhIonisationModel::SampleSecondaries ( const G4Track &  aTrack, G4bool  IsScatProjToProjCase, G4ParticleChange *  fParticleChange  )

virtual

Implements G4VEmAdjointModel.

Definition at line 84 of file G4AdjointhIonisationModel.cc.

{ 
 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
 
 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
  
 //Elastic inverse scattering 
 //---------------------------------------------------------
 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
 
 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
    return;
 }
 
 //Sample secondary energy
 //-----------------------
 G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
 CorrectPostStepWeight(fParticleChange, 
                  aTrack.GetWeight(), 
                  adjointPrimKinEnergy,
                  projectileKinEnergy,
                  IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
 
         
 //Kinematic: 
 //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest  and gives
 // him part of its  energy
 //----------------------------------------------------------------------------------------
 
 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;   
 
 
 
 //Companion
 //-----------
 G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
 if (IsScatProjToProjCase) {
    companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
 } 
 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;    
 
 
 //Projectile momentum
 //--------------------
 G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
 G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
 G4double phi =G4UniformRand()*2.*3.1415926;
 G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
 projectileMomentum.rotateUz(dir_parallel);
 
 
 
 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
    fParticleChange->ProposeTrackStatus(fStopAndKill);
    fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
    //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
 }
 else {
    fParticleChange->ProposeEnergy(projectileKinEnergy);
    fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
 }  
        
 
  
 
}

---

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

• G4AdjointhIonisationModel.hh
• G4AdjointhIonisationModel.cc