G4AdjointBremsstrahlungModel Class Reference

#include <G4AdjointBremsstrahlungModel.hh>

Inheritance diagram for G4AdjointBremsstrahlungModel:

G4VEmAdjointModel

Public Member Functions

 G4AdjointBremsstrahlungModel (G4VEmModel *aModel)
 G4AdjointBremsstrahlungModel ()
 ~G4AdjointBremsstrahlungModel ()
virtual void SampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
void RapidSampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
virtual G4double DiffCrossSectionPerVolumePrimToSecond (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
G4double DiffCrossSectionPerVolumePrimToSecondApproximated1 (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
G4double DiffCrossSectionPerVolumePrimToSecondApproximated2 (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
virtual G4double AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
virtual G4double GetAdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)

Detailed Description

Definition at line 63 of file G4AdjointBremsstrahlungModel.hh.


Constructor & Destructor Documentation

G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel ( G4VEmModel aModel  ) 

Definition at line 46 of file G4AdjointBremsstrahlungModel.cc.

References G4EmModelManager::AddEmModel(), G4AdjointElectron::AdjointElectron(), G4AdjointGamma::AdjointGamma(), G4Electron::Electron(), G4VEmAdjointModel::second_part_of_same_type, G4VEmAdjointModel::SetApplyCutInRange(), G4VEmAdjointModel::SetUseMatrix(), G4VEmAdjointModel::SetUseMatrixPerElement(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef, G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.

00046                                                                             :
00047  G4VEmAdjointModel("AdjointeBremModel")
00048 { 
00049   SetUseMatrix(false);
00050   SetUseMatrixPerElement(false);
00051   
00052   theDirectStdBremModel = aModel;
00053   theDirectEMModel=theDirectStdBremModel;
00054   theEmModelManagerForFwdModels = new G4EmModelManager();
00055   isDirectModelInitialised = false;
00056   G4VEmFluctuationModel* f=0;
00057   G4Region* r=0;
00058   theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r);
00059 
00060   SetApplyCutInRange(true);
00061   highKinEnergy= 100.*TeV;
00062   lowKinEnergy = 1.0*keV;
00063 
00064   lastCZ =0.;
00065 
00066   
00067   theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
00068   theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma();
00069   theDirectPrimaryPartDef=G4Electron::Electron();
00070   second_part_of_same_type=false;
00071   
00072   /*UsePenelopeModel=false;
00073   if (UsePenelopeModel) {
00074         G4PenelopeBremsstrahlungModel* thePenelopeModel = new G4PenelopeBremsstrahlungModel(G4Electron::Electron(),"PenelopeBrem");
00075         theEmModelManagerForFwdModels = new G4EmModelManager();
00076         isPenelopeModelInitialised = false;
00077         G4VEmFluctuationModel* f=0;
00078         G4Region* r=0;
00079         theDirectEMModel=thePenelopeModel;
00080         theEmModelManagerForFwdModels->AddEmModel(1, thePenelopeModel, f, r);
00081   }
00082   */    
00083   
00084 
00085   
00086 }

G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel (  ) 

Definition at line 89 of file G4AdjointBremsstrahlungModel.cc.

References G4EmModelManager::AddEmModel(), G4AdjointElectron::AdjointElectron(), G4AdjointGamma::AdjointGamma(), G4Electron::Electron(), G4VEmAdjointModel::second_part_of_same_type, G4VEmAdjointModel::SetApplyCutInRange(), G4VEmAdjointModel::SetUseMatrix(), G4VEmAdjointModel::SetUseMatrixPerElement(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef, G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.

00089                                                           :
00090  G4VEmAdjointModel("AdjointeBremModel")
00091 {
00092   SetUseMatrix(false);
00093   SetUseMatrixPerElement(false);
00094 
00095   theDirectStdBremModel = new G4SeltzerBergerModel();
00096   theDirectEMModel=theDirectStdBremModel;
00097   theEmModelManagerForFwdModels = new G4EmModelManager();
00098   isDirectModelInitialised = false;
00099   G4VEmFluctuationModel* f=0;
00100   G4Region* r=0;
00101   theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r);
00102  // theDirectPenelopeBremModel =0;
00103   SetApplyCutInRange(true);
00104   highKinEnergy= 1.*GeV;
00105   lowKinEnergy = 1.0*keV;
00106   lastCZ =0.;
00107   theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
00108   theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma();
00109   theDirectPrimaryPartDef=G4Electron::Electron();
00110   second_part_of_same_type=false;
00111 
00112 }

G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel (  ) 

Definition at line 115 of file G4AdjointBremsstrahlungModel.cc.

00116 {if (theDirectStdBremModel) delete theDirectStdBremModel;
00117  if (theEmModelManagerForFwdModels) delete theEmModelManagerForFwdModels;
00118 }


Member Function Documentation

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

Reimplemented from G4VEmAdjointModel.

Definition at line 399 of file G4AdjointBremsstrahlungModel.cc.

References G4VEmAdjointModel::AdjointCrossSection(), G4VEmModel::CrossSectionPerVolume(), G4VEmAdjointModel::currentTcutForDirectSecond, G4VEmAdjointModel::DefineCurrentMaterial(), G4Electron::Electron(), G4Gamma::Gamma(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(), G4EmModelManager::Initialise(), G4VEmAdjointModel::theDirectEMModel, G4VEmAdjointModel::theDirectPrimaryPartDef, and G4VEmAdjointModel::UseMatrix.

Referenced by GetAdjointCrossSection().

00402 { if (!isDirectModelInitialised) {
00403         theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
00404         isDirectModelInitialised =true;
00405   }
00406   if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
00407   DefineCurrentMaterial(aCouple);
00408   G4double Cross=0.;
00409   lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above
00410   
00411   if (!IsScatProjToProjCase ){
00412         G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
00413         G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
00414         if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= lastCZ*std::log(Emax_proj/Emin_proj);
00415   }
00416   else {
00417         G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
00418         G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
00419         if (Emax_proj>Emin_proj) Cross= lastCZ*std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy));
00420         
00421   }
00422   return Cross; 
00423 }                                            

G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond ( const G4Material aMaterial,
G4double  kinEnergyProj,
G4double  kinEnergyProd 
) [virtual]

Reimplemented from G4VEmAdjointModel.

Definition at line 310 of file G4AdjointBremsstrahlungModel.cc.

References DiffCrossSectionPerVolumePrimToSecondApproximated2(), G4Electron::Electron(), G4Gamma::Gamma(), and G4EmModelManager::Initialise().

Referenced by RapidSampleSecondaries().

00314 {if (!isDirectModelInitialised) {
00315         theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
00316         isDirectModelInitialised =true;
00317  }
00318 
00319  return  DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial,
00320                                                                  kinEnergyProj, 
00321                                                                  kinEnergyProd);
00322  /*return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
00323                                                                  kinEnergyProj, 
00324                                                                  kinEnergyProd);*/                                                                                                                               
00325 }                                     

G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1 ( const G4Material aMaterial,
G4double  kinEnergyProj,
G4double  kinEnergyProd 
)

Definition at line 329 of file G4AdjointBremsstrahlungModel.cc.

References G4VEmModel::CrossSectionPerVolume(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.

00334 {
00335  G4double dCrossEprod=0.;
00336  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00337  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00338  
00339  
00340  //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution
00341  //This is what is applied in the discrete standard model before the  rejection test  that make a correction
00342  //The application of the same rejection function is not possible here.
00343  //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the 
00344  // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and 
00345  // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one. 
00346  // In the future we plan to use the brem secondary spectra from the G4Penelope implementation 
00347  
00348  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
00349         G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV);
00350         dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV);
00351  }
00352  return dCrossEprod;
00353   
00354 }

G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2 ( const G4Material aMaterial,
G4double  kinEnergyProj,
G4double  kinEnergyProd 
)

Definition at line 358 of file G4AdjointBremsstrahlungModel.cc.

References C1, G4VEmModel::ComputeCrossSectionPerAtom(), G4Material::GetAtomicNumDensityVector(), G4Material::GetElementVector(), G4Material::GetNumberOfElements(), G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.

Referenced by DiffCrossSectionPerVolumePrimToSecond().

00363 {
00364  //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor 
00365   //used in the direct model
00366  
00367  G4double dCrossEprod=0.;
00368  
00369  const G4ElementVector* theElementVector = material->GetElementVector();
00370  const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
00371  G4double dum=0.;
00372  G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001;
00373  G4double dE=E2-E1;
00374  for (size_t i=0; i<material->GetNumberOfElements(); i++) { 
00375         G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1);
00376         G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2);
00377         dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE;
00378    
00379  }
00380  
00381  //Now the Migdal correction
00382 /*
00383  G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
00384  G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
00385                                              *(material->GetElectronDensity());
00386                                              
00387  
00388  G4double MigdalFactor = 1./(1.+kp2/(kinEnergyProd*kinEnergyProd)); // its seems that the factor used in the CS compuation i the direct
00389                                                                     //model is different than the one used in the secondary sampling by a
00390                                                                     //factor (1.+kp2) To be checked!
00391  
00392  dCrossEprod*=MigdalFactor;
00393  */
00394  return dCrossEprod;
00395   
00396 }

G4double G4AdjointBremsstrahlungModel::GetAdjointCrossSection ( const G4MaterialCutsCouple aCouple,
G4double  primEnergy,
G4bool  IsScatProjToProjCase 
) [virtual]

Reimplemented from G4VEmAdjointModel.

Definition at line 425 of file G4AdjointBremsstrahlungModel.cc.

References AdjointCrossSection(), G4VEmModel::CrossSectionPerVolume(), G4VEmAdjointModel::GetAdjointCrossSection(), G4MaterialCutsCouple::GetMaterial(), G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.

00428 { 
00429   return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
00430   lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above
00431   return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
00432         
00433 }

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

Definition at line 200 of file G4AdjointBremsstrahlungModel.cc.

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

Referenced by SampleSecondaries().

00203 { 
00204 
00205  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
00206  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
00207  
00208  
00209  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
00210  G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
00211  
00212  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
00213         return;
00214  }
00215   
00216  G4double projectileKinEnergy =0.;
00217  G4double gammaEnergy=0.;
00218  G4double diffCSUsed=0.; 
00219  if (!IsScatProjToProjCase){
00220         gammaEnergy=adjointPrimKinEnergy;
00221         G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
00222         G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
00223         if (Emin>=Emax) return;
00224         projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
00225         diffCSUsed=lastCZ/projectileKinEnergy;
00226         
00227  }
00228  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
00229         G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
00230         if (Emin>=Emax) return;
00231         G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
00232         G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
00233         //G4cout<<"f1 and f2 "<<f1<<'\t'<<f2<<G4endl;
00234         projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));
00235         gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy;
00236         diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy;
00237         
00238  }
00239   
00240   
00241   
00242                                                         
00243  //Weight correction
00244  //-----------------------
00245  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
00246  G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
00247 
00248  //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
00249  //Here we consider the true  diffCS as the one obtained by the numericla differentiation over Tcut of the direct CS, corrected by the Migdal term.
00250  //Basically any other differential CS   diffCS could be used here (example Penelope). 
00251  
00252  G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy);
00253  w_corr*=diffCS/diffCSUsed;
00254            
00255  G4double new_weight = aTrack.GetWeight()*w_corr;
00256  fParticleChange->SetParentWeightByProcess(false);
00257  fParticleChange->SetSecondaryWeightByProcess(false);
00258  fParticleChange->ProposeParentWeight(new_weight);
00259  
00260  //Kinematic
00261  //---------
00262  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00263  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
00264  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
00265  G4double projectileP = std::sqrt(projectileP2);
00266  
00267  
00268  //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
00269  //------------------------------------------------
00270   G4double u;
00271   const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
00272 
00273   if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
00274      else                          u = - std::log(G4UniformRand()*G4UniformRand())/a2;
00275 
00276   G4double theta = u*electron_mass_c2/projectileTotalEnergy;
00277 
00278   G4double sint = std::sin(theta);
00279   G4double cost = std::cos(theta);
00280 
00281   G4double phi = twopi * G4UniformRand() ;
00282   
00283   G4ThreeVector projectileMomentum;
00284   projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
00285   if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
00286         G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
00287         G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
00288         G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
00289         G4double sint1 =  std::sqrt(1.-cost1*cost1);
00290         projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
00291   
00292   }
00293   
00294   projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
00295  
00296  
00297  
00298   if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
00299         fParticleChange->ProposeTrackStatus(fStopAndKill);
00300         fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
00301   }
00302   else {
00303         fParticleChange->ProposeEnergy(projectileKinEnergy);
00304         fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
00305         
00306   }     
00307 } 

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

Implements G4VEmAdjointModel.

Definition at line 122 of file G4AdjointBremsstrahlungModel.cc.

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

00125 {
00126  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 
00127 
00128  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
00129  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
00130  
00131  
00132  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
00133  G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
00134  
00135  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
00136         return;
00137  }
00138   
00139   G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
00140                                                         IsScatProjToProjCase);
00141  //Weight correction
00142  //-----------------------                                         
00143  CorrectPostStepWeight(fParticleChange, 
00144                        aTrack.GetWeight(), 
00145                        adjointPrimKinEnergy,
00146                        projectileKinEnergy,
00147                        IsScatProjToProjCase);   
00148  
00149  
00150  //Kinematic
00151  //---------
00152  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00153  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
00154  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
00155  G4double projectileP = std::sqrt(projectileP2);
00156  
00157  
00158  //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
00159  //------------------------------------------------
00160   G4double u;
00161   const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
00162 
00163   if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
00164      else                          u = - std::log(G4UniformRand()*G4UniformRand())/a2;
00165 
00166   G4double theta = u*electron_mass_c2/projectileTotalEnergy;
00167 
00168   G4double sint = std::sin(theta);
00169   G4double cost = std::cos(theta);
00170 
00171   G4double phi = twopi * G4UniformRand() ;
00172   
00173   G4ThreeVector projectileMomentum;
00174   projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
00175   if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
00176         G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
00177         G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
00178         G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
00179         G4double sint1 =  std::sqrt(1.-cost1*cost1);
00180         projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
00181   
00182   }
00183   
00184   projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
00185  
00186  
00187  
00188   if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
00189         fParticleChange->ProposeTrackStatus(fStopAndKill);
00190         fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
00191   }
00192   else {
00193         fParticleChange->ProposeEnergy(projectileKinEnergy);
00194         fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
00195         
00196   }     
00197 } 


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