Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4AdjointComptonModel Class Reference

#include <G4AdjointComptonModel.hh>

Inheritance diagram for G4AdjointComptonModel:
G4VEmAdjointModel

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)
 

Additional Inherited Members

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

Detailed Description

Definition at line 54 of file G4AdjointComptonModel.hh.

Constructor & Destructor Documentation

G4AdjointComptonModel::G4AdjointComptonModel ( )

Definition at line 44 of file G4AdjointComptonModel.cc.

References G4AdjointElectron::AdjointElectron(), G4AdjointGamma::AdjointGamma(), G4Gamma::Gamma(), G4VEmAdjointModel::second_part_of_same_type, G4VEmAdjointModel::SetApplyCutInRange(), G4VEmAdjointModel::SetUseMatrix(), G4VEmAdjointModel::SetUseMatrixPerElement(), G4VEmAdjointModel::SetUseOnlyOneMatrixForAllElements(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef, G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.

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:298
G4AdjointElectron::AdjointElectron
static G4AdjointElectron * AdjointElectron()
Definition: G4AdjointElectron.cc:100
G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
Definition: G4VEmAdjointModel.hh:296
G4VEmAdjointModel::SetUseMatrixPerElement
void SetUseMatrixPerElement(G4bool aBool)
Definition: G4VEmAdjointModel.hh:207
G4VEmAdjointModel::theDirectEMModel
G4VEmModel * theDirectEMModel
Definition: G4VEmAdjointModel.hh:255
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:297
G4VEmAdjointModel::SetApplyCutInRange
void SetApplyCutInRange(G4bool aBool)
Definition: G4VEmAdjointModel.hh:210
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 376 of file G4AdjointComptonModel.cc.

References G4VEmAdjointModel::AdjointCrossSection(), G4VEmAdjointModel::currentMaterial, G4VEmAdjointModel::DefineCurrentMaterial(), python.hepunit::electron_mass_c2, G4Material::GetElectronDensity(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(), GetSecondAdjEnergyMaxForScatProjToProjCase(), GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(), G4VEmAdjointModel::lastCS, python.hepunit::twopi_mc2_rcl2, and G4VEmAdjointModel::UseMatrix.

Referenced by GetAdjointCrossSection().

379 {
380  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
381  DefineCurrentMaterial(aCouple);
382 
383 
384  float Cross=0.;
385  float Emax_proj =0.;
386  float Emin_proj =0.;
387  if (!IsScatProjToProjCase ){
388  Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
389  Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
390  if (Emax_proj>Emin_proj ){
391  Cross= 0.1*std::log((Emax_proj-float (primEnergy))*Emin_proj/Emax_proj/(Emin_proj-primEnergy))
392  *(1.+2.*std::log(float(1.+electron_mass_c2/primEnergy)));
393  }
394  }
395  else {
396  Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
397  Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.);
398  if (Emax_proj>Emin_proj) {
399  Cross = 0.1*std::log(Emax_proj/Emin_proj);
400  //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment
401  }
402 
403 
404  }
405 
406  Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
407  lastCS=Cross;
408  return double(Cross);
409 }
G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:671
G4VEmAdjointModel::currentMaterial
G4Material * currentMaterial
Definition: G4VEmAdjointModel.hh:308
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:215
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointComptonModel.cc:360
G4VEmAdjointModel::DefineCurrentMaterial
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
Definition: G4VEmAdjointModel.cc:683
G4VEmAdjointModel::AdjointCrossSection
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
Definition: G4VEmAdjointModel.cc:58
G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Definition: G4VEmAdjointModel.cc:664
G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointComptonModel.cc:368
G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim ( G4double  kinEnergyProj,
G4double  kinEnergyScatProj,
G4double  Z,
G4double  A = 0. 
)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 293 of file G4AdjointComptonModel.cc.

References G4VEmModel::ComputeCrossSectionPerAtom(), python.hepunit::electron_mass_c2, G4Gamma::Gamma(), G4VEmAdjointModel::theDirectEMModel, and test::v.

Referenced by DiffCrossSectionPerAtomPrimToSecond(), and RapidSampleSecondaries().

298 { //Based on Klein Nishina formula
299  // In the forward case (see G4KleinNishinaModel) the cross section is parametrised while
300  // the secondaries are sampled from the
301  // Klein Nishida differential cross section
302  // The used diffrential cross section here is therefore the cross section multiplied by the normalised
303  //differential Klein Nishida cross section
304 
305 
306  //Klein Nishida Cross Section
307  //-----------------------------
308  G4double epsilon = gamEnergy0 / electron_mass_c2 ;
309  G4double one_plus_two_epsi =1.+2.*epsilon;
310  G4double gamEnergy1_max = gamEnergy0;
311  G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi;
312  if (gamEnergy1 >gamEnergy1_max || gamEnergy1<gamEnergy1_min) {
313  /*G4cout<<"the differential CS is null"<<G4endl;
314  G4cout<<gamEnergy0<<G4endl;
315  G4cout<<gamEnergy1<<G4endl;
316  G4cout<<gamEnergy1_min<<G4endl;*/
317  return 0.;
318  }
319 
320 
321  G4double epsi2 = epsilon *epsilon ;
322  G4double one_plus_two_epsi_2=one_plus_two_epsi*one_plus_two_epsi;
323 
324 
325  G4double CS=std::log(one_plus_two_epsi)*(1.- 2.*(1.+epsilon)/epsi2);
326  CS+=4./epsilon +0.5*(1.-1./one_plus_two_epsi_2);
327  CS/=epsilon;
328  //Note that the pi*re2*Z factor is neglected because it is suppresed when computing dCS_dE1/CS;
329  // in the differential cross section
330 
331 
332  //Klein Nishida Differential Cross Section
333  //-----------------------------------------
334  G4double epsilon1 = gamEnergy1 / electron_mass_c2 ;
335  G4double v= epsilon1/epsilon;
336  G4double term1 =1.+ 1./epsilon -1/epsilon1;
337  G4double dCS_dE1= 1./v +v + term1*term1 -1.;
338  dCS_dE1 *=1./epsilon/gamEnergy0;
339 
340 
341  //Normalised to the CS used in G4
342  //-------------------------------
343 
344  G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
345  gamEnergy0,
346  Z, 0., 0.,0.);
347 
348  dCS_dE1 *= G4direct_CS/CS;
349 /* G4cout<<"the differential CS is not null"<<G4endl;
350  G4cout<<gamEnergy0<<G4endl;
351  G4cout<<gamEnergy1<<G4endl;*/
352 
353  return dCS_dE1;
354 
355 
356 }
G4VEmAdjointModel::theDirectEMModel
G4VEmModel * theDirectEMModel
Definition: G4VEmAdjointModel.hh:255
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
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:300
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToSecond ( G4double  kinEnergyProj,
G4double  kinEnergyProd,
G4double  Z,
G4double  A = 0. 
)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 278 of file G4AdjointComptonModel.cc.

References DiffCrossSectionPerAtomPrimToScatPrim().

283 {
284  G4double gamEnergy1 = gamEnergy0 - kinEnergyElec;
285  G4double dSigmadEprod=0.;
286  if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A);
287  return dSigmadEprod;
288 }
G4double
double G4double
Definition: G4Types.hh:76
G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
Definition: G4AdjointComptonModel.cc:293
G4double G4AdjointComptonModel::GetAdjointCrossSection ( const G4MaterialCutsCouple *  aCouple,
G4double  primEnergy,
G4bool  IsScatProjToProjCase 
)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 412 of file G4AdjointComptonModel.cc.

References AdjointCrossSection().

415 { return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
416  //return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
417 }
G4AdjointComptonModel::AdjointCrossSection
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
Definition: G4AdjointComptonModel.cc:376
G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 360 of file G4AdjointComptonModel.cc.

References python.hepunit::electron_mass_c2, G4VEmAdjointModel::HighEnergyLimit, and G4INCL::Math::min().

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

361 { G4double inv_e_max = 1./PrimAdjEnergy - 2./electron_mass_c2;
362  G4double e_max = HighEnergyLimit;
363  if (inv_e_max > 0. ) e_max=std::min(1./inv_e_max,HighEnergyLimit);
364  return e_max;
365 }
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 368 of file G4AdjointComptonModel.cc.

References python.hepunit::electron_mass_c2.

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

369 { G4double half_e=PrimAdjEnergy/2.;
370  G4double term=std::sqrt(half_e*(electron_mass_c2+half_e));
371  G4double emin=half_e+term;
372  return emin;
373 }
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4double
double G4double
Definition: G4Types.hh:76
void G4AdjointComptonModel::RapidSampleSecondaries ( const G4Track &  aTrack,
G4bool  IsScatProjToProjCase,
G4ParticleChange *  fParticleChange 
)

Definition at line 156 of file G4AdjointComptonModel.cc.

References G4ParticleChange::AddSecondary(), G4VEmAdjointModel::currentCouple, G4VEmAdjointModel::currentMaterial, G4VEmAdjointModel::currentTcutForDirectSecond, G4VEmAdjointModel::DefineCurrentMaterial(), DiffCrossSectionPerAtomPrimToScatPrim(), python.hepunit::electron_mass_c2, fStopAndKill, G4UniformRand, G4AdjointCSManager::GetAdjointCSManager(), G4Track::GetDynamicParticle(), G4Material::GetElectronDensity(), G4DynamicParticle::GetKineticEnergy(), G4VEmProcess::GetLambda(), G4Track::GetMaterialCutsCouple(), G4DynamicParticle::GetMomentumDirection(), G4AdjointCSManager::GetPostStepWeightCorrection(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(), GetSecondAdjEnergyMaxForScatProjToProjCase(), GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetWeight(), G4VEmAdjointModel::HighEnergyLimit, G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeParentWeight(), G4VParticleChange::ProposeTrackStatus(), CLHEP::Hep3Vector::rotateUz(), G4VParticleChange::SetParentWeightByProcess(), G4VParticleChange::SetSecondaryWeightByProcess(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, python.hepunit::twopi_mc2_rcl2, and CLHEP::Hep3Vector::unit().

Referenced by SampleSecondaries().

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 
202  //Weight correction
203  //-----------------------
204  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
205  G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
206 
207  //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the
208  //one consistent with the direct model
209 
210 
211  G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.);
212  if (diffCS >0) diffCS /=G4direct_CS; // here we have the normalised diffCS
213  diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
214  //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1);
215  //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl;
216 
217  w_corr*=diffCS/diffCSUsed;
218 
219  G4double new_weight = aTrack.GetWeight()*w_corr;
220  fParticleChange->SetParentWeightByProcess(false);
221  fParticleChange->SetSecondaryWeightByProcess(false);
222  fParticleChange->ProposeParentWeight(new_weight);
223 
224 
225 
226  //Cos th
227  //-------
228 
229  G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
230  if (!IsScatProjToProjCase) {
231  G4double p_elec=theAdjointPrimary->GetTotalMomentum();
232  cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
233  }
234  G4double sin_th = 0.;
235  if (std::abs(cos_th)>1){
236  //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
237  if (cos_th>0) {
238  cos_th=1.;
239  }
240  else cos_th=-1.;
241  sin_th=0.;
242  }
243  else sin_th = std::sqrt(1.-cos_th*cos_th);
244 
245 
246 
247 
248  //gamma0 momentum
249  //--------------------
250 
251 
252  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
253  G4double phi =G4UniformRand()*2.*3.1415926;
254  G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
255  gammaMomentum1.rotateUz(dir_parallel);
256 
257 
258 
259 
260  if (!IsScatProjToProjCase){ //kill the primary and add a secondary
261  fParticleChange->ProposeTrackStatus(fStopAndKill);
262  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
263  //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
264  }
265  else {
266  fParticleChange->ProposeEnergy(gammaE1);
267  fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
268  }
269 
270 
271 
272 }
G4VEmProcess::GetLambda
G4double GetLambda(G4double &kinEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VEmProcess.hh:501
G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:671
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4AdjointCSManager::GetPostStepWeightCorrection
G4double GetPostStepWeightCorrection()
Definition: G4AdjointCSManager.cc:525
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4VParticleChange::ProposeParentWeight
void ProposeParentWeight(G4double finalWeight)
G4Track::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4VEmAdjointModel::currentMaterial
G4Material * currentMaterial
Definition: G4VEmAdjointModel.hh:308
G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
Definition: G4VEmAdjointModel.hh:296
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:87
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:215
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointComptonModel.cc:360
G4VEmAdjointModel::currentCouple
G4MaterialCutsCouple * currentCouple
Definition: G4VEmAdjointModel.hh:309
G4VEmAdjointModel::DefineCurrentMaterial
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
Definition: G4VEmAdjointModel.cc:683
G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Definition: G4VEmAdjointModel.cc:664
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
G4VEmAdjointModel::currentTcutForDirectSecond
G4double currentTcutForDirectSecond
Definition: G4VEmAdjointModel.hh:313
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:368
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:293
void G4AdjointComptonModel::SampleSecondaries ( const G4Track &  aTrack,
G4bool  IsScatProjToProjCase,
G4ParticleChange *  fParticleChange 
)
virtual

Implements G4VEmAdjointModel.

Definition at line 64 of file G4AdjointComptonModel.cc.

References G4ParticleChange::AddSecondary(), G4VEmAdjointModel::CorrectPostStepWeight(), python.hepunit::electron_mass_c2, fStopAndKill, G4UniformRand, G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetWeight(), G4VEmAdjointModel::HighEnergyLimit, G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), RapidSampleSecondaries(), CLHEP::Hep3Vector::rotateUz(), G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, CLHEP::Hep3Vector::unit(), and G4VEmAdjointModel::UseMatrix.

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:454
G4AdjointComptonModel::RapidSampleSecondaries
void RapidSampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
Definition: G4AdjointComptonModel.cc:156
G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
Definition: G4VEmAdjointModel.hh:296
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
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:624
void G4AdjointComptonModel::SetDirectProcess ( G4VEmProcess *  aProcess)
inline

Definition at line 95 of file G4AdjointComptonModel.hh.

Referenced by G4AdjointPhysicsList::ConstructEM().

95 {theDirectEMProcess = aProcess;};
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:10 for Geant4.10 by   doxygen 1.8.7