G4AdjointhIonisationModel Class Reference

#include <G4AdjointhIonisationModel.hh>

Inheritance diagram for G4AdjointhIonisationModel:

G4VEmAdjointModel

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)

Detailed Description

Definition at line 72 of file G4AdjointhIonisationModel.hh.


Constructor & Destructor Documentation

G4AdjointhIonisationModel::G4AdjointhIonisationModel ( G4ParticleDefinition projectileDefinition  ) 

Definition at line 46 of file G4AdjointhIonisationModel.cc.

References G4AdjointElectron::AdjointElectron(), G4AdjointProton::AdjointProton(), G4VEmAdjointModel::ApplyCutInRange, G4VEmAdjointModel::CS_biasing_factor, G4Proton::Proton(), G4VEmAdjointModel::second_part_of_same_type, G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef, G4VEmAdjointModel::theDirectEMModel, G4VEmAdjointModel::theDirectPrimaryPartDef, G4VEmAdjointModel::UseMatrix, G4VEmAdjointModel::UseMatrixPerElement, and G4VEmAdjointModel::UseOnlyOneMatrixForAllElements.

00046                                                                                               :
00047   G4VEmAdjointModel("Adjoint_hIonisation")
00048 { 
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   
00059   //The direct EM Modfel is taken has BetheBloch it is only used for the computation 
00060   // of the differential cross section.
00061   //The Bragg model could be used  as an alternative as  it offers the same differential cross section
00062   
00063   theDirectEMModel = new G4BetheBlochModel(projectileDefinition);
00064   theBraggDirectEMModel = new G4BraggModel(projectileDefinition); 
00065   theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
00066   
00067   theDirectPrimaryPartDef = projectileDefinition;
00068   theAdjEquivOfDirectPrimPartDef = 0;
00069   if (projectileDefinition == G4Proton::Proton()) {
00070         theAdjEquivOfDirectPrimPartDef = G4AdjointProton::AdjointProton();
00071         
00072   }
00073   
00074   DefineProjectileProperty();
00075 }

G4AdjointhIonisationModel::~G4AdjointhIonisationModel (  )  [virtual]

Definition at line 78 of file G4AdjointhIonisationModel.cc.

00079 {;}


Member Function Documentation

G4double G4AdjointhIonisationModel::AdjointCrossSection ( const G4MaterialCutsCouple aCouple,
G4double  primEnergy,
G4bool  IsScatProjToProjCase 
) [virtual]

Reimplemented from G4VEmAdjointModel.

Definition at line 438 of file G4AdjointhIonisationModel.cc.

References G4VEmAdjointModel::AdjointCrossSection(), G4VEmAdjointModel::currentMaterial, G4VEmAdjointModel::currentTcutForDirectSecond, G4VEmAdjointModel::DefineCurrentMaterial(), G4Material::GetElectronDensity(), GetSecondAdjEnergyMaxForProdToProjCase(), GetSecondAdjEnergyMaxForScatProjToProjCase(), GetSecondAdjEnergyMinForProdToProjCase(), GetSecondAdjEnergyMinForScatProjToProjCase(), G4VEmAdjointModel::lastCS, and G4VEmAdjointModel::UseMatrix.

00441 { 
00442   if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
00443   DefineCurrentMaterial(aCouple);
00444         
00445   
00446   G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
00447   
00448   if (!IsScatProjToProjCase ){
00449         G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
00450         G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
00451         if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
00452                 Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
00453         } 
00454         else Cross=0.;
00455                 
00456         
00457 
00458         
00459         
00460         
00461   }
00462   else {
00463         G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
00464         G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
00465         G4double diff1=Emin_proj-primEnergy;
00466         G4double diff2=Emax_proj-primEnergy;
00467         G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
00468         //G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
00469         G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy;
00470         Cross*=(t1+t2);
00471 
00472         
00473   }
00474   lastCS =Cross;
00475   return Cross; 
00476 }

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

Reimplemented from G4VEmAdjointModel.

Definition at line 306 of file G4AdjointhIonisationModel.cc.

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

Referenced by RapidSampleSecondaries().

00311 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
00312   
00313 
00314  
00315  G4double dSigmadEprod=0;
00316  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00317  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00318  
00319  
00320  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile 
00321         G4double Tmax=kinEnergyProj;
00322         
00323         G4double E1=kinEnergyProd;
00324         G4double E2=kinEnergyProd*1.000001;
00325         G4double dE=(E2-E1);
00326         G4double sigma1,sigma2;
00327         if (kinEnergyProj >2.*MeV){
00328             sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
00329             sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
00330         }
00331         else {
00332             sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
00333             sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
00334         }
00335         
00336         
00337         dSigmadEprod=(sigma1-sigma2)/dE;
00338         if (dSigmadEprod>1.) {
00339                 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
00340                 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
00341                 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
00342                 
00343         }
00344         
00345          
00346         
00347          //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
00348          //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
00349          //to test the rejection of a secondary
00350          //-------------------------
00351          
00352          //Source code taken from    G4BetheBlochModel::SampleSecondaries
00353          
00354          G4double deltaKinEnergy = kinEnergyProd;
00355          
00356          //Part of the taken code
00357          //----------------------
00358          
00359          
00360          
00361          // projectile formfactor - suppresion of high energy
00362          // delta-electron production at high energy
00363          G4double x = formfact*deltaKinEnergy;
00364           if(x > 1.e-6) {
00365                 
00366                 
00367                 G4double totEnergy     = kinEnergyProj + mass;
00368                 G4double etot2         = totEnergy*totEnergy;
00369                 G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
00370                 G4double f;
00371                 G4double f1 = 0.0;
00372                 f = 1.0 - beta2*deltaKinEnergy/Tmax;
00373                 if( 0.5 == spin ) {
00374                         f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
00375                         f += f1;
00376                 }
00377                 G4double x1 = 1.0 + x;
00378                 G4double gg  = 1.0/(x1*x1);
00379                 if( 0.5 == spin ) {
00380                         G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
00381                         gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
00382                 }
00383                 if(gg > 1.0) {
00384                         G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
00385                                 << G4endl;
00386                         gg=1.;
00387                 }
00388                 //G4cout<<"gg"<<gg<<G4endl;
00389                 dSigmadEprod*=gg;               
00390          }
00391          
00392   }
00393 
00394  return dSigmadEprod;   
00395 }

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

Reimplemented from G4VEmAdjointModel.

Definition at line 491 of file G4AdjointhIonisationModel.cc.

References G4VEmAdjointModel::HighEnergyLimit.

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

00492 { return HighEnergyLimit;
00493 }

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

Reimplemented from G4VEmAdjointModel.

Definition at line 479 of file G4AdjointhIonisationModel.cc.

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

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

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

Reimplemented from G4VEmAdjointModel.

Definition at line 496 of file G4AdjointhIonisationModel.cc.

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

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

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

Reimplemented from G4VEmAdjointModel.

Definition at line 486 of file G4AdjointhIonisationModel.cc.

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

00487 { return PrimAdjEnergy+Tcut;
00488 }

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

Definition at line 160 of file G4AdjointhIonisationModel.cc.

References G4ParticleChange::AddSecondary(), G4VEmAdjointModel::currentMaterial, G4VEmAdjointModel::currentTcutForDirectSecond, G4VEmAdjointModel::DefineCurrentMaterial(), DiffCrossSectionPerAtomPrimToSecond(), fStopAndKill, G4UniformRand, G4AdjointCSManager::GetAdjointCSManager(), G4Track::GetDynamicParticle(), G4Material::GetElectronDensity(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterialCutsCouple(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGMass(), G4AdjointCSManager::GetPostStepWeightCorrection(), GetSecondAdjEnergyMaxForProdToProjCase(), GetSecondAdjEnergyMaxForScatProjToProjCase(), GetSecondAdjEnergyMinForProdToProjCase(), GetSecondAdjEnergyMinForScatProjToProjCase(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetWeight(), G4VEmAdjointModel::HighEnergyLimit, G4VEmAdjointModel::lastCS, G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeParentWeight(), G4VParticleChange::ProposeTrackStatus(), G4VParticleChange::SetParentWeightByProcess(), G4VParticleChange::SetSecondaryWeightByProcess(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, and G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef.

Referenced by SampleSecondaries().

00163 { 
00164 
00165  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
00166  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
00167  
00168  
00169  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
00170  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
00171  
00172  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
00173         return;
00174  }
00175   
00176  G4double projectileKinEnergy =0.;
00177  G4double eEnergy=0.;
00178  G4double newCS=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
00179  if (!IsScatProjToProjCase){//1/E^2 distribution
00180         
00181         eEnergy=adjointPrimKinEnergy;
00182         G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
00183         G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
00184         if (Emin>=Emax) return;
00185         G4double a=1./Emax;
00186         G4double b=1./Emin;
00187         newCS=newCS*(b-a)/eEnergy;
00188         
00189         projectileKinEnergy =1./(b- (b-a)*G4UniformRand()); 
00190         
00191         
00192  }
00193  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
00194         G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
00195         if (Emin>=Emax) return;
00196         G4double diff1=Emin-adjointPrimKinEnergy;
00197         G4double diff2=Emax-adjointPrimKinEnergy;
00198         
00199         G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
00200         G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
00201         /*G4double f31=diff1/Emin;
00202         G4double f32=diff2/Emax/f31;*/
00203         G4double t3=2.*std::log(Emax/Emin);
00204         G4double sum_t=t1+t2+t3;
00205         newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
00206         G4double t=G4UniformRand()*sum_t;
00207         if (t <=t1 ){
00208                 G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ;
00209                 projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
00210                 
00211         }
00212         else if (t <=t2 )  {
00213                 G4double q= G4UniformRand()*t2/adjointPrimKinEnergy;
00214                 projectileKinEnergy =1./(1./Emin-q);
00215         }
00216         else {  
00217                 projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
00218                 
00219         }
00220         eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
00221         
00222         
00223  }
00224 
00225  
00226 
00227  G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy; 
00228   
00229   
00230                                                         
00231  //Weight correction
00232  //-----------------------
00233  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
00234  G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
00235  
00236  //G4cout<<w_corr<<G4endl;
00237  w_corr*=newCS/lastCS;
00238  //G4cout<<w_corr<<G4endl;
00239  //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
00240  //Here we consider the true  diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS
00241  
00242  G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1);
00243  w_corr*=diffCS/diffCS_perAtom_Used;
00244  //G4cout<<w_corr<<G4endl;
00245            
00246  G4double new_weight = aTrack.GetWeight()*w_corr;
00247  fParticleChange->SetParentWeightByProcess(false);
00248  fParticleChange->SetSecondaryWeightByProcess(false);
00249  fParticleChange->ProposeParentWeight(new_weight);
00250  
00251  
00252  
00253  
00254  //Kinematic: 
00255  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest  and gives
00256  // him part of its  energy
00257  //----------------------------------------------------------------------------------------
00258  
00259  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00260  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
00261  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
00262  
00263  
00264  
00265  //Companion
00266  //-----------
00267  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00268  if (IsScatProjToProjCase) {
00269         companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
00270  } 
00271  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
00272  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;    
00273  
00274  
00275  //Projectile momentum
00276  //--------------------
00277  G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
00278  G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
00279  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
00280  G4double phi =G4UniformRand()*2.*3.1415926;
00281  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
00282  projectileMomentum.rotateUz(dir_parallel);
00283  
00284  
00285  
00286  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
00287         fParticleChange->ProposeTrackStatus(fStopAndKill);
00288         fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
00289         //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
00290  }
00291  else {
00292         fParticleChange->ProposeEnergy(projectileKinEnergy);
00293         fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
00294  }      
00295         
00296 
00297  
00298  
00299 
00300 
00301 
00302 } 

void G4AdjointhIonisationModel::SampleSecondaries ( const G4Track aTrack,
G4bool  IsScatProjToProjCase,
G4ParticleChange fParticleChange 
) [virtual]

Implements G4VEmAdjointModel.

Definition at line 84 of file G4AdjointhIonisationModel.cc.

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

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


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:24 2013 for Geant4 by  doxygen 1.4.7