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

#include <G4AdjointhIonisationModel.hh>

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

Public Member Functions

 G4AdjointhIonisationModel (G4ParticleDefinition *projectileDefinition)
 
virtual ~G4AdjointhIonisationModel ()
 
virtual void SampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
 
void RapidSampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
 
virtual G4double DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
 
virtual G4double AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
 
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase (G4double PrimAdjEnergy)
 
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase (G4double PrimAdjEnergy, G4double Tcut=0)
 
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase (G4double PrimAdjEnergy)
 
virtual G4double GetSecondAdjEnergyMinForProdToProjCase (G4double PrimAdjEnergy)
 
- Public Member Functions inherited from G4VEmAdjointModel
 G4VEmAdjointModel (const G4String &nam)
 
virtual ~G4VEmAdjointModel ()
 
virtual G4double GetAdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
 
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)
 
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 72 of file G4AdjointhIonisationModel.hh.

Constructor & Destructor Documentation

G4AdjointhIonisationModel::G4AdjointhIonisationModel ( G4ParticleDefinition *  projectileDefinition)

Definition at line 46 of file G4AdjointhIonisationModel.cc.

46  :
47  G4VEmAdjointModel("Adjoint_hIonisation")
48 {
49 
50 
51 
52  UseMatrix =true;
53  UseMatrixPerElement = true;
54  ApplyCutInRange = true;
55  UseOnlyOneMatrixForAllElements = true;
56  CS_biasing_factor =1.;
57  second_part_of_same_type =false;
58 
59  //The direct EM Modfel is taken has BetheBloch it is only used for the computation
60  // of the differential cross section.
61  //The Bragg model could be used as an alternative as it offers the same differential cross section
62 
63  theDirectEMModel = new G4BetheBlochModel(projectileDefinition);
64  theBraggDirectEMModel = new G4BraggModel(projectileDefinition);
65  theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
66 
67  theDirectPrimaryPartDef = projectileDefinition;
68  theAdjEquivOfDirectPrimPartDef = 0;
69  if (projectileDefinition == G4Proton::Proton()) {
70  theAdjEquivOfDirectPrimPartDef = G4AdjointProton::AdjointProton();
71 
72  }
73 
74  DefineProjectileProperty();
75 }
G4VEmAdjointModel::theDirectPrimaryPartDef
G4ParticleDefinition * theDirectPrimaryPartDef
Definition: G4VEmAdjointModel.hh:297
G4AdjointElectron::AdjointElectron
static G4AdjointElectron * AdjointElectron()
Definition: G4AdjointElectron.cc:100
G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
Definition: G4VEmAdjointModel.hh:295
G4VEmAdjointModel::theDirectEMModel
G4VEmModel * theDirectEMModel
Definition: G4VEmAdjointModel.hh:258
G4VEmAdjointModel::G4VEmAdjointModel
G4VEmAdjointModel(const G4String &nam)
Definition: G4VEmAdjointModel.cc:41
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4VEmAdjointModel::UseOnlyOneMatrixForAllElements
G4bool UseOnlyOneMatrixForAllElements
Definition: G4VEmAdjointModel.hh:333
G4AdjointProton::AdjointProton
static G4AdjointProton * AdjointProton()
Definition: G4AdjointProton.cc:84
G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
Definition: G4VEmAdjointModel.hh:296

Here is the call graph for this function:

G4AdjointhIonisationModel::~G4AdjointhIonisationModel ( )
virtual

Definition at line 78 of file G4AdjointhIonisationModel.cc.

79 {;}

Member Function Documentation

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

Reimplemented from G4VEmAdjointModel.

Definition at line 438 of file G4AdjointhIonisationModel.cc.

441 {
442  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
443  DefineCurrentMaterial(aCouple);
444 
445 
446  G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
447 
448  if (!IsScatProjToProjCase ){
449  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
450  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
451  if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
452  Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
453  }
454  else Cross=0.;
455 
456 
457 
458 
459 
460 
461  }
462  else {
463  G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
464  G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
465  G4double diff1=Emin_proj-primEnergy;
466  G4double diff2=Emax_proj-primEnergy;
467  G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
468  //G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
469  G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy;
470  Cross*=(t1+t2);
471 
472 
473  }
474  lastCS =Cross;
475  return Cross;
476 }
G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointhIonisationModel.cc:479
G4VEmAdjointModel::currentMaterial
G4Material * currentMaterial
Definition: G4VEmAdjointModel.hh:306
G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Definition: G4AdjointhIonisationModel.cc:486
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4VEmAdjointModel::DefineCurrentMaterial
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
Definition: G4VEmAdjointModel.cc:693
G4VEmAdjointModel::AdjointCrossSection
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
Definition: G4VEmAdjointModel.cc:59
plottest35.t1
tuple t1
Definition: plottest35.py:33
G4VEmAdjointModel::currentTcutForDirectSecond
G4double currentTcutForDirectSecond
Definition: G4VEmAdjointModel.hh:311
G4double
double G4double
Definition: G4Types.hh:76
G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointhIonisationModel.cc:496
G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointhIonisationModel.cc:491

Here is the call graph for this function:

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

Reimplemented from G4VEmAdjointModel.

Definition at line 306 of file G4AdjointhIonisationModel.cc.

311 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
312 
313 
314 
315  G4double dSigmadEprod=0;
316  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
317  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
318 
319 
320  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
321  G4double Tmax=kinEnergyProj;
322 
323  G4double E1=kinEnergyProd;
324  G4double E2=kinEnergyProd*1.000001;
325  G4double dE=(E2-E1);
326  G4double sigma1,sigma2;
327  if (kinEnergyProj >2.*MeV){
328  sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
329  sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
330  }
331  else {
332  sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
333  sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
334  }
335 
336 
337  dSigmadEprod=(sigma1-sigma2)/dE;
338  if (dSigmadEprod>1.) {
339  G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
340  G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
341  G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
342 
343  }
344 
345 
346 
347  //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
348  //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
349  //to test the rejection of a secondary
350  //-------------------------
351 
352  //Source code taken from G4BetheBlochModel::SampleSecondaries
353 
354  G4double deltaKinEnergy = kinEnergyProd;
355 
356  //Part of the taken code
357  //----------------------
358 
359 
360 
361  // projectile formfactor - suppresion of high energy
362  // delta-electron production at high energy
363  G4double x = formfact*deltaKinEnergy;
364  if(x > 1.e-6) {
365 
366 
367  G4double totEnergy = kinEnergyProj + mass;
368  G4double etot2 = totEnergy*totEnergy;
369  G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
370  G4double f;
371  G4double f1 = 0.0;
372  f = 1.0 - beta2*deltaKinEnergy/Tmax;
373  if( 0.5 == spin ) {
374  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
375  f += f1;
376  }
377  G4double x1 = 1.0 + x;
378  G4double gg = 1.0/(x1*x1);
379  if( 0.5 == spin ) {
380  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
381  gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
382  }
383  if(gg > 1.0) {
384  G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
385  << G4endl;
386  gg=1.;
387  }
388  //G4cout<<"gg"<<gg<<G4endl;
389  dSigmadEprod*=gg;
390  }
391 
392  }
393 
394  return dSigmadEprod;
395 }
G4VEmAdjointModel::theDirectPrimaryPartDef
G4ParticleDefinition * theDirectPrimaryPartDef
Definition: G4VEmAdjointModel.hh:297
test.x
tuple x
Definition: test.py:50
g
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
dE
static const G4double dE
Definition: G4PhotoNuclearCrossSection.cc:47
G4VEmAdjointModel::theDirectEMModel
G4VEmModel * theDirectEMModel
Definition: G4VEmAdjointModel.hh:258
G4cout
G4GLOB_DLL std::ostream G4cout
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:60
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4VEmModel::ComputeCrossSectionPerAtom
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:321
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointhIonisationModel.cc:496
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262
G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointhIonisationModel.cc:491

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 491 of file G4AdjointhIonisationModel.cc.

492 { return HighEnergyLimit;
493 }

Here is the caller graph for this function:

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 479 of file G4AdjointhIonisationModel.cc.

480 {
481  G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
482  return Tmax;
483 }
G4double
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 496 of file G4AdjointhIonisationModel.cc.

497 { G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
498  return Tmin;
499 }
G4double
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase ( G4double  PrimAdjEnergy,
G4double  Tcut = 0 
)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 486 of file G4AdjointhIonisationModel.cc.

487 { return PrimAdjEnergy+Tcut;
488 }

Here is the caller graph for this function:

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

Definition at line 160 of file G4AdjointhIonisationModel.cc.

163 {
164 
165  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
166  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
167 
168 
169  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
170  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
171 
172  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
173  return;
174  }
175 
176  G4double projectileKinEnergy =0.;
177  G4double eEnergy=0.;
178  G4double newCS=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
179  if (!IsScatProjToProjCase){//1/E^2 distribution
180 
181  eEnergy=adjointPrimKinEnergy;
182  G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
183  G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
184  if (Emin>=Emax) return;
185  G4double a=1./Emax;
186  G4double b=1./Emin;
187  newCS=newCS*(b-a)/eEnergy;
188 
189  projectileKinEnergy =1./(b- (b-a)*G4UniformRand());
190 
191 
192  }
193  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
194  G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
195  if (Emin>=Emax) return;
196  G4double diff1=Emin-adjointPrimKinEnergy;
197  G4double diff2=Emax-adjointPrimKinEnergy;
198 
199  G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
200  G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
201  /*G4double f31=diff1/Emin;
202  G4double f32=diff2/Emax/f31;*/
203  G4double t3=2.*std::log(Emax/Emin);
204  G4double sum_t=t1+t2+t3;
205  newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
206  G4double t=G4UniformRand()*sum_t;
207  if (t <=t1 ){
208  G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ;
209  projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
210 
211  }
212  else if (t <=t2 ) {
213  G4double q= G4UniformRand()*t2/adjointPrimKinEnergy;
214  projectileKinEnergy =1./(1./Emin-q);
215  }
216  else {
217  projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
218 
219  }
220  eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
221 
222 
223  }
224 
225 
226 
227  G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy;
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  G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
235 
236  //G4cout<<w_corr<<G4endl;
237  w_corr*=newCS/lastCS;
238  //G4cout<<w_corr<<G4endl;
239  //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
240  //Here we consider the true diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS
241 
242  G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1);
243  w_corr*=diffCS/diffCS_perAtom_Used;
244  //G4cout<<w_corr<<G4endl;
245 
246  G4double new_weight = aTrack.GetWeight()*w_corr;
247  fParticleChange->SetParentWeightByProcess(false);
248  fParticleChange->SetSecondaryWeightByProcess(false);
249  fParticleChange->ProposeParentWeight(new_weight);
250 
251 
252 
253 
254  //Kinematic:
255  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
256  // him part of its energy
257  //----------------------------------------------------------------------------------------
258 
259  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
260  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
261  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
262 
263 
264 
265  //Companion
266  //-----------
267  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
268  if (IsScatProjToProjCase) {
269  companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
270  }
271  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
272  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
273 
274 
275  //Projectile momentum
276  //--------------------
277  G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
278  G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
279  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
280  G4double phi =G4UniformRand()*2.*3.1415926;
281  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
282  projectileMomentum.rotateUz(dir_parallel);
283 
284 
285 
286  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
287  fParticleChange->ProposeTrackStatus(fStopAndKill);
288  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
289  //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
290  }
291  else {
292  fParticleChange->ProposeEnergy(projectileKinEnergy);
293  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
294  }
295 
296 
297 
298 
299 
300 
301 
302 }
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4AdjointCSManager::GetPostStepWeightCorrection
G4double GetPostStepWeightCorrection()
Definition: G4AdjointCSManager.cc:527
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4VParticleChange::ProposeParentWeight
void ProposeParentWeight(G4double finalWeight)
G4Track::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointhIonisationModel.cc:479
G4VEmAdjointModel::currentMaterial
G4Material * currentMaterial
Definition: G4VEmAdjointModel.hh:306
G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
Definition: G4VEmAdjointModel.hh:295
G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Definition: G4AdjointhIonisationModel.cc:486
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
test.b
tuple b
Definition: test.py:12
G4VParticleChange::SetSecondaryWeightByProcess
void SetSecondaryWeightByProcess(G4bool)
G4VParticleChange::SetParentWeightByProcess
void SetParentWeightByProcess(G4bool)
Definition: G4VParticleChange.cc:522
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4VEmAdjointModel::DefineCurrentMaterial
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
Definition: G4VEmAdjointModel.cc:693
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
plottest35.t1
tuple t1
Definition: plottest35.py:33
G4ParticleChange::ProposeEnergy
void ProposeEnergy(G4double finalEnergy)
fStopAndKill
Definition: G4TrackStatus.hh:56
G4ParticleChange::AddSecondary
void AddSecondary(G4Track *aSecondary)
Definition: G4ParticleChange.cc:218
G4Track::GetWeight
G4double GetWeight() const
Emin
static const G4double Emin
Definition: G4PhotoNuclearCrossSection.cc:50
Emax
static const G4double Emax
Definition: G4PhotoNuclearCrossSection.cc:51
G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
Definition: G4VEmAdjointModel.hh:296
G4VEmAdjointModel::currentTcutForDirectSecond
G4double currentTcutForDirectSecond
Definition: G4VEmAdjointModel.hh:311
G4ParticleChange::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double
double G4double
Definition: G4Types.hh:76
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointhIonisationModel.cc:496
G4AdjointCSManager::GetAdjointCSManager
static G4AdjointCSManager * GetAdjointCSManager()
Definition: G4AdjointCSManager.cc:61
G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
Definition: G4AdjointhIonisationModel.cc:306
G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4AdjointhIonisationModel.cc:491

Here is the call graph for this function:

Here is the caller graph for this function:

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

Implements G4VEmAdjointModel.

Definition at line 84 of file G4AdjointhIonisationModel.cc.

87 {
88  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
89 
90  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
91 
92  //Elastic inverse scattering
93  //---------------------------------------------------------
94  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
95  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
96 
97  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
98  return;
99  }
100 
101  //Sample secondary energy
102  //-----------------------
103  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
104  CorrectPostStepWeight(fParticleChange,
105  aTrack.GetWeight(),
106  adjointPrimKinEnergy,
107  projectileKinEnergy,
108  IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
109 
110 
111  //Kinematic:
112  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
113  // him part of its energy
114  //----------------------------------------------------------------------------------------
115 
116  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
117  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
118  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
119 
120 
121 
122  //Companion
123  //-----------
124  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
125  if (IsScatProjToProjCase) {
126  companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
127  }
128  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
129  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
130 
131 
132  //Projectile momentum
133  //--------------------
134  G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
135  G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
136  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
137  G4double phi =G4UniformRand()*2.*3.1415926;
138  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
139  projectileMomentum.rotateUz(dir_parallel);
140 
141 
142 
143  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
144  fParticleChange->ProposeTrackStatus(fStopAndKill);
145  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
146  //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
147  }
148  else {
149  fParticleChange->ProposeEnergy(projectileKinEnergy);
150  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
151  }
152 
153 
154 
155 
156 }
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
Definition: G4VEmAdjointModel.cc:459
G4AdjointhIonisationModel::RapidSampleSecondaries
void RapidSampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
Definition: G4AdjointhIonisationModel.cc:160
G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
Definition: G4VEmAdjointModel.hh:295
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
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::theAdjEquivOfDirectSecondPartDef
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
Definition: G4VEmAdjointModel.hh:296
G4ParticleChange::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double
double G4double
Definition: G4Types.hh:76
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VEmAdjointModel::CorrectPostStepWeight
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
Definition: G4VEmAdjointModel.cc:634

Here is the call graph for this function:

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