52 post_step_AdjointCS =0.;
53 pre_step_AdjointCS =0.;
54 totBiasedAdjointCS =0.;
72 G4bool IsScatProjToProjCase,
74 {
if (IsScatProjToProjCase) return ;
82 pre_step_AdjointCS = totAdjointCS;
84 post_step_AdjointCS = totAdjointCS;
94 for (index_element=0; index_element<nelm-1; index_element++){
95 if (rand_CS<xsec[index_element])
break;
100 G4int nShells = (*theElementVector)[index_element]->GetNbOfAtomicShells();
101 rand_CS= shell_prob[index_element][nShells-1]*
G4UniformRand();
103 for (i=0; i<nShells-1; i++){
104 if (rand_CS<shell_prob[index_element][i])
break;
106 G4double gammaEnergy= electronEnergy+(*theElementVector)[index_element]->GetAtomicShell(i);
117 G4double beta = std::sqrt(gamma*gamma-1.)/gamma;
118 G4double b = 0.5*gamma*(gamma-1.)*(gamma-2);
121 if (gamma < 2.) grejsup = gamma*gamma*(1.+b-beta*
b);
122 else grejsup = gamma*gamma*(1.+b+beta*
b);
125 cos_theta = (rndm+beta)/(rndm*beta+1.);
126 term = 1.-beta*cos_theta;
127 greject = (1.-cos_theta*cos_theta)*(1.+b*term)/(term*term);
135 G4double sin_theta = std::sqrt(1.-cos_theta*cos_theta);
137 G4double dirx = sin_theta*std::cos(Phi),diry = sin_theta*std::sin(Phi),dirz = cos_theta;
139 adjoint_gammaDirection.
rotateUz(electronDirection);
177 w_corr*=post_step_AdjointCS/pre_step_AdjointCS;
181 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
192 G4bool IsScatProjToProjCase)
196 if (IsScatProjToProjCase)
return 0.;
199 if (aCouple !=
currentCouple || current_eEnergy !=electronEnergy) {
201 DefineCurrentMaterialAndElectronEnergy(aCouple, electronEnergy);
205 for (index_element=0;index_element<nelm;index_element++){
207 totAdjointCS +=
AdjointCrossSectionPerAtom((*theElementVector)[index_element],electronEnergy)*theAtomNumDensityVector[index_element];
208 xsec[index_element] = totAdjointCS;
211 totBiasedAdjointCS=
std::min(totAdjointCS,0.01);
213 factorCSBiasing = totBiasedAdjointCS/totAdjointCS;
214 lastCS=totBiasedAdjointCS;
218 return totBiasedAdjointCS;
227 G4bool IsScatProjToProjCase)
239 G4double gammaEnergy = electronEnergy+B0;
242 if (CS >0) adjointCS += CS/gammaEnergy;
243 shell_prob[index_element][0] = adjointCS;
244 for (i=1;i<nShells;i++){
249 if (electronEnergy <Bi_-Bi) {
250 gammaEnergy = electronEnergy+Bi;
253 if (CS>0) adjointCS +=CS/gammaEnergy;
255 shell_prob[index_element][i] = adjointCS;
258 adjointCS*=electronEnergy;
270 current_eEnergy = anEnergy;
static G4AdjointGamma * AdjointGamma()
G4int GetNbOfAtomicShells() const
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
std::vector< G4Element * > G4ElementVector
G4double GetKineticEnergy() const
G4double GetPostStepWeightCorrection()
const G4DynamicParticle * GetDynamicParticle() const
G4ParticleDefinition * theDirectPrimaryPartDef
DOUBLE PRECISION function rndm(RDUMMY)
~G4AdjointPhotoElectricModel()
void ProposeParentWeight(G4double finalWeight)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
static G4AdjointElectron * AdjointElectron()
virtual G4double GetAdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
const G4ElementVector * GetElementVector() const
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
const G4double * GetVecNbOfAtomsPerVolume() const
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double, G4double)
G4double AdjointCrossSectionPerAtom(const G4Element *anElement, G4double electronEnergy)
void SetUseMatrix(G4bool aBool)
G4MaterialCutsCouple * currentCouple
virtual void SampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
void SetCurrentCouple(const G4MaterialCutsCouple *)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4AdjointPhotoElectricModel()
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
size_t currentCoupleIndex
size_t GetNumberOfElements() const
G4double GetAtomicShell(G4int index) const
void SetApplyCutInRange(G4bool aBool)
void ProposeTrackStatus(G4TrackStatus status)
static G4AdjointCSManager * GetAdjointCSManager()
size_t currentMaterialIndex
const G4Material * GetMaterial() const
G4bool second_part_of_same_type