Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointComptonModel Class Reference

#include <G4AdjointComptonModel.hh>

Inheritance diagram for G4AdjointComptonModel:
Collaboration diagram for G4AdjointComptonModel:

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)
 
G4ParticleDefinitionGetAdjointEquivalentOfDirectPrimaryParticleDefinition ()
 
G4ParticleDefinitionGetAdjointEquivalentOfDirectSecondaryParticleDefinition ()
 
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
G4VEmModeltheDirectEMModel
 
G4VParticleChangepParticleChange
 
const G4String name
 
G4int ASelectedNucleus
 
G4int ZSelectedNucleus
 
G4MaterialSelectedMaterial
 
G4double kinEnergyProdForIntegration
 
G4double kinEnergyScatProjForIntegration
 
G4double kinEnergyProjForIntegration
 
std::vector< G4AdjointCSMatrix * > * pOnCSMatrixForProdToProjBackwardScattering
 
std::vector< G4AdjointCSMatrix * > * pOnCSMatrixForScatProjToProjBackwardScattering
 
std::vector< G4doubleCS_Vs_ElementForScatProjToProjCase
 
std::vector< G4doubleCS_Vs_ElementForProdToProjCase
 
G4double lastCS
 
G4double lastAdjointCSForScatProjToProjCase
 
G4double lastAdjointCSForProdToProjCase
 
G4ParticleDefinitiontheAdjEquivOfDirectPrimPartDef
 
G4ParticleDefinitiontheAdjEquivOfDirectSecondPartDef
 
G4ParticleDefinitiontheDirectPrimaryPartDef
 
G4bool second_part_of_same_type
 
G4double preStepEnergy
 
G4MaterialcurrentMaterial
 
G4MaterialCutsCouplecurrentCouple
 
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);
55  theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel");
56  G4direct_CS = 0.;
57 }
static G4AdjointGamma * AdjointGamma()
G4ParticleDefinition * theDirectPrimaryPartDef
static G4AdjointElectron * AdjointElectron()
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
void SetUseMatrixPerElement(G4bool aBool)
G4VEmModel * theDirectEMModel
G4VEmAdjointModel(const G4String &nam)
void SetUseMatrix(G4bool aBool)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
void SetUseOnlyOneMatrixForAllElements(G4bool aBool)
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
void SetApplyCutInRange(G4bool aBool)

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 
409  lastCS=Cross;
410  return double(Cross);
411 }
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
G4Material * currentMaterial
G4double GetElectronDensity() const
Definition: G4Material.hh:217
float electron_mass_c2
Definition: hepunit.py:274
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)

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 
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 }
G4VEmModel * theDirectEMModel
float electron_mass_c2
Definition: hepunit.py:274
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:321
tuple v
Definition: test.py:18
double G4double
Definition: G4Types.hh:76
double epsilon(double density, double temperature)

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 }
double A(double temperature)
double G4double
Definition: G4Types.hh:76
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)

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 }
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)

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 }
float electron_mass_c2
Definition: hepunit.py:274
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
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 }
float electron_mass_c2
Definition: hepunit.py:274
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();
163 
164 
165  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
166 
167 
168  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
169  return;
170  }
171 
172 
173 
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);
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
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 }
G4double GetLambda(G4double &kinEnergy, const G4MaterialCutsCouple *couple)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetPostStepWeightCorrection()
const G4DynamicParticle * GetDynamicParticle() const
void ProposeParentWeight(G4double finalWeight)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
G4double GetTotalMomentum() const
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
#define G4UniformRand()
Definition: Randomize.hh:97
G4double GetElectronDensity() const
Definition: G4Material.hh:217
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
float electron_mass_c2
Definition: hepunit.py:274
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
G4MaterialCutsCouple * currentCouple
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4bool correct_weight_for_post_step_in_model
G4double additional_weight_correction_factor_for_post_step_outside_model
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Hep3Vector unit() const
void ProposeEnergy(G4double finalEnergy)
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
static const G4double Emin
static const G4double Emax
G4double currentTcutForDirectSecond
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
static G4AdjointCSManager * GetAdjointCSManager()
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)

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 }
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
void RapidSampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
G4double GetTotalMomentum() const
#define G4UniformRand()
Definition: Randomize.hh:97
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
float electron_mass_c2
Definition: hepunit.py:274
Hep3Vector unit() const
void ProposeEnergy(G4double finalEnergy)
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)

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;};

Here is the caller graph for this function:


The documentation for this class was generated from the following files: