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

#include <G4AdjointBremsstrahlungModel.hh>

Inheritance diagram for G4AdjointBremsstrahlungModel:
Collaboration diagram for G4AdjointBremsstrahlungModel:

Public Member Functions

 G4AdjointBremsstrahlungModel (G4VEmModel *aModel)
 
 G4AdjointBremsstrahlungModel ()
 
 ~G4AdjointBremsstrahlungModel ()
 
virtual void SampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
 
void RapidSampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
 
virtual G4double DiffCrossSectionPerVolumePrimToSecond (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
 
G4double DiffCrossSectionPerVolumePrimToSecondApproximated1 (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
 
G4double DiffCrossSectionPerVolumePrimToSecondApproximated2 (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
 
virtual G4double AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
 
virtual G4double GetAdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
 
- Public Member Functions inherited from G4VEmAdjointModel
 G4VEmAdjointModel (const G4String &nam)
 
virtual ~G4VEmAdjointModel ()
 
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 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)
 
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 62 of file G4AdjointBremsstrahlungModel.hh.

Constructor & Destructor Documentation

G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel ( G4VEmModel aModel)

Definition at line 46 of file G4AdjointBremsstrahlungModel.cc.

46  :
47  G4VEmAdjointModel("AdjointeBremModel")
48 {
49  SetUseMatrix(false);
51 
52  theDirectStdBremModel = aModel;
53  theDirectEMModel=theDirectStdBremModel;
54  theEmModelManagerForFwdModels = new G4EmModelManager();
55  isDirectModelInitialised = false;
57  G4Region* r=0;
58  theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r);
59 
60  SetApplyCutInRange(true);
61  highKinEnergy= 1.*GeV;
62  lowKinEnergy = 1.0*keV;
63 
64  lastCZ =0.;
65 
66 
71 
72 
74 
75 
76 }
static G4AdjointGamma * AdjointGamma()
G4ParticleDefinition * theDirectPrimaryPartDef
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
static G4AdjointElectron * AdjointElectron()
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
void SetUseMatrixPerElement(G4bool aBool)
G4VEmModel * theDirectEMModel
G4VEmAdjointModel(const G4String &nam)
void SetUseMatrix(G4bool aBool)
static constexpr double GeV
Definition: G4SIunits.hh:217
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
void SetApplyCutInRange(G4bool aBool)
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel ( )

Definition at line 79 of file G4AdjointBremsstrahlungModel.cc.

79  :
80  G4VEmAdjointModel("AdjointeBremModel")
81 {
82  SetUseMatrix(false);
84 
85  theDirectStdBremModel = new G4SeltzerBergerModel();
86  theDirectEMModel=theDirectStdBremModel;
87  theEmModelManagerForFwdModels = new G4EmModelManager();
88  isDirectModelInitialised = false;
90  G4Region* r=0;
91  theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r);
92  // theDirectPenelopeBremModel =0;
93  SetApplyCutInRange(true);
94  highKinEnergy= 1.*GeV;
95  lowKinEnergy = 1.0*keV;
96  lastCZ =0.;
101 }
static G4AdjointGamma * AdjointGamma()
G4ParticleDefinition * theDirectPrimaryPartDef
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
static G4AdjointElectron * AdjointElectron()
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
void SetUseMatrixPerElement(G4bool aBool)
G4VEmModel * theDirectEMModel
G4VEmAdjointModel(const G4String &nam)
void SetUseMatrix(G4bool aBool)
static constexpr double GeV
Definition: G4SIunits.hh:217
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
void SetApplyCutInRange(G4bool aBool)
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel ( )

Definition at line 104 of file G4AdjointBremsstrahlungModel.cc.

105 {if (theDirectStdBremModel) delete theDirectStdBremModel;
106  if (theEmModelManagerForFwdModels) delete theEmModelManagerForFwdModels;
107 }

Member Function Documentation

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

Reimplemented from G4VEmAdjointModel.

Definition at line 400 of file G4AdjointBremsstrahlungModel.cc.

403 { if (!isDirectModelInitialised) {
404  theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
405  isDirectModelInitialised =true;
406  }
407  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
408  DefineCurrentMaterial(aCouple);
409  G4double Cross=0.;
410  lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above
411 
412  if (!IsScatProjToProjCase ){
413  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
414  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
415  if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= CS_biasing_factor*lastCZ*std::log(Emax_proj/Emin_proj);
416  }
417  else {
420  if (Emax_proj>Emin_proj) Cross= lastCZ*std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy));
421 
422  }
423  return Cross;
424 }
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:257
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
G4ParticleDefinition * theDirectPrimaryPartDef
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
G4VEmModel * theDirectEMModel
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
const G4DataVector * Initialise(const G4ParticleDefinition *part, const G4ParticleDefinition *secPart, G4double minSubRange, G4int verb)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
static G4Electron * Electron()
Definition: G4Electron.cc:94
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double currentTcutForDirectSecond
double G4double
Definition: G4Types.hh:76
const G4Material * GetMaterial() const

Here is the call graph for this function:

Here is the caller graph for this function:

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

Reimplemented from G4VEmAdjointModel.

Definition at line 324 of file G4AdjointBremsstrahlungModel.cc.

328 {if (!isDirectModelInitialised) {
329  theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
330  isDirectModelInitialised =true;
331  }
332 /*
333  return DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial,
334  kinEnergyProj,
335  kinEnergyProd);
336 */
338  kinEnergyProj,
339  kinEnergyProd);
340 }
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4DataVector * Initialise(const G4ParticleDefinition *part, const G4ParticleDefinition *secPart, G4double minSubRange, G4int verb)
static G4Electron * Electron()
Definition: G4Electron.cc:94

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 344 of file G4AdjointBremsstrahlungModel.cc.

349 {
350  G4double dCrossEprod=0.;
351  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
352  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
353 
354 
355  //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution
356  //This is what is applied in the discrete standard model before the rejection test that make a correction
357  //The application of the same rejection function is not possible here.
358  //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the
359  // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and
360  // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one.
361  // In the future we plan to use the brem secondary spectra from the G4Penelope implementation
362 
363  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
365  dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV);
366  }
367  return dCrossEprod;
368 
369 }
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:257
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
G4ParticleDefinition * theDirectPrimaryPartDef
G4VEmModel * theDirectEMModel
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
double G4double
Definition: G4Types.hh:76
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

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

Definition at line 373 of file G4AdjointBremsstrahlungModel.cc.

378 {
379  //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor
380  //used in the direct model
381 
382  G4double dCrossEprod=0.;
383 
384  const G4ElementVector* theElementVector = material->GetElementVector();
385  const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
386  G4double dum=0.;
387  G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001;
388  G4double dE=E2-E1;
389  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
390  G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1);
391  G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2);
392  dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE;
393  }
394 
395  return dCrossEprod;
396 
397 }
const double C2
const double C1
std::vector< G4Element * > G4ElementVector
G4ParticleDefinition * theDirectPrimaryPartDef
string material
Definition: eplot.py:19
static const G4double dE
G4VEmModel * theDirectEMModel
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:321
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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

Reimplemented from G4VEmAdjointModel.

Definition at line 426 of file G4AdjointBremsstrahlungModel.cc.

429 {
430  return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
431  lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above
432  return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
433 
434 }
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:257
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4ParticleDefinition * theDirectPrimaryPartDef
virtual G4double GetAdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4VEmModel * theDirectEMModel
static constexpr double MeV
Definition: G4SIunits.hh:214
const G4Material * GetMaterial() const

Here is the call graph for this function:

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

Definition at line 189 of file G4AdjointBremsstrahlungModel.cc.

192 {
193 
194  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
196 
197 
198  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
199  G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
200 
201  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
202  return;
203  }
204 
205  G4double projectileKinEnergy =0.;
206  G4double gammaEnergy=0.;
207  G4double diffCSUsed=0.;
208  if (!IsScatProjToProjCase){
209  gammaEnergy=adjointPrimKinEnergy;
210  G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
211  G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
212  if (Emin>=Emax) return;
213  projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
214  diffCSUsed=CS_biasing_factor*lastCZ/projectileKinEnergy;
215 
216  }
217  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
219  if (Emin>=Emax) return;
220  G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
221  G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
222  projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));
223  gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy;
224  diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy;
225 
226  }
227 
228 
229 
230 
231  //Weight correction
232  //-----------------------
233  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
234  //if this has to be done in the model
235  //For the case of forced interaction this will be done in the PostStepDoIt of the
236  //forced interaction
237  //It is important to set the weight before the vreation of the secondary
238  //
242  }
243  //G4cout<<"Correction factor start in brem model "<<w_corr<<std::endl;
244 
245 
246  //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the one consistent with the direct model
247  //Here we consider the true diffCS as the one obtained by the numericla differentiation over Tcut of the direct CS, corrected by the Migdal term.
248  //Basically any other differential CS diffCS could be used here (example Penelope).
249 
250  G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy);
251  /*G4cout<<"diffCS "<<diffCS <<std::endl;
252  G4cout<<"diffCS_Used "<<diffCSUsed <<std::endl;*/
253  w_corr*=diffCS/diffCSUsed;
254 
255 
256  G4double new_weight = aTrack.GetWeight()*w_corr;
257  /*G4cout<<"New weight brem "<<new_weight<<std::endl;
258  G4cout<<"Weight correction brem "<<w_corr<<std::endl;*/
259  fParticleChange->SetParentWeightByProcess(false);
260  fParticleChange->SetSecondaryWeightByProcess(false);
261  fParticleChange->ProposeParentWeight(new_weight);
262 
263  //Kinematic
264  //---------
266  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
267  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
268  G4double projectileP = std::sqrt(projectileP2);
269 
270 
271  //Use the angular model of the forward model to generate the gamma direction
272  //---------------------------------------------------------------------------
273 //Dum dynamic particle to use the model
274  G4DynamicParticle * aDynPart = new G4DynamicParticle(G4Electron::Electron(),G4ThreeVector(0.,0.,1.)*projectileP);
275 
276  //Get the element from the direct model
278  projectileKinEnergy,currentTcutForDirectSecond);
279  G4int Z=elm->GetZasInt();
280  G4double energy = aDynPart->GetTotalEnergy()-gammaEnergy;
281  G4ThreeVector projectileMomentum =
282  theDirectEMModel->GetAngularDistribution()->SampleDirection(aDynPart,energy,Z,currentMaterial)*projectileP;
283  G4double phi = projectileMomentum.getPhi();
284 
285 /*
286  //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
287  //------------------------------------------------
288  G4double u;
289  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
290 
291  if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
292  else u = - std::log(G4UniformRand()*G4UniformRand())/a2;
293 
294  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
295 
296  G4double sint = std::sin(theta);
297  G4double cost = std::cos(theta);
298 
299  G4double phi = twopi * G4UniformRand() ;
300  G4ThreeVector projectileMomentum;
301  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
302 */
303  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
304  G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
305  G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
306  G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
307  G4double sint1 = std::sqrt(1.-cost1*cost1);
308  projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
309  }
310 
311  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
312 
313  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
314  fParticleChange->ProposeTrackStatus(fStopAndKill);
315  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
316  }
317  else {
318  fParticleChange->ProposeEnergy(projectileKinEnergy);
319  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
320  }
321 }
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
double angle(const Hep3Vector &) const
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetTotalEnergy() const
G4double GetPostStepWeightCorrection()
const G4DynamicParticle * GetDynamicParticle() const
void ProposeParentWeight(G4double finalWeight)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:619
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
int G4int
Definition: G4Types.hh:78
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
G4VEmModel * theDirectEMModel
#define G4UniformRand()
Definition: Randomize.hh:97
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4MaterialCutsCouple * currentCouple
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4bool correct_weight_for_post_step_in_model
G4double GetPDGMass() const
G4double additional_weight_correction_factor_for_post_step_outside_model
G4double energy(const ThreeVector &p, const G4double m)
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Hep3Vector unit() const
double getPhi() const
G4int GetZasInt() const
Definition: G4Element.hh:132
void ProposeEnergy(G4double finalEnergy)
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
static const G4double Emin
static G4Electron * Electron()
Definition: G4Electron.cc:94
static const G4double Emax
G4double currentTcutForDirectSecond
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
static G4AdjointCSManager * GetAdjointCSManager()
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544

Here is the call graph for this function:

Here is the caller graph for this function:

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

Implements G4VEmAdjointModel.

Definition at line 111 of file G4AdjointBremsstrahlungModel.cc.

114 {
115  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
116 
117  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
119 
120 
121  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
122  G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
123 
124  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
125  return;
126  }
127 
128  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
129  IsScatProjToProjCase);
130  //Weight correction
131  //-----------------------
132  CorrectPostStepWeight(fParticleChange,
133  aTrack.GetWeight(),
134  adjointPrimKinEnergy,
135  projectileKinEnergy,
136  IsScatProjToProjCase);
137 
138 
139  //Kinematic
140  //---------
142  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
143  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
144  G4double projectileP = std::sqrt(projectileP2);
145 
146 
147  //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
148  //------------------------------------------------
149  G4double u;
150  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
151 
152  if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
153  else u = - std::log(G4UniformRand()*G4UniformRand())/a2;
154 
155  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
156 
157  G4double sint = std::sin(theta);
158  G4double cost = std::cos(theta);
159 
160  G4double phi = twopi * G4UniformRand() ;
161 
162  G4ThreeVector projectileMomentum;
163  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
164  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
165  G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
166  G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
167  G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
168  G4double sint1 = std::sqrt(1.-cost1*cost1);
169  projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
170 
171  }
172 
173  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
174 
175 
176 
177  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
178  fParticleChange->ProposeTrackStatus(fStopAndKill);
179  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
180  }
181  else {
182  fParticleChange->ProposeEnergy(projectileKinEnergy);
183  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
184 
185  }
186 }
double angle(const Hep3Vector &) const
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetTotalEnergy() const
const G4DynamicParticle * GetDynamicParticle() const
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
static constexpr double twopi
Definition: G4SIunits.hh:76
#define G4UniformRand()
Definition: Randomize.hh:97
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
void RapidSampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
float electron_mass_c2
Definition: hepunit.py:274
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4double GetPDGMass() const
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:


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