| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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::G4AdjointPhotoEle
40 : G4VEmAdjointModel("AdjointPEEffect")
41
42 {
43 SetUseMatrix(false);
44 SetApplyCutInRange(false);
45
46 fAdjEquivDirectPrimPart = G4AdjointGamma::
47 fAdjEquivDirectSecondPart = G4AdjointElectro
48 fDirectPrimaryPart = G4Gamma::Gamma()
49 fSecondPartSameType = false;
50 fDirectModel = new G4PEEffectFl
51 }
52
53 //////////////////////////////////////////////
54 G4AdjointPhotoElectricModel::~G4AdjointPhotoEl
55
56 //////////////////////////////////////////////
57 void G4AdjointPhotoElectricModel::SampleSecond
58 const G4Track& aTrack, G4bool isScatProjToPr
59 G4ParticleChange* fParticleChange)
60 {
61 if(isScatProjToProj)
62 return;
63
64 // Compute the fTotAdjointCS vectors if not
65 // couple and electron energy
66 const G4DynamicParticle* aDynPart = aTrack.G
67 G4double electronEnergy = aDynPart
68 G4ThreeVector electronDirection = aDynPart
69 fPreStepAdjointCS =
70 fTotAdjointCS; // The last computed CS wa
71 AdjointCrossSection(aTrack.GetMaterialCutsCo
72 isScatProjToProj);
73 fPostStepAdjointCS = fTotAdjointCS;
74
75 // Sample element
76 const G4ElementVector* theElementVector =
77 fCurrentMaterial->GetElementVector();
78 size_t nelm = fCurrentMaterial->GetNumb
79 G4double rand_CS = G4UniformRand() * fXsec[n
80 for(fIndexElement = 0; fIndexElement < nelm
81 {
82 if(rand_CS < fXsec[fIndexElement])
83 break;
84 }
85
86 // Sample shell and binding energy
87 G4int nShells = (*theElementVector)[fIndexEl
88 rand_CS = fShellProb[fIndexElement][nS
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)[fInde
97
98 // Sample cos theta
99 // Copy of the G4PEEfectFluoModel cos theta
100 // ElecCosThetaDistribution. This method can
101 // G4PEEffectFluoModel because it is a frien
102 G4double cos_theta = 1.;
103 G4double gamma = 1. + electronEnergy / e
104 if(gamma <= 5.)
105 {
106 G4double beta = std::sqrt(gamma * gamma -
107 G4double b = 0.5 * gamma * (gamma - 1.)
108
109 G4double rndm, term, greject, grejsup;
110 if(gamma < 2.)
111 grejsup = gamma * gamma * (1. + b - beta
112 else
113 grejsup = gamma * gamma * (1. + b + beta
114
115 do
116 {
117 rndm = 1. - 2. * G4UniformRand();
118 cos_theta = (rndm + beta) / (rndm * beta
119 term = 1. - beta * cos_theta;
120 greject = (1. - cos_theta * cos_theta) *
121 // Loop checking, 07-Aug-2015, Vladimir
122 } while(greject < G4UniformRand() * grejsu
123 }
124
125 // direction of the adjoint gamma electron
126 G4double sin_theta = std::sqrt(1. - cos_thet
127 G4double phi = twopi * G4UniformRand()
128 G4double dirx = sin_theta * std::cos(ph
129 G4double diry = sin_theta * std::sin(ph
130 G4double dirz = cos_theta;
131 G4ThreeVector adjoint_gammaDirection(dirx, d
132 adjoint_gammaDirection.rotateUz(electronDire
133
134 // Weight correction
135 CorrectPostStepWeight(fParticleChange, aTrac
136 gammaEnergy, isScatPro
137
138 // Create secondary and modify fParticleChan
139 G4DynamicParticle* anAdjointGamma = new G4Dy
140 G4AdjointGamma::AdjointGamma(), adjoint_ga
141
142 fParticleChange->ProposeTrackStatus(fStopAnd
143 fParticleChange->AddSecondary(anAdjointGamma
144 }
145
146 //////////////////////////////////////////////
147 void G4AdjointPhotoElectricModel::CorrectPostS
148 G4ParticleChange* fParticleChange, G4double
149 G4double adjointPrimKinEnergy, G4double proj
150 {
151 G4double new_weight = old_weight;
152
153 G4double w_corr =
154 G4AdjointCSManager::GetAdjointCSManager()-
155 fFactorCSBiasing;
156 w_corr *= fPostStepAdjointCS / fPreStepAdjoi
157
158 new_weight *= w_corr * projectileKinEnergy /
159 fParticleChange->SetParentWeightByProcess(fa
160 fParticleChange->SetSecondaryWeightByProcess
161 fParticleChange->ProposeParentWeight(new_wei
162 }
163
164 //////////////////////////////////////////////
165 G4double G4AdjointPhotoElectricModel::AdjointC
166 const G4MaterialCutsCouple* aCouple, G4doubl
167 G4bool isScatProjToProj)
168 {
169 if(isScatProjToProj)
170 return 0.;
171
172 G4double totBiasedAdjointCS = 0.;
173 if(aCouple != fCurrentCouple || fCurrenteEne
174 {
175 fTotAdjointCS = 0.;
176 DefineCurrentMaterialAndElectronEnergy(aCo
177 const G4ElementVector* theElementVector =
178 fCurrentMaterial->GetElementVector();
179 const G4double* theAtomNumDensityVector =
180 fCurrentMaterial->GetVecNbOfAtomsPerVolu
181 size_t nelm = fCurrentMaterial->GetNumberO
182 for(fIndexElement = 0; fIndexElement < nel
183 {
184 fTotAdjointCS += AdjointCrossSectionPerA
185 (*theElementVector)[f
186 theAtomNumDensityVector
187 fXsec[fIndexElement] = fTotAdjointCS;
188 }
189
190 totBiasedAdjointCS = std::min(fTotAdjointC
191 fFactorCSBiasing = totBiasedAdjointCS /
192 }
193 return totBiasedAdjointCS;
194 }
195
196 //////////////////////////////////////////////
197 G4double G4AdjointPhotoElectricModel::AdjointC
198 const G4Element* anElement, G4double electro
199 {
200 G4int nShells = anElement->GetNbOfAto
201 G4double Z = anElement->GetZ();
202 G4double gammaEnergy = electronEnergy + anEl
203 G4double CS = fDirectModel->Compute
204 G4Gamma::Gamma(), gammaEnergy, Z, 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
212 G4double Bi = anElement->GetAtomicShell(i
213 if(electronEnergy < Bi1 - Bi)
214 {
215 gammaEnergy = electronEnergy + Bi;
216 CS = fDirectModel->ComputeCross
217
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::DefineCurren
229 const G4MaterialCutsCouple* couple, G4double
230 {
231 fCurrentCouple = const_cast<G4MaterialCuts
232 fCurrentMaterial = const_cast<G4Material*>(c
233 fCurrenteEnergy = anEnergy;
234 fDirectModel->SetCurrentCouple(couple);
235 }
236