Geant4: G4AdjointIonIonisationModel Class Reference

G4AdjointIonIonisationModel::G4AdjointIonIonisationModel ( )

Definition at line 47 of file G4AdjointIonIonisationModel.cc.

References G4AdjointElectron::AdjointElectron(), G4VEmAdjointModel::ApplyCutInRange, G4VEmAdjointModel::CS_biasing_factor, G4GenericIon::GenericIon(), G4VEmAdjointModel::second_part_of_same_type, G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef, G4VEmAdjointModel::theDirectPrimaryPartDef, G4VEmAdjointModel::UseMatrix, G4VEmAdjointModel::UseMatrixPerElement, and G4VEmAdjointModel::UseOnlyOneMatrixForAllElements.

00047                                                         :
00048   G4VEmAdjointModel("Adjoint_IonIonisation")
00049 { 
00050 
00051 
00052   UseMatrix =true;
00053   UseMatrixPerElement = true;
00054   ApplyCutInRange = true;
00055   UseOnlyOneMatrixForAllElements = true;
00056   CS_biasing_factor =1.;
00057   second_part_of_same_type =false;
00058   use_only_bragg = false; // for the Ion ionisation using the parametrised table model the cross sections and the sample of secondaries is done
00059                                 // as in the BraggIonModel, Therefore the use of this flag;  
00060   
00061   //The direct EM Model is taken has BetheBloch it is only used for the computation 
00062   // of the differential cross section.
00063   //The Bragg model could be used  as an alternative as  it offers the same differential cross section
00064   
00065   theBetheBlochDirectEMModel = new G4BetheBlochModel(G4GenericIon::GenericIon());
00066   theBraggIonDirectEMModel = new G4BraggIonModel(G4GenericIon::GenericIon()); 
00067   theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
00068   theDirectPrimaryPartDef =0;
00069   theAdjEquivOfDirectPrimPartDef =0;
00070  /* theDirectPrimaryPartDef =fwd_ion;
00071   theAdjEquivOfDirectPrimPartDef =adj_ion;
00072  
00073   DefineProjectileProperty();*/
00074 
00075 
00076 
00077 
00078 }

G4AdjointIonIonisationModel::~G4AdjointIonIonisationModel ( ) [virtual]

Definition at line 81 of file G4AdjointIonIonisationModel.cc.

00082 {;}

void G4AdjointIonIonisationModel::CorrectPostStepWeight	(	G4ParticleChange *	fParticleChange,
		G4double	old_weight,
		G4double	adjointPrimKinEnergy,
		G4double	projectileKinEnergy,
		G4bool	IsScatProjToProjCase
	)			`[virtual]`

Reimplemented from G4VEmAdjointModel.

Definition at line 266 of file G4AdjointIonIonisationModel.cc.

References G4VEmModel::ComputeCrossSectionPerAtom(), G4VEmAdjointModel::CS_biasing_factor, G4VEmAdjointModel::currentMaterial, G4VEmAdjointModel::currentTcutForDirectSecond, G4GenericIon::GenericIon(), G4AdjointCSManager::GetAdjointCSManager(), G4VEmModel::GetChargeSquareRatio(), G4AdjointCSManager::GetPostStepWeightCorrection(), G4VParticleChange::ProposeParentWeight(), G4VParticleChange::SetParentWeightByProcess(), G4VParticleChange::SetSecondaryWeightByProcess(), G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.

Referenced by SampleSecondaries().

00268 {
00269  //It is needed because the direct cross section used to compute the differential cross section is not the one used in 
00270  // the direct model where the GenericIon stuff is considered with correction of effective charge.  In the direct model the samnepl of secondaries does
00271  // not reflect the integral cross section. The integral fwd cross section that we used to compute the differential CS
00272  // match the sample of secondaries in the forward case despite the fact that its is not the same total CS than in the FWD case. For this reasion an extra
00273  // weight correction is needed at the end.
00274  
00275 
00276  G4double new_weight=old_weight;
00277  
00278  //the correction of CS due to the problem explained above
00279  G4double kinEnergyProjScaled = massRatio*projectileKinEnergy;
00280  theDirectEMModel =theBraggIonDirectEMModel;
00281  if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
00282  G4double UsedFwdCS=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,projectileKinEnergy,1,1 ,currentTcutForDirectSecond,1.e20);
00283  G4double chargeSqRatio =1.;
00284  if (chargeSquare>1.) chargeSqRatio =  theDirectEMModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,projectileKinEnergy);
00285  G4double  CorrectFwdCS = chargeSqRatio*theDirectEMModel->ComputeCrossSectionPerAtom(G4GenericIon::GenericIon(),kinEnergyProjScaled,1,1 ,currentTcutForDirectSecond,1.e20);
00286  if (UsedFwdCS >0)  new_weight*= CorrectFwdCS/UsedFwdCS;//May be some check is needed if UsedFwdCS ==0 probably that then we should avoid a secondary to be produced, 
00287  
00288  
00289  //additional CS crorrection  needed for cross section biasing in general. 
00290  //May be wrong for ions!!! Most of the time not used!
00291   G4double w_corr =1./CS_biasing_factor;
00292   w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
00293   new_weight*=w_corr;
00294  
00295   new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
00296  
00297  fParticleChange->SetParentWeightByProcess(false);
00298  fParticleChange->SetSecondaryWeightByProcess(false);
00299  fParticleChange->ProposeParentWeight(new_weight);
00300 }

G4double G4AdjointIonIonisationModel::DiffCrossSectionPerAtomPrimToSecond	(	G4double	kinEnergyProj,
		G4double	kinEnergyProd,
		G4double	Z,
		G4double	A = `0.`
	)			`[virtual]`

Reimplemented from G4VEmAdjointModel.

Definition at line 159 of file G4AdjointIonIonisationModel.cc.

References G4VEmModel::ComputeCrossSectionPerAtom(), G4cout, G4endl, GetSecondAdjEnergyMaxForProdToProjCase(), GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.

00164 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
00165   
00166 
00167  
00168  G4double dSigmadEprod=0;
00169  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00170  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00171  
00172  G4double kinEnergyProjScaled = massRatio*kinEnergyProj;
00173  
00174  
00175  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile 
00176         G4double Tmax=kinEnergyProj;
00177         
00178         G4double E1=kinEnergyProd;
00179         G4double E2=kinEnergyProd*1.000001;
00180         G4double dE=(E2-E1);
00181         G4double sigma1,sigma2;
00182         theDirectEMModel =theBraggIonDirectEMModel;
00183         if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
00184         sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
00185         sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
00186         
00187         dSigmadEprod=(sigma1-sigma2)/dE;
00188         
00189         //G4double chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,E);
00190         
00191         
00192         
00193         if (dSigmadEprod>1.) {
00194                 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
00195                 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
00196                 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
00197                 
00198         }
00199         
00200          
00201         
00202          
00203          
00204          
00205          if (theDirectEMModel == theBetheBlochDirectEMModel ){
00206          //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
00207          //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
00208          //to test the rejection of a secondary
00209          //-------------------------
00210          
00211          //Source code taken from    G4BetheBlochModel::SampleSecondaries
00212          
00213                 G4double deltaKinEnergy = kinEnergyProd;
00214          
00215          //Part of the taken code
00216          //----------------------
00217          
00218          
00219          
00220          // projectile formfactor - suppresion of high energy
00221          // delta-electron production at high energy
00222          
00223          
00224                 G4double x = formfact*deltaKinEnergy;
00225                 if(x > 1.e-6) {
00226                         G4double totEnergy     = kinEnergyProj + mass;
00227                         G4double etot2         = totEnergy*totEnergy;
00228                         G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
00229                         G4double f;
00230                         G4double f1 = 0.0;
00231                         f = 1.0 - beta2*deltaKinEnergy/Tmax;
00232                         if( 0.5 == spin ) {
00233                                 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
00234                                 f += f1;
00235                         }
00236                         G4double x1 = 1.0 + x;
00237                         G4double gg  = 1.0/(x1*x1);
00238                         if( 0.5 == spin ) {
00239                                 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
00240                                 gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
00241                         }
00242                         if(gg > 1.0) {
00243                                 G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: gg= " << gg
00244                                 << G4endl;
00245                                 gg=1.;
00246                         }
00247                         //G4cout<<"gg"<<gg<<G4endl;
00248                         dSigmadEprod*=gg;               
00249                 }
00250         }       
00251          
00252   }
00253 
00254  return dSigmadEprod;   
00255 }

G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase ( G4double PrimAdjEnergy ) [virtual]

Reimplemented from G4VEmAdjointModel.

Definition at line 356 of file G4AdjointIonIonisationModel.cc.

References G4VEmAdjointModel::HighEnergyLimit.

Referenced by DiffCrossSectionPerAtomPrimToSecond().

00357 { return HighEnergyLimit;
00358 }

G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase ( G4double PrimAdjEnergy ) [virtual]

Reimplemented from G4VEmAdjointModel.

Definition at line 344 of file G4AdjointIonIonisationModel.cc.

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

G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForProdToProjCase ( G4double PrimAdjEnergy ) [virtual]

Reimplemented from G4VEmAdjointModel.

Definition at line 361 of file G4AdjointIonIonisationModel.cc.

Referenced by DiffCrossSectionPerAtomPrimToSecond().

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

G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase	(	G4double	PrimAdjEnergy,
		G4double	Tcut = `0`
	)			`[virtual]`

Reimplemented from G4VEmAdjointModel.

Definition at line 351 of file G4AdjointIonIonisationModel.cc.

00352 { return PrimAdjEnergy+Tcut;
00353 }

void G4AdjointIonIonisationModel::SampleSecondaries	(	const G4Track &	aTrack,
		G4bool	IsScatProjToProjCase,
		G4ParticleChange *	fParticleChange
	)			`[virtual]`

Implements G4VEmAdjointModel.

Definition at line 85 of file G4AdjointIonIonisationModel.cc.

References G4ParticleChange::AddSecondary(), CorrectPostStepWeight(), fStopAndKill, G4UniformRand, G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGMass(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetWeight(), G4VEmAdjointModel::HighEnergyLimit, G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, and G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef.

00088 { 
00089  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
00090   
00091  //Elastic inverse scattering 
00092  //---------------------------------------------------------
00093  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
00094  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
00095 
00096  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
00097         return;
00098  }
00099  
00100  //Sample secondary energy
00101  //-----------------------
00102  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
00103  CorrectPostStepWeight(fParticleChange, 
00104                               aTrack.GetWeight(), 
00105                               adjointPrimKinEnergy,
00106                               projectileKinEnergy,
00107                               IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
00108 
00109                  
00110  //Kinematic: 
00111  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest  and gives
00112  // him part of its  energy
00113  //----------------------------------------------------------------------------------------
00114  
00115  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00116  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
00117  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
00118  
00119  
00120  
00121  //Companion
00122  //-----------
00123  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00124  if (IsScatProjToProjCase) {
00125         companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
00126  } 
00127  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
00128  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;    
00129  
00130  
00131  //Projectile momentum
00132  //--------------------
00133  G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
00134  G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
00135  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
00136  G4double phi =G4UniformRand()*2.*3.1415926;
00137  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
00138  projectileMomentum.rotateUz(dir_parallel);
00139  
00140  
00141  
00142  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
00143         fParticleChange->ProposeTrackStatus(fStopAndKill);
00144         fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
00145         //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
00146  }
00147  else {
00148         fParticleChange->ProposeEnergy(projectileKinEnergy);
00149         fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
00150  }      
00151         
00152 
00153   
00154 
00155 }

void G4AdjointIonIonisationModel::SetIon	(	G4ParticleDefinition *	adj_ion,
		G4ParticleDefinition *	fwd_ion
	)

Definition at line 258 of file G4AdjointIonIonisationModel.cc.

References G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, and G4VEmAdjointModel::theDirectPrimaryPartDef.

00259 { theDirectPrimaryPartDef =fwd_ion;
00260   theAdjEquivOfDirectPrimPartDef =adj_ion;
00261  
00262   DefineProjectileProperty();
00263 }

void G4AdjointIonIonisationModel::SetUseOnlyBragg ( G4bool aBool ) [inline]

Definition at line 106 of file G4AdjointIonIonisationModel.hh.

00106 {use_only_bragg =aBool;}


Public Member Functions
	G4AdjointIonIonisationModel ()
virtual	~G4AdjointIonIonisationModel ()
virtual void	SampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
virtual G4double	DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
virtual void	CorrectPostStepWeight (G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, 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)
void	SetUseOnlyBragg (G4bool aBool)
void	SetIon (G4ParticleDefinition adj_ion, G4ParticleDefinition fwd_ion)

G4AdjointIonIonisationModel Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation