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