Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4AdjointComptonModel Class Reference

#include <G4AdjointComptonModel.hh>

Inheritance diagram for G4AdjointComptonModel:
Inheritance graph
[legend]
Collaboration diagram for G4AdjointComptonModel:
Collaboration graph
[legend]

Public Member Functions

 G4AdjointComptonModel ()
 
 ~G4AdjointComptonModel ()
 
virtual void SampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
 
void RapidSampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
 
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim (G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
 
virtual G4double DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
 
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase (G4double PrimAdjEnergy)
 
virtual G4double GetSecondAdjEnergyMinForProdToProjCase (G4double PrimAdjEnergy)
 
virtual G4double AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
 
virtual G4double GetAdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
 
void SetDirectProcess (G4VEmProcess *aProcess)
 
- Public Member Functions inherited from G4VEmAdjointModel
 G4VEmAdjointModel (const G4String &nam)
 
virtual ~G4VEmAdjointModel ()
 
virtual G4double DiffCrossSectionPerVolumePrimToSecond (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
 
virtual G4double DiffCrossSectionPerVolumePrimToScatPrim (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
 
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase (G4double PrimAdjEnergy, G4double Tcut=0)
 
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase (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)
 
void SetCorrectWeightForPostStepInModel (G4bool aBool)
 
void SetAdditionalWeightCorrectionFactorForPostStepOutsideModel (G4double factor)
 

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
 
G4bool correct_weight_for_post_step_in_model
 
G4double additional_weight_correction_factor_for_post_step_outside_model
 

Detailed Description

Definition at line 54 of file G4AdjointComptonModel.hh.

Constructor & Destructor Documentation

G4AdjointComptonModel::G4AdjointComptonModel ( )

Definition at line 44 of file G4AdjointComptonModel.cc.

44  :
45  G4VEmAdjointModel("AdjointCompton")
46 
47 { SetApplyCutInRange(false);
48  SetUseMatrix(false);
49  SetUseMatrixPerElement(true);
50  SetUseOnlyOneMatrixForAllElements(true);
51  theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
52  theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
53  theDirectPrimaryPartDef=G4Gamma::Gamma();
54  second_part_of_same_type=false;
55  theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel");
56  G4direct_CS = 0.;
57 }
G4AdjointGamma::AdjointGamma
static G4AdjointGamma * AdjointGamma()
Definition: G4AdjointGamma.cc:84
G4VEmAdjointModel::theDirectPrimaryPartDef
G4ParticleDefinition * theDirectPrimaryPartDef
Definition: G4VEmAdjointModel.hh:297
G4AdjointElectron::AdjointElectron
static G4AdjointElectron * AdjointElectron()
Definition: G4AdjointElectron.cc:100
G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
Definition: G4VEmAdjointModel.hh:295
G4VEmAdjointModel::SetUseMatrixPerElement
void SetUseMatrixPerElement(G4bool aBool)
Definition: G4VEmAdjointModel.hh:207
G4VEmAdjointModel::theDirectEMModel
G4VEmModel * theDirectEMModel
Definition: G4VEmAdjointModel.hh:258
G4VEmAdjointModel::G4VEmAdjointModel
G4VEmAdjointModel(const G4String &nam)
Definition: G4VEmAdjointModel.cc:41
G4VEmAdjointModel::SetUseMatrix
void SetUseMatrix(G4bool aBool)
Definition: G4VEmAdjointModel.hh:205
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4VEmAdjointModel::SetUseOnlyOneMatrixForAllElements
void SetUseOnlyOneMatrixForAllElements(G4bool aBool)
Definition: G4VEmAdjointModel.hh:208
G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
Definition: G4VEmAdjointModel.hh:296
G4VEmAdjointModel::SetApplyCutInRange
void SetApplyCutInRange(G4bool aBool)
Definition: G4VEmAdjointModel.hh:210

Here is the call graph for this function:

G4AdjointComptonModel::~G4AdjointComptonModel ( )

Definition at line 60 of file G4AdjointComptonModel.cc.

61 {;}

Member Function Documentation

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

Reimplemented from G4VEmAdjointModel.

Definition at line 378 of file G4AdjointComptonModel.cc.

381 {
382  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
383  DefineCurrentMaterial(aCouple);
384 
385 
386  float Cross=0.;
387  float Emax_proj =0.;
388  float Emin_proj =0.;
389  if (!IsScatProjToProjCase ){
390  Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
391  Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
392  if (Emax_proj>Emin_proj ){
393  Cross= 0.1*std::log((Emax_proj-float (primEnergy))*Emin_proj/Emax_proj/(Emin_proj-primEnergy))
394  *(1.+2.*std::log(float(1.+electron_mass_c2/primEnergy)));
395  }
396  }
397  else {
398  Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
399  Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.);
400  if (Emax_proj>Emin_proj) {
401  Cross = 0.1*std::log(Emax_proj/Emin_proj);
402  //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment
403  }
404 
405 
406  }
407 
408  Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
409  lastCS=Cross;
410  return double(Cross);
411 }
G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:681
CLHEP::twopi_mc2_rcl2
static constexpr double twopi_mc2_rcl2
Definition: PhysicalConstants.h:105
G4VEmAdjointModel::currentMaterial
G4Material * currentMaterial
Definition: G4VEmAdjointModel.hh:306
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointComptonModel.cc:362
G4VEmAdjointModel::DefineCurrentMaterial
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
Definition: G4VEmAdjointModel.cc:693
G4VEmAdjointModel::AdjointCrossSection
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
Definition: G4VEmAdjointModel.cc:59
G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Definition: G4VEmAdjointModel.cc:674
G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointComptonModel.cc:370

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim ( G4double  kinEnergyProj,
G4double  kinEnergyScatProj,
G4double  Z,
G4double  A = 0. 
)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 295 of file G4AdjointComptonModel.cc.

300 { //Based on Klein Nishina formula
301  // In the forward case (see G4KleinNishinaCompton) the cross section is parametrised
302  // but the secondaries are sampled from the
303  // Klein Nishida differential cross section
304  // The used diffrential cross section here is therefore the cross section multiplied by the normalised
305  //differential Klein Nishida cross section
306 
307 
308  //Klein Nishida Cross Section
309  //-----------------------------
310  G4double epsilon = gamEnergy0 / electron_mass_c2 ;
311  G4double one_plus_two_epsi =1.+2.*epsilon;
312  G4double gamEnergy1_max = gamEnergy0;
313  G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi;
314  if (gamEnergy1 >gamEnergy1_max || gamEnergy1<gamEnergy1_min) {
315  /*G4cout<<"the differential CS is null"<<G4endl;
316  G4cout<<gamEnergy0<<G4endl;
317  G4cout<<gamEnergy1<<G4endl;
318  G4cout<<gamEnergy1_min<<G4endl;*/
319  return 0.;
320  }
321 
322 
323  G4double epsi2 = epsilon *epsilon ;
324  G4double one_plus_two_epsi_2=one_plus_two_epsi*one_plus_two_epsi;
325 
326 
327  G4double CS=std::log(one_plus_two_epsi)*(1.- 2.*(1.+epsilon)/epsi2);
328  CS+=4./epsilon +0.5*(1.-1./one_plus_two_epsi_2);
329  CS/=epsilon;
330  //Note that the pi*re2*Z factor is neglected because it is suppresed when computing dCS_dE1/CS;
331  // in the differential cross section
332 
333 
334  //Klein Nishida Differential Cross Section
335  //-----------------------------------------
336  G4double epsilon1 = gamEnergy1 / electron_mass_c2 ;
337  G4double v= epsilon1/epsilon;
338  G4double term1 =1.+ 1./epsilon -1/epsilon1;
339  G4double dCS_dE1= 1./v +v + term1*term1 -1.;
340  dCS_dE1 *=1./epsilon/gamEnergy0;
341 
342 
343  //Normalised to the CS used in G4
344  //-------------------------------
345 
346  G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
347  gamEnergy0,
348  Z, 0., 0.,0.);
349 
350  dCS_dE1 *= G4direct_CS/CS;
351 /* G4cout<<"the differential CS is not null"<<G4endl;
352  G4cout<<gamEnergy0<<G4endl;
353  G4cout<<gamEnergy1<<G4endl;*/
354 
355  return dCS_dE1;
356 
357 
358 }
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4VEmAdjointModel::theDirectEMModel
G4VEmModel * theDirectEMModel
Definition: G4VEmAdjointModel.hh:258
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4VEmModel::ComputeCrossSectionPerAtom
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:321
G4double
double G4double
Definition: G4Types.hh:76
epsilon
double epsilon(double density, double temperature)
Definition: G4DNAElectronHoleRecombination.cc:96
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

Here is the caller graph for this function:

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

Reimplemented from G4VEmAdjointModel.

Definition at line 280 of file G4AdjointComptonModel.cc.

285 {
286  G4double gamEnergy1 = gamEnergy0 - kinEnergyElec;
287  G4double dSigmadEprod=0.;
288  if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A);
289  return dSigmadEprod;
290 }
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:60
G4double
double G4double
Definition: G4Types.hh:76
G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
Definition: G4AdjointComptonModel.cc:295
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

G4double G4AdjointComptonModel::GetAdjointCrossSection ( const G4MaterialCutsCouple *  aCouple,
G4double  primEnergy,
G4bool  IsScatProjToProjCase 
)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 414 of file G4AdjointComptonModel.cc.

417 { return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
418  //return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
419 }
G4AdjointComptonModel::AdjointCrossSection
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
Definition: G4AdjointComptonModel.cc:378

Here is the call graph for this function:

G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 362 of file G4AdjointComptonModel.cc.

363 { G4double inv_e_max = 1./PrimAdjEnergy - 2./electron_mass_c2;
364  G4double e_max = HighEnergyLimit;
365  if (inv_e_max > 0. ) e_max=std::min(1./inv_e_max,HighEnergyLimit);
366  return e_max;
367 }
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 370 of file G4AdjointComptonModel.cc.

371 { G4double half_e=PrimAdjEnergy/2.;
372  G4double term=std::sqrt(half_e*(electron_mass_c2+half_e));
373  G4double emin=half_e+term;
374  return emin;
375 }
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4double
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

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

Definition at line 156 of file G4AdjointComptonModel.cc.

159 {
160 
161  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
162  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
163 
164 
165  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
166 
167 
168  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
169  return;
170  }
171 
172 
173 
174  G4double diffCSUsed=0.1*currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
175  G4double gammaE1=0.;
176  G4double gammaE2=0.;
177  if (!IsScatProjToProjCase){
178 
179  G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
180  G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
181  if (Emin>=Emax) return;
182  G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
183  G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
184  gammaE1=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));;
185  gammaE2=gammaE1-adjointPrimKinEnergy;
186  diffCSUsed= diffCSUsed*(1.+2.*std::log(1.+electron_mass_c2/adjointPrimKinEnergy))*adjointPrimKinEnergy/gammaE1/gammaE2;
187 
188 
189  }
190  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
191  G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
192  if (Emin>=Emax) return;
193  gammaE2 =adjointPrimKinEnergy;
194  gammaE1=Emin*std::pow(Emax/Emin,G4UniformRand());
195  diffCSUsed= diffCSUsed/gammaE1;
196  }
197 
198 
199 
200 
201  //Weight correction
202  //-----------------------
203  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
204  G4double w_corr=additional_weight_correction_factor_for_post_step_outside_model;
205  if (correct_weight_for_post_step_in_model) {
206  w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
207  }
208  //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the
209  //one consistent with the direct model
210 
211 
212  G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.);
213  if (diffCS >0) diffCS /=G4direct_CS; // here we have the normalised diffCS
214  //An we remultiply by the lambda of the forward process
215  diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
216  //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1);
217  //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl;
218 
219  w_corr*=diffCS/diffCSUsed;
220 
221  G4double new_weight = aTrack.GetWeight()*w_corr;
222  fParticleChange->SetParentWeightByProcess(false);
223  fParticleChange->SetSecondaryWeightByProcess(false);
224  fParticleChange->ProposeParentWeight(new_weight);
225 
226 
227 
228  //Cos th
229  //-------
230 
231  G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
232  if (!IsScatProjToProjCase) {
233  G4double p_elec=theAdjointPrimary->GetTotalMomentum();
234  cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
235  }
236  G4double sin_th = 0.;
237  if (std::abs(cos_th)>1){
238  //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
239  if (cos_th>0) {
240  cos_th=1.;
241  }
242  else cos_th=-1.;
243  sin_th=0.;
244  }
245  else sin_th = std::sqrt(1.-cos_th*cos_th);
246 
247 
248 
249 
250  //gamma0 momentum
251  //--------------------
252 
253 
254  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
255  G4double phi =G4UniformRand()*2.*3.1415926;
256  G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
257  gammaMomentum1.rotateUz(dir_parallel);
258 
259 
260 
261 
262  if (!IsScatProjToProjCase){ //kill the primary and add a secondary
263  fParticleChange->ProposeTrackStatus(fStopAndKill);
264  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
265  //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
266  }
267  else {
268  fParticleChange->ProposeEnergy(gammaE1);
269  fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
270  }
271 
272 
273 
274 }
G4VEmProcess::GetLambda
G4double GetLambda(G4double &kinEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VEmProcess.hh:519
G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:681
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4AdjointCSManager::GetPostStepWeightCorrection
G4double GetPostStepWeightCorrection()
Definition: G4AdjointCSManager.cc:527
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
CLHEP::twopi_mc2_rcl2
static constexpr double twopi_mc2_rcl2
Definition: PhysicalConstants.h:105
nf_Legendre_from_ptwXY_callback_s::f2
double f2
Definition: nf_Legendre.cc:15
G4VParticleChange::ProposeParentWeight
void ProposeParentWeight(G4double finalWeight)
G4Track::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4VEmAdjointModel::currentMaterial
G4Material * currentMaterial
Definition: G4VEmAdjointModel.hh:306
G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
Definition: G4VEmAdjointModel.hh:295
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
G4VParticleChange::SetSecondaryWeightByProcess
void SetSecondaryWeightByProcess(G4bool)
G4VParticleChange::SetParentWeightByProcess
void SetParentWeightByProcess(G4bool)
Definition: G4VParticleChange.cc:522
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:217
nf_Legendre_from_ptwXY_callback_s::f1
double f1
Definition: nf_Legendre.cc:15
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointComptonModel.cc:362
G4VEmAdjointModel::currentCouple
G4MaterialCutsCouple * currentCouple
Definition: G4VEmAdjointModel.hh:307
G4VEmAdjointModel::DefineCurrentMaterial
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
Definition: G4VEmAdjointModel.cc:693
G4VEmAdjointModel::correct_weight_for_post_step_in_model
G4bool correct_weight_for_post_step_in_model
Definition: G4VEmAdjointModel.hh:344
G4VEmAdjointModel::additional_weight_correction_factor_for_post_step_outside_model
G4double additional_weight_correction_factor_for_post_step_outside_model
Definition: G4VEmAdjointModel.hh:345
G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Definition: G4VEmAdjointModel.cc:674
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4ParticleChange::ProposeEnergy
void ProposeEnergy(G4double finalEnergy)
fStopAndKill
Definition: G4TrackStatus.hh:56
G4ParticleChange::AddSecondary
void AddSecondary(G4Track *aSecondary)
Definition: G4ParticleChange.cc:218
G4Track::GetWeight
G4double GetWeight() const
Emin
static const G4double Emin
Definition: G4PhotoNuclearCrossSection.cc:50
Emax
static const G4double Emax
Definition: G4PhotoNuclearCrossSection.cc:51
G4VEmAdjointModel::currentTcutForDirectSecond
G4double currentTcutForDirectSecond
Definition: G4VEmAdjointModel.hh:311
G4ParticleChange::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double
double G4double
Definition: G4Types.hh:76
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointComptonModel.cc:370
G4AdjointCSManager::GetAdjointCSManager
static G4AdjointCSManager * GetAdjointCSManager()
Definition: G4AdjointCSManager.cc:61
G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
Definition: G4AdjointComptonModel.cc:295

Here is the call graph for this function:

Here is the caller graph for this function:

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

Implements G4VEmAdjointModel.

Definition at line 64 of file G4AdjointComptonModel.cc.

67 {
68  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
69 
70  //A recall of the compton scattering law is
71  //Egamma2=Egamma1/(1+(Egamma1/E0_electron)(1.-cos_th))
72  //Therefore Egamma2_max= Egamma2(cos_th=1) = Egamma1
73  //Therefore Egamma2_min= Egamma2(cos_th=-1) = Egamma1/(1+2.(Egamma1/E0_electron))
74 
75 
76  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
77  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
78  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
79  return;
80  }
81 
82 
83  //Sample secondary energy
84  //-----------------------
85  G4double gammaE1;
86  gammaE1 = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
87  IsScatProjToProjCase);
88 
89 
90  //gammaE2
91  //-----------
92 
93  G4double gammaE2 = adjointPrimKinEnergy;
94  if (!IsScatProjToProjCase) gammaE2 = gammaE1 - adjointPrimKinEnergy;
95 
96 
97 
98 
99 
100 
101  //Cos th
102  //-------
103 // G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<G4endl;
104  G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
105  if (!IsScatProjToProjCase) {
106  G4double p_elec=theAdjointPrimary->GetTotalMomentum();
107  cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
108  }
109  G4double sin_th = 0.;
110  if (std::abs(cos_th)>1){
111  //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
112  if (cos_th>0) {
113  cos_th=1.;
114  }
115  else cos_th=-1.;
116  sin_th=0.;
117  }
118  else sin_th = std::sqrt(1.-cos_th*cos_th);
119 
120 
121 
122 
123  //gamma0 momentum
124  //--------------------
125 
126 
127  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
128  G4double phi =G4UniformRand()*2.*3.1415926;
129  G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
130  gammaMomentum1.rotateUz(dir_parallel);
131 // G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<G4endl;
132 
133 
134  //It is important to correct the weight of particles before adding the secondary
135  //------------------------------------------------------------------------------
136  CorrectPostStepWeight(fParticleChange,
137  aTrack.GetWeight(),
138  adjointPrimKinEnergy,
139  gammaE1,
140  IsScatProjToProjCase);
141 
142  if (!IsScatProjToProjCase){ //kill the primary and add a secondary
143  fParticleChange->ProposeTrackStatus(fStopAndKill);
144  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
145  //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
146  }
147  else {
148  fParticleChange->ProposeEnergy(gammaE1);
149  fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
150  }
151 
152 
153 }
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
Definition: G4VEmAdjointModel.cc:459
G4AdjointComptonModel::RapidSampleSecondaries
void RapidSampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
Definition: G4AdjointComptonModel.cc:156
G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
Definition: G4VEmAdjointModel.hh:295
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4ParticleChange::ProposeEnergy
void ProposeEnergy(G4double finalEnergy)
fStopAndKill
Definition: G4TrackStatus.hh:56
G4ParticleChange::AddSecondary
void AddSecondary(G4Track *aSecondary)
Definition: G4ParticleChange.cc:218
G4Track::GetWeight
G4double GetWeight() const
G4ParticleChange::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double
double G4double
Definition: G4Types.hh:76
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VEmAdjointModel::CorrectPostStepWeight
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
Definition: G4VEmAdjointModel.cc:634

Here is the call graph for this function:

void G4AdjointComptonModel::SetDirectProcess ( G4VEmProcess *  aProcess)
inline

Definition at line 95 of file G4AdjointComptonModel.hh.

95 {theDirectEMProcess = aProcess;};
The documentation for this class was generated from the following files: