Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4AdjointComptonModel Class Reference

#include <G4AdjointComptonModel.hh>

+ Inheritance 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 void SampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)=0
 
virtual G4double AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
 
virtual G4double GetAdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
 
virtual G4double DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
 
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim (G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
 
virtual G4double DiffCrossSectionPerVolumePrimToSecond (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
 
virtual G4double DiffCrossSectionPerVolumePrimToScatPrim (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
 
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase (G4double PrimAdjEnergy)
 
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase (G4double PrimAdjEnergy, G4double Tcut=0)
 
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase (G4double PrimAdjEnergy)
 
virtual G4double GetSecondAdjEnergyMinForProdToProjCase (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 53 of file G4AdjointComptonModel.hh.

Constructor & Destructor Documentation

◆ G4AdjointComptonModel()

G4AdjointComptonModel::G4AdjointComptonModel ( )

Definition at line 43 of file G4AdjointComptonModel.cc.

43 :
44 G4VEmAdjointModel("AdjointCompton")
45
46{ SetApplyCutInRange(false);
47 SetUseMatrix(false);
48 SetUseMatrixPerElement(true);
49 SetUseOnlyOneMatrixForAllElements(true);
50 theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
51 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
52 theDirectPrimaryPartDef=G4Gamma::Gamma();
53 second_part_of_same_type=false;
54 theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel");
55 G4direct_CS = 0.;
56}
G4AdjointElectron::AdjointElectron
static G4AdjointElectron * AdjointElectron()
Definition: G4AdjointElectron.cc:99
G4AdjointGamma::AdjointGamma
static G4AdjointGamma * AdjointGamma()
Definition: G4AdjointGamma.cc:83
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:85
G4VEmAdjointModel::SetUseMatrixPerElement
void SetUseMatrixPerElement(G4bool aBool)
Definition: G4VEmAdjointModel.hh:206
G4VEmAdjointModel::theDirectEMModel
G4VEmModel * theDirectEMModel
Definition: G4VEmAdjointModel.hh:257
G4VEmAdjointModel::SetUseMatrix
void SetUseMatrix(G4bool aBool)
Definition: G4VEmAdjointModel.hh:204
G4VEmAdjointModel::theDirectPrimaryPartDef
G4ParticleDefinition * theDirectPrimaryPartDef
Definition: G4VEmAdjointModel.hh:296
G4VEmAdjointModel::SetUseOnlyOneMatrixForAllElements
void SetUseOnlyOneMatrixForAllElements(G4bool aBool)
Definition: G4VEmAdjointModel.hh:207
G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
Definition: G4VEmAdjointModel.hh:295
G4VEmAdjointModel::SetApplyCutInRange
void SetApplyCutInRange(G4bool aBool)
Definition: G4VEmAdjointModel.hh:209
G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
Definition: G4VEmAdjointModel.hh:294

◆ ~G4AdjointComptonModel()

G4AdjointComptonModel::~G4AdjointComptonModel ( )

Definition at line 59 of file G4AdjointComptonModel.cc.

60{;}

Member Function Documentation

◆ AdjointCrossSection()

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

Reimplemented from G4VEmAdjointModel.

Definition at line 377 of file G4AdjointComptonModel.cc.

380{
381 if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
382 DefineCurrentMaterial(aCouple);
383
384
385 float Cross=0.;
386 float Emax_proj =0.;
387 float Emin_proj =0.;
388 if (!IsScatProjToProjCase ){
389 Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
390 Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
391 if (Emax_proj>Emin_proj ){
392 Cross= 0.1*std::log((Emax_proj-float (primEnergy))*Emin_proj/Emax_proj/(Emin_proj-primEnergy))
393 *(1.+2.*std::log(float(1.+electron_mass_c2/primEnergy)));
394 }
395 }
396 else {
397 Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
398 Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.);
399 if (Emax_proj>Emin_proj) {
400 Cross = 0.1*std::log(Emax_proj/Emin_proj);
401 //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment
402 }
403
404
405 }
406
407 Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
408 lastCS=Cross;
409 return double(Cross);
410}
G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointComptonModel.cc:361
G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointComptonModel.cc:369
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:215
G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:680
G4VEmAdjointModel::DefineCurrentMaterial
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
Definition: G4VEmAdjointModel.cc:692
G4VEmAdjointModel::currentMaterial
G4Material * currentMaterial
Definition: G4VEmAdjointModel.hh:305
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:673

Referenced by GetAdjointCrossSection().

◆ DiffCrossSectionPerAtomPrimToScatPrim()

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

Reimplemented from G4VEmAdjointModel.

Definition at line 294 of file G4AdjointComptonModel.cc.

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

Referenced by DiffCrossSectionPerAtomPrimToSecond(), and RapidSampleSecondaries().

◆ DiffCrossSectionPerAtomPrimToSecond()

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

Reimplemented from G4VEmAdjointModel.

Definition at line 279 of file G4AdjointComptonModel.cc.

284{
285 G4double gamEnergy1 = gamEnergy0 - kinEnergyElec;
286 G4double dSigmadEprod=0.;
287 if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A);
288 return dSigmadEprod;
289}
G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
Definition: G4AdjointComptonModel.cc:294

◆ GetAdjointCrossSection()

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

Reimplemented from G4VEmAdjointModel.

Definition at line 413 of file G4AdjointComptonModel.cc.

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

◆ GetSecondAdjEnergyMaxForScatProjToProjCase()

G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 361 of file G4AdjointComptonModel.cc.

362{ G4double inv_e_max = 1./PrimAdjEnergy - 2./electron_mass_c2;
363 G4double e_max = HighEnergyLimit;
364 if (inv_e_max > 0. ) e_max=std::min(1./inv_e_max,HighEnergyLimit);
365 return e_max;
366}

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

◆ GetSecondAdjEnergyMinForProdToProjCase()

G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 369 of file G4AdjointComptonModel.cc.

370{ G4double half_e=PrimAdjEnergy/2.;
371 G4double term=std::sqrt(half_e*(electron_mass_c2+half_e));
372 G4double emin=half_e+term;
373 return emin;
374}

Referenced by AdjointCrossSection(), and RapidSampleSecondaries().

◆ RapidSampleSecondaries()

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

Definition at line 155 of file G4AdjointComptonModel.cc.

158{
159
160 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
161 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
162
163
164 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
165
166
167 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
168 return;
169 }
170
171
172
173 G4double diffCSUsed=0.1*currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
174 G4double gammaE1=0.;
175 G4double gammaE2=0.;
176 if (!IsScatProjToProjCase){
177
178 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
179 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
180 if (Emin>=Emax) return;
181 G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
182 G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
183 gammaE1=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));;
184 gammaE2=gammaE1-adjointPrimKinEnergy;
185 diffCSUsed= diffCSUsed*(1.+2.*std::log(1.+electron_mass_c2/adjointPrimKinEnergy))*adjointPrimKinEnergy/gammaE1/gammaE2;
186
187
188 }
189 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
190 G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
191 if (Emin>=Emax) return;
192 gammaE2 =adjointPrimKinEnergy;
193 gammaE1=Emin*std::pow(Emax/Emin,G4UniformRand());
194 diffCSUsed= diffCSUsed/gammaE1;
195 }
196
197
198
199
200 //Weight correction
201 //-----------------------
202 //First w_corr is set to the ratio between adjoint total CS and fwd total CS
203 G4double w_corr=additional_weight_correction_factor_for_post_step_outside_model;
204 if (correct_weight_for_post_step_in_model) {
205 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 //An we remultiply by the lambda of the forward process
214 diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
215 //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1);
216 //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl;
217
218 w_corr*=diffCS/diffCSUsed;
219
220 G4double new_weight = aTrack.GetWeight()*w_corr;
221 fParticleChange->SetParentWeightByProcess(false);
222 fParticleChange->SetSecondaryWeightByProcess(false);
223 fParticleChange->ProposeParentWeight(new_weight);
224
225
226
227 //Cos th
228 //-------
229
230 G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
231 if (!IsScatProjToProjCase) {
232 G4double p_elec=theAdjointPrimary->GetTotalMomentum();
233 cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
234 }
235 G4double sin_th = 0.;
236 if (std::abs(cos_th)>1){
237 //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
238 if (cos_th>0) {
239 cos_th=1.;
240 }
241 else cos_th=-1.;
242 sin_th=0.;
243 }
244 else sin_th = std::sqrt(1.-cos_th*cos_th);
245
246
247
248
249 //gamma0 momentum
250 //--------------------
251
252
253 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
254 G4double phi =G4UniformRand()*2.*3.1415926;
255 G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
256 gammaMomentum1.rotateUz(dir_parallel);
257
258
259
260
261 if (!IsScatProjToProjCase){ //kill the primary and add a secondary
262 fParticleChange->ProposeTrackStatus(fStopAndKill);
263 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
264 //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
265 }
266 else {
267 fParticleChange->ProposeEnergy(gammaE1);
268 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
269 }
270
271
272
273}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
fStopAndKill
@ fStopAndKill
Definition: G4TrackStatus.hh:46
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:33
G4AdjointCSManager::GetPostStepWeightCorrection
G4double GetPostStepWeightCorrection()
Definition: G4AdjointCSManager.cc:526
G4AdjointCSManager::GetAdjointCSManager
static G4AdjointCSManager * GetAdjointCSManager()
Definition: G4AdjointCSManager.cc:60
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
G4ParticleChange::AddSecondary
void AddSecondary(G4Track *aSecondary)
Definition: G4ParticleChange.cc:187
G4ParticleChange::ProposeEnergy
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4Track::GetWeight
G4double GetWeight() const
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4Track::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4VEmAdjointModel::additional_weight_correction_factor_for_post_step_outside_model
G4double additional_weight_correction_factor_for_post_step_outside_model
Definition: G4VEmAdjointModel.hh:344
G4VEmAdjointModel::correct_weight_for_post_step_in_model
G4bool correct_weight_for_post_step_in_model
Definition: G4VEmAdjointModel.hh:343
G4VEmAdjointModel::currentTcutForDirectSecond
G4double currentTcutForDirectSecond
Definition: G4VEmAdjointModel.hh:310
G4VEmAdjointModel::currentCouple
G4MaterialCutsCouple * currentCouple
Definition: G4VEmAdjointModel.hh:306
G4VEmProcess::GetLambda
G4double GetLambda(G4double kinEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VEmProcess.cc:1259
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VParticleChange::SetSecondaryWeightByProcess
void SetSecondaryWeightByProcess(G4bool)
G4VParticleChange::SetParentWeightByProcess
void SetParentWeightByProcess(G4bool)
Definition: G4VParticleChange.cc:515
G4VParticleChange::ProposeParentWeight
void ProposeParentWeight(G4double finalWeight)

Referenced by SampleSecondaries().

◆ SampleSecondaries()

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

Implements G4VEmAdjointModel.

Definition at line 63 of file G4AdjointComptonModel.cc.

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

◆ SetDirectProcess()

void G4AdjointComptonModel::SetDirectProcess ( G4VEmProcess *  aProcess)
inline

Definition at line 94 of file G4AdjointComptonModel.hh.

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