| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26
27 #include "G4AdjointPhotoElectricModel.hh"
28
29 #include "G4AdjointCSManager.hh"
30 #include "G4AdjointElectron.hh"
31 #include "G4AdjointGamma.hh"
32 #include "G4Gamma.hh"
33 #include "G4ParticleChange.hh"
34 #include "G4PEEffectFluoModel.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "G4TrackStatus.hh"
37
38 ////////////////////////////////////////////////////////////////////////////////
39 G4AdjointPhotoElectricModel::G4AdjointPhotoElectricModel()
40 : G4VEmAdjointModel("AdjointPEEffect")
41
42 {
43 SetUseMatrix(false);
44 SetApplyCutInRange(false);
45
46 fAdjEquivDirectPrimPart = G4AdjointGamma::AdjointGamma();
47 fAdjEquivDirectSecondPart = G4AdjointElectron::AdjointElectron();
48 fDirectPrimaryPart = G4Gamma::Gamma();
49 fSecondPartSameType = false;
50 fDirectModel = new G4PEEffectFluoModel();
51 }
52
53 ////////////////////////////////////////////////////////////////////////////////
54 G4AdjointPhotoElectricModel::~G4AdjointPhotoElectricModel() {}
55
56 ////////////////////////////////////////////////////////////////////////////////
57 void G4AdjointPhotoElectricModel::SampleSecondaries(
58 const G4Track& aTrack, G4bool isScatProjToProj,
59 G4ParticleChange* fParticleChange)
60 {
61 if(isScatProjToProj)
62 return;
63
64 // Compute the fTotAdjointCS vectors if not already done for the current
65 // couple and electron energy
66 const G4DynamicParticle* aDynPart = aTrack.GetDynamicParticle();
67 G4double electronEnergy = aDynPart->GetKineticEnergy();
68 G4ThreeVector electronDirection = aDynPart->GetMomentumDirection();
69 fPreStepAdjointCS =
70 fTotAdjointCS; // The last computed CS was at pre step point
71 AdjointCrossSection(aTrack.GetMaterialCutsCouple(), electronEnergy,
72 isScatProjToProj);
73 fPostStepAdjointCS = fTotAdjointCS;
74
75 // Sample element
76 const G4ElementVector* theElementVector =
77 fCurrentMaterial->GetElementVector();
78 size_t nelm = fCurrentMaterial->GetNumberOfElements();
79 G4double rand_CS = G4UniformRand() * fXsec[nelm - 1];
80 for(fIndexElement = 0; fIndexElement < nelm - 1; ++fIndexElement)
81 {
82 if(rand_CS < fXsec[fIndexElement])
83 break;
84 }
85
86 // Sample shell and binding energy
87 G4int nShells = (*theElementVector)[fIndexElement]->GetNbOfAtomicShells();
88 rand_CS = fShellProb[fIndexElement][nShells - 1] * G4UniformRand();
89 G4int i;
90 for(i = 0; i < nShells - 1; ++i)
91 {
92 if(rand_CS < fShellProb[fIndexElement][i])
93 break;
94 }
95 G4double gammaEnergy =
96 electronEnergy + (*theElementVector)[fIndexElement]->GetAtomicShell(i);
97
98 // Sample cos theta
99 // Copy of the G4PEEfectFluoModel cos theta sampling method
100 // ElecCosThetaDistribution. This method cannot be used directly from
101 // G4PEEffectFluoModel because it is a friend method.
102 G4double cos_theta = 1.;
103 G4double gamma = 1. + electronEnergy / electron_mass_c2;
104 if(gamma <= 5.)
105 {
106 G4double beta = std::sqrt(gamma * gamma - 1.) / gamma;
107 G4double b = 0.5 * gamma * (gamma - 1.) * (gamma - 2.);
108
109 G4double rndm, term, greject, grejsup;
110 if(gamma < 2.)
111 grejsup = gamma * gamma * (1. + b - beta * b);
112 else
113 grejsup = gamma * gamma * (1. + b + beta * b);
114
115 do
116 {
117 rndm = 1. - 2. * G4UniformRand();
118 cos_theta = (rndm + beta) / (rndm * beta + 1.);
119 term = 1. - beta * cos_theta;
120 greject = (1. - cos_theta * cos_theta) * (1. + b * term) / (term * term);
121 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
122 } while(greject < G4UniformRand() * grejsup);
123 }
124
125 // direction of the adjoint gamma electron
126 G4double sin_theta = std::sqrt(1. - cos_theta * cos_theta);
127 G4double phi = twopi * G4UniformRand();
128 G4double dirx = sin_theta * std::cos(phi);
129 G4double diry = sin_theta * std::sin(phi);
130 G4double dirz = cos_theta;
131 G4ThreeVector adjoint_gammaDirection(dirx, diry, dirz);
132 adjoint_gammaDirection.rotateUz(electronDirection);
133
134 // Weight correction
135 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,
136 gammaEnergy, isScatProjToProj);
137
138 // Create secondary and modify fParticleChange
139 G4DynamicParticle* anAdjointGamma = new G4DynamicParticle(
140 G4AdjointGamma::AdjointGamma(), adjoint_gammaDirection, gammaEnergy);
141
142 fParticleChange->ProposeTrackStatus(fStopAndKill);
143 fParticleChange->AddSecondary(anAdjointGamma);
144 }
145
146 ////////////////////////////////////////////////////////////////////////////////
147 void G4AdjointPhotoElectricModel::CorrectPostStepWeight(
148 G4ParticleChange* fParticleChange, G4double old_weight,
149 G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool)
150 {
151 G4double new_weight = old_weight;
152
153 G4double w_corr =
154 G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection() /
155 fFactorCSBiasing;
156 w_corr *= fPostStepAdjointCS / fPreStepAdjointCS;
157
158 new_weight *= w_corr * projectileKinEnergy / adjointPrimKinEnergy;
159 fParticleChange->SetParentWeightByProcess(false);
160 fParticleChange->SetSecondaryWeightByProcess(false);
161 fParticleChange->ProposeParentWeight(new_weight);
162 }
163
164 ////////////////////////////////////////////////////////////////////////////////
165 G4double G4AdjointPhotoElectricModel::AdjointCrossSection(
166 const G4MaterialCutsCouple* aCouple, G4double electronEnergy,
167 G4bool isScatProjToProj)
168 {
169 if(isScatProjToProj)
170 return 0.;
171
172 G4double totBiasedAdjointCS = 0.;
173 if(aCouple != fCurrentCouple || fCurrenteEnergy != electronEnergy)
174 {
175 fTotAdjointCS = 0.;
176 DefineCurrentMaterialAndElectronEnergy(aCouple, electronEnergy);
177 const G4ElementVector* theElementVector =
178 fCurrentMaterial->GetElementVector();
179 const G4double* theAtomNumDensityVector =
180 fCurrentMaterial->GetVecNbOfAtomsPerVolume();
181 size_t nelm = fCurrentMaterial->GetNumberOfElements();
182 for(fIndexElement = 0; fIndexElement < nelm; ++fIndexElement)
183 {
184 fTotAdjointCS += AdjointCrossSectionPerAtom(
185 (*theElementVector)[fIndexElement], electronEnergy) *
186 theAtomNumDensityVector[fIndexElement];
187 fXsec[fIndexElement] = fTotAdjointCS;
188 }
189
190 totBiasedAdjointCS = std::min(fTotAdjointCS, 0.01);
191 fFactorCSBiasing = totBiasedAdjointCS / fTotAdjointCS;
192 }
193 return totBiasedAdjointCS;
194 }
195
196 ////////////////////////////////////////////////////////////////////////////////
197 G4double G4AdjointPhotoElectricModel::AdjointCrossSectionPerAtom(
198 const G4Element* anElement, G4double electronEnergy)
199 {
200 G4int nShells = anElement->GetNbOfAtomicShells();
201 G4double Z = anElement->GetZ();
202 G4double gammaEnergy = electronEnergy + anElement->GetAtomicShell(0);
203 G4double CS = fDirectModel->ComputeCrossSectionPerAtom(
204 G4Gamma::Gamma(), gammaEnergy, Z, 0., 0., 0.);
205 G4double adjointCS = 0.;
206 if(CS > 0.)
207 adjointCS += CS / gammaEnergy;
208 fShellProb[fIndexElement][0] = adjointCS;
209 for(G4int i = 1; i < nShells; ++i)
210 {
211 G4double Bi1 = anElement->GetAtomicShell(i - 1);
212 G4double Bi = anElement->GetAtomicShell(i);
213 if(electronEnergy < Bi1 - Bi)
214 {
215 gammaEnergy = electronEnergy + Bi;
216 CS = fDirectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
217 gammaEnergy, Z, 0., 0., 0.);
218 if(CS > 0.)
219 adjointCS += CS / gammaEnergy;
220 }
221 fShellProb[fIndexElement][i] = adjointCS;
222 }
223 adjointCS *= electronEnergy;
224 return adjointCS;
225 }
226
227 ////////////////////////////////////////////////////////////////////////////////
228 void G4AdjointPhotoElectricModel::DefineCurrentMaterialAndElectronEnergy(
229 const G4MaterialCutsCouple* couple, G4double anEnergy)
230 {
231 fCurrentCouple = const_cast<G4MaterialCutsCouple*>(couple);
232 fCurrentMaterial = const_cast<G4Material*>(couple->GetMaterial());
233 fCurrenteEnergy = anEnergy;
234 fDirectModel->SetCurrentCouple(couple);
235 }
236