| 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 // $Id: G4LowEnergyPolarizedCompton.cc 107396 2017-11-10 08:28:08Z gcosmo $
>> 28 // GEANT4 tag $Name: $
27 // 29 //
28 // ------------------------------------------- 30 // ------------------------------------------------------------
29 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file
30 // CERN Geneva Switzerland 32 // CERN Geneva Switzerland
31 // 33 //
32 34
33 // --------- G4LowEnergyPolarizedCompton class 35 // --------- G4LowEnergyPolarizedCompton class -----
34 // 36 //
35 // by G.Depaola & F.Longo (21 may 20 37 // by G.Depaola & F.Longo (21 may 2001)
36 // 38 //
37 // 21 May 2001 - MGP Modified to inherit 39 // 21 May 2001 - MGP Modified to inherit from G4VDiscreteProcess
38 // Applies same algorit 40 // Applies same algorithm as LowEnergyCompton
39 // if the incoming phot 41 // if the incoming photon is not polarised
40 // Temporary protection 42 // Temporary protection to avoid crash in the case
41 // of polarisation || i 43 // of polarisation || incident photon direction
42 // 44 //
43 // 17 October 2001 - F.Longo - Revised accordi 45 // 17 October 2001 - F.Longo - Revised according to a design iteration
44 // 46 //
45 // 21 February 2002 - F.Longo Revisions with A 47 // 21 February 2002 - F.Longo Revisions with A.Zoglauer and G.Depaola
46 // - better descrip 48 // - better description of parallelism
47 // - system of ref 49 // - system of ref change method improved
48 // 22 January 2003 - V.Ivanchenko Cut per reg 50 // 22 January 2003 - V.Ivanchenko Cut per region
49 // 24 April 2003 - V.Ivanchenko Cut per reg 51 // 24 April 2003 - V.Ivanchenko Cut per region mfpt
50 // 52 //
51 // 53 //
52 // ******************************************* 54 // ************************************************************
53 // 55 //
54 // Corrections by Rui Curado da Silva (2000) 56 // Corrections by Rui Curado da Silva (2000)
55 // New Implementation by G.Depaola & F.Longo 57 // New Implementation by G.Depaola & F.Longo
56 // 58 //
57 // - sampling of phi 59 // - sampling of phi
58 // - polarization of scattered photon 60 // - polarization of scattered photon
59 // 61 //
60 // ------------------------------------------- 62 // --------------------------------------------------------------
61 63
62 #include "G4LowEnergyPolarizedCompton.hh" 64 #include "G4LowEnergyPolarizedCompton.hh"
63 #include "Randomize.hh" 65 #include "Randomize.hh"
64 #include "G4PhysicalConstants.hh" 66 #include "G4PhysicalConstants.hh"
65 #include "G4SystemOfUnits.hh" 67 #include "G4SystemOfUnits.hh"
66 #include "G4ParticleDefinition.hh" 68 #include "G4ParticleDefinition.hh"
67 #include "G4Track.hh" 69 #include "G4Track.hh"
68 #include "G4Step.hh" 70 #include "G4Step.hh"
69 #include "G4ForceCondition.hh" 71 #include "G4ForceCondition.hh"
70 #include "G4Gamma.hh" 72 #include "G4Gamma.hh"
71 #include "G4Electron.hh" 73 #include "G4Electron.hh"
72 #include "G4DynamicParticle.hh" 74 #include "G4DynamicParticle.hh"
73 #include "G4VParticleChange.hh" 75 #include "G4VParticleChange.hh"
74 #include "G4ThreeVector.hh" 76 #include "G4ThreeVector.hh"
75 #include "G4RDVCrossSectionHandler.hh" 77 #include "G4RDVCrossSectionHandler.hh"
76 #include "G4RDCrossSectionHandler.hh" 78 #include "G4RDCrossSectionHandler.hh"
77 #include "G4RDVEMDataSet.hh" 79 #include "G4RDVEMDataSet.hh"
78 #include "G4RDCompositeEMDataSet.hh" 80 #include "G4RDCompositeEMDataSet.hh"
79 #include "G4RDVDataSetAlgorithm.hh" 81 #include "G4RDVDataSetAlgorithm.hh"
80 #include "G4RDLogLogInterpolation.hh" 82 #include "G4RDLogLogInterpolation.hh"
81 #include "G4RDVRangeTest.hh" 83 #include "G4RDVRangeTest.hh"
82 #include "G4RDRangeTest.hh" 84 #include "G4RDRangeTest.hh"
83 #include "G4MaterialCutsCouple.hh" 85 #include "G4MaterialCutsCouple.hh"
84 86
85 // constructor 87 // constructor
86 88
87 G4LowEnergyPolarizedCompton::G4LowEnergyPolari 89 G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton(const G4String& processName)
88 : G4VDiscreteProcess(processName), 90 : G4VDiscreteProcess(processName),
89 lowEnergyLimit (250*eV), // i 91 lowEnergyLimit (250*eV), // initialization
90 highEnergyLimit(100*GeV), 92 highEnergyLimit(100*GeV),
91 intrinsicLowEnergyLimit(10*eV), 93 intrinsicLowEnergyLimit(10*eV),
92 intrinsicHighEnergyLimit(100*GeV) 94 intrinsicHighEnergyLimit(100*GeV)
93 { 95 {
94 if (lowEnergyLimit < intrinsicLowEnergyLimit 96 if (lowEnergyLimit < intrinsicLowEnergyLimit ||
95 highEnergyLimit > intrinsicHighEnergyLim 97 highEnergyLimit > intrinsicHighEnergyLimit)
96 { 98 {
97 G4Exception("G4LowEnergyPolarizedCompton 99 G4Exception("G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton()",
98 "OutOfRange", FatalException 100 "OutOfRange", FatalException,
99 "Energy outside intrinsic pr 101 "Energy outside intrinsic process validity range!");
100 } 102 }
101 103
102 crossSectionHandler = new G4RDCrossSectionHa 104 crossSectionHandler = new G4RDCrossSectionHandler;
103 105
104 106
105 G4RDVDataSetAlgorithm* scatterInterpolation 107 G4RDVDataSetAlgorithm* scatterInterpolation = new G4RDLogLogInterpolation;
106 G4String scatterFile = "comp/ce-sf-"; 108 G4String scatterFile = "comp/ce-sf-";
107 scatterFunctionData = new G4RDCompositeEMDat 109 scatterFunctionData = new G4RDCompositeEMDataSet(scatterInterpolation,1.,1.);
108 scatterFunctionData->LoadData(scatterFile); 110 scatterFunctionData->LoadData(scatterFile);
109 111
110 meanFreePathTable = 0; 112 meanFreePathTable = 0;
111 113
112 rangeTest = new G4RDRangeTest; 114 rangeTest = new G4RDRangeTest;
113 115
114 // For Doppler broadening 116 // For Doppler broadening
115 shellData.SetOccupancyData(); 117 shellData.SetOccupancyData();
116 118
117 119
118 if (verboseLevel > 0) 120 if (verboseLevel > 0)
119 { 121 {
120 G4cout << GetProcessName() << " is crea 122 G4cout << GetProcessName() << " is created " << G4endl
121 << "Energy range: " 123 << "Energy range: "
122 << lowEnergyLimit / keV << " keV - " 124 << lowEnergyLimit / keV << " keV - "
123 << highEnergyLimit / GeV << " GeV" 125 << highEnergyLimit / GeV << " GeV"
124 << G4endl; 126 << G4endl;
125 } 127 }
126 } 128 }
127 129
128 // destructor 130 // destructor
129 131
130 G4LowEnergyPolarizedCompton::~G4LowEnergyPolar 132 G4LowEnergyPolarizedCompton::~G4LowEnergyPolarizedCompton()
131 { 133 {
132 delete meanFreePathTable; 134 delete meanFreePathTable;
133 delete crossSectionHandler; 135 delete crossSectionHandler;
134 delete scatterFunctionData; 136 delete scatterFunctionData;
135 delete rangeTest; 137 delete rangeTest;
136 } 138 }
137 139
138 140
139 void G4LowEnergyPolarizedCompton::BuildPhysics 141 void G4LowEnergyPolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& )
140 { 142 {
141 143
142 crossSectionHandler->Clear(); 144 crossSectionHandler->Clear();
143 G4String crossSectionFile = "comp/ce-cs-"; 145 G4String crossSectionFile = "comp/ce-cs-";
144 crossSectionHandler->LoadData(crossSectionFi 146 crossSectionHandler->LoadData(crossSectionFile);
145 delete meanFreePathTable; 147 delete meanFreePathTable;
146 meanFreePathTable = crossSectionHandler->Bui 148 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
147 149
148 // For Doppler broadening 150 // For Doppler broadening
149 G4String file = "/doppler/shell-doppler"; 151 G4String file = "/doppler/shell-doppler";
150 shellData.LoadData(file); 152 shellData.LoadData(file);
151 153
152 } 154 }
153 155
154 G4VParticleChange* G4LowEnergyPolarizedCompton 156 G4VParticleChange* G4LowEnergyPolarizedCompton::PostStepDoIt(const G4Track& aTrack,
155 const G4Step& aStep) 157 const G4Step& aStep)
156 { 158 {
157 // The scattered gamma energy is sampled acc 159 // The scattered gamma energy is sampled according to Klein - Nishina formula.
158 // The random number techniques of Butcher & 160 // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
159 // GEANT4 internal units 161 // GEANT4 internal units
160 // 162 //
161 // Note : Effects due to binding of atomic e 163 // Note : Effects due to binding of atomic electrons are negliged.
162 164
163 aParticleChange.Initialize(aTrack); 165 aParticleChange.Initialize(aTrack);
164 166
165 // Dynamic particle quantities 167 // Dynamic particle quantities
166 const G4DynamicParticle* incidentPhoton = aT 168 const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
167 G4double gammaEnergy0 = incidentPhoton->GetK 169 G4double gammaEnergy0 = incidentPhoton->GetKineticEnergy();
168 G4ThreeVector gammaPolarization0 = incidentP 170 G4ThreeVector gammaPolarization0 = incidentPhoton->GetPolarization();
169 171
170 // gammaPolarization0 = gammaPolarization0. 172 // gammaPolarization0 = gammaPolarization0.unit(); //
171 173
172 // Protection: a polarisation parallel to th 174 // Protection: a polarisation parallel to the
173 // direction causes problems; 175 // direction causes problems;
174 // in that case find a random polarization 176 // in that case find a random polarization
175 177
176 G4ThreeVector gammaDirection0 = incidentPhot 178 G4ThreeVector gammaDirection0 = incidentPhoton->GetMomentumDirection();
177 // ---- MGP ---- Next two lines commented ou 179 // ---- MGP ---- Next two lines commented out to remove compilation warnings
178 // G4double scalarproduct = gammaPolarizatio 180 // G4double scalarproduct = gammaPolarization0.dot(gammaDirection0);
179 // G4double angle = gammaPolarization0.angle 181 // G4double angle = gammaPolarization0.angle(gammaDirection0);
180 182
181 // Make sure that the polarization vector is 183 // Make sure that the polarization vector is perpendicular to the
182 // gamma direction. If not 184 // gamma direction. If not
183 185
184 if(!(gammaPolarization0.isOrthogonal(gammaDi 186 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
185 { // only for testing now 187 { // only for testing now
186 gammaPolarization0 = GetRandomPolarizati 188 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
187 } 189 }
188 else 190 else
189 { 191 {
190 if ( gammaPolarization0.howOrthogonal(ga 192 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
191 { 193 {
192 gammaPolarization0 = GetPerpendicularPolar 194 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
193 } 195 }
194 } 196 }
195 197
196 // End of Protection 198 // End of Protection
197 199
198 // Within energy limit? 200 // Within energy limit?
199 201
200 if(gammaEnergy0 <= lowEnergyLimit) 202 if(gammaEnergy0 <= lowEnergyLimit)
201 { 203 {
202 aParticleChange.ProposeTrackStatus(fStop 204 aParticleChange.ProposeTrackStatus(fStopAndKill);
203 aParticleChange.ProposeEnergy(0.); 205 aParticleChange.ProposeEnergy(0.);
204 aParticleChange.ProposeLocalEnergyDeposi 206 aParticleChange.ProposeLocalEnergyDeposit(gammaEnergy0);
205 return G4VDiscreteProcess::PostStepDoIt( 207 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
206 } 208 }
207 209
208 G4double E0_m = gammaEnergy0 / electron_mass 210 G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
209 211
210 // Select randomly one element in the curren 212 // Select randomly one element in the current material
211 213
212 const G4MaterialCutsCouple* couple = aTrack. 214 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
213 G4int Z = crossSectionHandler->SelectRandomA 215 G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
214 216
215 // Sample the energy and the polarization of 217 // Sample the energy and the polarization of the scattered photon
216 218
217 G4double epsilon, epsilonSq, onecost, sinThe 219 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
218 220
219 G4double epsilon0 = 1./(1. + 2*E0_m); 221 G4double epsilon0 = 1./(1. + 2*E0_m);
220 G4double epsilon0Sq = epsilon0*epsilon0; 222 G4double epsilon0Sq = epsilon0*epsilon0;
221 G4double alpha1 = - std::log(epsilon0); 223 G4double alpha1 = - std::log(epsilon0);
222 G4double alpha2 = 0.5*(1.- epsilon0Sq); 224 G4double alpha2 = 0.5*(1.- epsilon0Sq);
223 225
224 G4double wlGamma = h_Planck*c_light/gammaEne 226 G4double wlGamma = h_Planck*c_light/gammaEnergy0;
225 G4double gammaEnergy1; 227 G4double gammaEnergy1;
226 G4ThreeVector gammaDirection1; 228 G4ThreeVector gammaDirection1;
227 229
228 do { 230 do {
229 if ( alpha1/(alpha1+alpha2) > G4UniformRan 231 if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
230 { 232 {
231 epsilon = std::exp(-alpha1*G4UniformRand() 233 epsilon = std::exp(-alpha1*G4UniformRand());
232 epsilonSq = epsilon*epsilon; 234 epsilonSq = epsilon*epsilon;
233 } 235 }
234 else 236 else
235 { 237 {
236 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4 238 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
237 epsilon = std::sqrt(epsilonSq); 239 epsilon = std::sqrt(epsilonSq);
238 } 240 }
239 241
240 onecost = (1.- epsilon)/(epsilon*E0_m); 242 onecost = (1.- epsilon)/(epsilon*E0_m);
241 sinThetaSqr = onecost*(2.-onecost); 243 sinThetaSqr = onecost*(2.-onecost);
242 244
243 // Protection 245 // Protection
244 if (sinThetaSqr > 1.) 246 if (sinThetaSqr > 1.)
245 { 247 {
246 if (verboseLevel>0) G4cout 248 if (verboseLevel>0) G4cout
247 << " -- Warning -- G4LowEnergyPolarizedCom 249 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
248 << "sin(theta)**2 = " 250 << "sin(theta)**2 = "
249 << sinThetaSqr 251 << sinThetaSqr
250 << "; set to 1" 252 << "; set to 1"
251 << G4endl; 253 << G4endl;
252 sinThetaSqr = 1.; 254 sinThetaSqr = 1.;
253 } 255 }
254 if (sinThetaSqr < 0.) 256 if (sinThetaSqr < 0.)
255 { 257 {
256 if (verboseLevel>0) G4cout 258 if (verboseLevel>0) G4cout
257 << " -- Warning -- G4LowEnergyPolarizedCom 259 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
258 << "sin(theta)**2 = " 260 << "sin(theta)**2 = "
259 << sinThetaSqr 261 << sinThetaSqr
260 << "; set to 0" 262 << "; set to 0"
261 << G4endl; 263 << G4endl;
262 sinThetaSqr = 0.; 264 sinThetaSqr = 0.;
263 } 265 }
264 // End protection 266 // End protection
265 267
266 G4double x = std::sqrt(onecost/2.) / (wlG 268 G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);;
267 G4double scatteringFunction = scatterFunct 269 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
268 greject = (1. - epsilon*sinThetaSqr/(1.+ e 270 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
269 271
270 } while(greject < G4UniformRand()*Z); 272 } while(greject < G4UniformRand()*Z);
271 273
272 274
273 // ***************************************** 275 // ****************************************************
274 // Phi determination 276 // Phi determination
275 // ***************************************** 277 // ****************************************************
276 278
277 G4double phi = SetPhi(epsilon,sinThetaSqr); 279 G4double phi = SetPhi(epsilon,sinThetaSqr);
278 280
279 // 281 //
280 // scattered gamma angles. ( Z - axis along 282 // scattered gamma angles. ( Z - axis along the parent gamma)
281 // 283 //
282 284
283 G4double cosTheta = 1. - onecost; 285 G4double cosTheta = 1. - onecost;
284 286
285 // Protection 287 // Protection
286 288
287 if (cosTheta > 1.) 289 if (cosTheta > 1.)
288 { 290 {
289 if (verboseLevel>0) G4cout 291 if (verboseLevel>0) G4cout
290 << " -- Warning -- G4LowEnergyPolarizedCompt 292 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
291 << "cosTheta = " 293 << "cosTheta = "
292 << cosTheta 294 << cosTheta
293 << "; set to 1" 295 << "; set to 1"
294 << G4endl; 296 << G4endl;
295 cosTheta = 1.; 297 cosTheta = 1.;
296 } 298 }
297 if (cosTheta < -1.) 299 if (cosTheta < -1.)
298 { 300 {
299 if (verboseLevel>0) G4cout 301 if (verboseLevel>0) G4cout
300 << " -- Warning -- G4LowEnergyPolarizedCompt 302 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
301 << "cosTheta = " 303 << "cosTheta = "
302 << cosTheta 304 << cosTheta
303 << "; set to -1" 305 << "; set to -1"
304 << G4endl; 306 << G4endl;
305 cosTheta = -1.; 307 cosTheta = -1.;
306 } 308 }
307 // End protection 309 // End protection
308 310
309 311
310 G4double sinTheta = std::sqrt (sinThetaSqr); 312 G4double sinTheta = std::sqrt (sinThetaSqr);
311 313
312 // Protection 314 // Protection
313 if (sinTheta > 1.) 315 if (sinTheta > 1.)
314 { 316 {
315 if (verboseLevel>0) G4cout 317 if (verboseLevel>0) G4cout
316 << " -- Warning -- G4LowEnergyPolarizedCompt 318 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
317 << "sinTheta = " 319 << "sinTheta = "
318 << sinTheta 320 << sinTheta
319 << "; set to 1" 321 << "; set to 1"
320 << G4endl; 322 << G4endl;
321 sinTheta = 1.; 323 sinTheta = 1.;
322 } 324 }
323 if (sinTheta < -1.) 325 if (sinTheta < -1.)
324 { 326 {
325 if (verboseLevel>0) G4cout 327 if (verboseLevel>0) G4cout
326 << " -- Warning -- G4LowEnergyPolarizedCompt 328 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
327 << "sinTheta = " 329 << "sinTheta = "
328 << sinTheta 330 << sinTheta
329 << "; set to -1" 331 << "; set to -1"
330 << G4endl; 332 << G4endl;
331 sinTheta = -1.; 333 sinTheta = -1.;
332 } 334 }
333 // End protection 335 // End protection
334 336
335 337
336 G4double dirx = sinTheta*std::cos(phi); 338 G4double dirx = sinTheta*std::cos(phi);
337 G4double diry = sinTheta*std::sin(phi); 339 G4double diry = sinTheta*std::sin(phi);
338 G4double dirz = cosTheta ; 340 G4double dirz = cosTheta ;
339 341
340 342
341 // oneCosT , eom 343 // oneCosT , eom
342 344
343 345
344 346
345 // Doppler broadening - Method based on: 347 // Doppler broadening - Method based on:
346 // Y. Namito, S. Ban and H. Hirayama, 348 // Y. Namito, S. Ban and H. Hirayama,
347 // "Implementation of the Doppler Broadening 349 // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
348 // NIM A 349, pp. 489-494, 1994 350 // NIM A 349, pp. 489-494, 1994
349 351
350 // Maximum number of sampling iterations 352 // Maximum number of sampling iterations
351 353
352 G4int maxDopplerIterations = 1000; 354 G4int maxDopplerIterations = 1000;
353 G4double bindingE = 0.; 355 G4double bindingE = 0.;
354 G4double photonEoriginal = epsilon * gammaEn 356 G4double photonEoriginal = epsilon * gammaEnergy0;
355 G4double photonE = -1.; 357 G4double photonE = -1.;
356 G4int iteration = 0; 358 G4int iteration = 0;
357 G4double eMax = gammaEnergy0; 359 G4double eMax = gammaEnergy0;
358 360
359 do 361 do
360 { 362 {
361 iteration++; 363 iteration++;
362 // Select shell based on shell occupancy 364 // Select shell based on shell occupancy
363 G4int shell = shellData.SelectRandomShel 365 G4int shell = shellData.SelectRandomShell(Z);
364 bindingE = shellData.BindingEnergy(Z,she 366 bindingE = shellData.BindingEnergy(Z,shell);
365 367
366 eMax = gammaEnergy0 - bindingE; 368 eMax = gammaEnergy0 - bindingE;
367 369
368 // Randomly sample bound electron moment 370 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
369 G4double pSample = profileData.RandomSel 371 G4double pSample = profileData.RandomSelectMomentum(Z,shell);
370 // Rescale from atomic units 372 // Rescale from atomic units
371 G4double pDoppler = pSample * fine_struc 373 G4double pDoppler = pSample * fine_structure_const;
372 G4double pDoppler2 = pDoppler * pDoppler 374 G4double pDoppler2 = pDoppler * pDoppler;
373 G4double var2 = 1. + onecost * E0_m; 375 G4double var2 = 1. + onecost * E0_m;
374 G4double var3 = var2*var2 - pDoppler2; 376 G4double var3 = var2*var2 - pDoppler2;
375 G4double var4 = var2 - pDoppler2 * cosTh 377 G4double var4 = var2 - pDoppler2 * cosTheta;
376 G4double var = var4*var4 - var3 + pDoppl 378 G4double var = var4*var4 - var3 + pDoppler2 * var3;
377 if (var > 0.) 379 if (var > 0.)
378 { 380 {
379 G4double varSqrt = std::sqrt(var); 381 G4double varSqrt = std::sqrt(var);
380 G4double scale = gammaEnergy0 / var3; 382 G4double scale = gammaEnergy0 / var3;
381 // Random select either root 383 // Random select either root
382 if (G4UniformRand() < 0.5) photonE = (var4 384 if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
383 else photonE = (var4 + varSqrt) * scale; 385 else photonE = (var4 + varSqrt) * scale;
384 } 386 }
385 else 387 else
386 { 388 {
387 photonE = -1.; 389 photonE = -1.;
388 } 390 }
389 } while ( iteration <= maxDopplerIterations 391 } while ( iteration <= maxDopplerIterations &&
390 (photonE < 0. || photonE > eMax || phot 392 (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
391 393
392 // End of recalculation of photon energy wit 394 // End of recalculation of photon energy with Doppler broadening
393 // Revert to original if maximum number of i 395 // Revert to original if maximum number of iterations threshold has been reached
394 if (iteration >= maxDopplerIterations) 396 if (iteration >= maxDopplerIterations)
395 { 397 {
396 photonE = photonEoriginal; 398 photonE = photonEoriginal;
397 bindingE = 0.; 399 bindingE = 0.;
398 } 400 }
399 401
400 gammaEnergy1 = photonE; 402 gammaEnergy1 = photonE;
401 403
402 // G4cout << "--> PHOTONENERGY1 = " << photo 404 // G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl;
403 405
404 406
405 /// Doppler Broadeing 407 /// Doppler Broadeing
406 408
407 409
408 410
409 411
410 // 412 //
411 // update G4VParticleChange for the scattere 413 // update G4VParticleChange for the scattered photon
412 // 414 //
413 415
414 // gammaEnergy1 = epsilon*gammaEnergy0; 416 // gammaEnergy1 = epsilon*gammaEnergy0;
415 417
416 418
417 // New polarization 419 // New polarization
418 420
419 G4ThreeVector gammaPolarization1 = SetNewPol 421 G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
420 sinThetaSqr, 422 sinThetaSqr,
421 phi, 423 phi,
422 cosTheta); 424 cosTheta);
423 425
424 // Set new direction 426 // Set new direction
425 G4ThreeVector tmpDirection1( dirx,diry,dirz 427 G4ThreeVector tmpDirection1( dirx,diry,dirz );
426 gammaDirection1 = tmpDirection1; 428 gammaDirection1 = tmpDirection1;
427 429
428 // Change reference frame. 430 // Change reference frame.
429 431
430 SystemOfRefChange(gammaDirection0,gammaDirec 432 SystemOfRefChange(gammaDirection0,gammaDirection1,
431 gammaPolarization0,gammaPolarization1) 433 gammaPolarization0,gammaPolarization1);
432 434
433 if (gammaEnergy1 > 0.) 435 if (gammaEnergy1 > 0.)
434 { 436 {
435 aParticleChange.ProposeEnergy( gammaEner 437 aParticleChange.ProposeEnergy( gammaEnergy1 ) ;
436 aParticleChange.ProposeMomentumDirection 438 aParticleChange.ProposeMomentumDirection( gammaDirection1 );
437 aParticleChange.ProposePolarization( gam 439 aParticleChange.ProposePolarization( gammaPolarization1 );
438 } 440 }
439 else 441 else
440 { 442 {
441 aParticleChange.ProposeEnergy(0.) ; 443 aParticleChange.ProposeEnergy(0.) ;
442 aParticleChange.ProposeTrackStatus(fStop 444 aParticleChange.ProposeTrackStatus(fStopAndKill);
443 } 445 }
444 446
445 // 447 //
446 // kinematic of the scattered electron 448 // kinematic of the scattered electron
447 // 449 //
448 450
449 G4double ElecKineEnergy = gammaEnergy0 - gam 451 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
450 452
451 453
452 // Generate the electron only if with large 454 // Generate the electron only if with large enough range w.r.t. cuts and safety
453 455
454 G4double safety = aStep.GetPostStepPoint()-> 456 G4double safety = aStep.GetPostStepPoint()->GetSafety();
455 457
456 458
457 if (rangeTest->Escape(G4Electron::Electron() 459 if (rangeTest->Escape(G4Electron::Electron(),couple,ElecKineEnergy,safety))
458 { 460 {
459 G4double ElecMomentum = std::sqrt(ElecKi 461 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
460 G4ThreeVector ElecDirection((gammaEnergy 462 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
461 gammaEnergy1 * gammaDirection1) * ( 463 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
462 G4DynamicParticle* electron = new G4Dyna 464 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
463 aParticleChange.SetNumberOfSecondaries(1 465 aParticleChange.SetNumberOfSecondaries(1);
464 aParticleChange.AddSecondary(electron); 466 aParticleChange.AddSecondary(electron);
465 // aParticleChange.ProposeLocalEner 467 // aParticleChange.ProposeLocalEnergyDeposit(0.);
466 aParticleChange.ProposeLocalEnergyDeposi 468 aParticleChange.ProposeLocalEnergyDeposit(bindingE);
467 } 469 }
468 else 470 else
469 { 471 {
470 aParticleChange.SetNumberOfSecondaries(0 472 aParticleChange.SetNumberOfSecondaries(0);
471 aParticleChange.ProposeLocalEnergyDeposi 473 aParticleChange.ProposeLocalEnergyDeposit(ElecKineEnergy+bindingE);
472 } 474 }
473 475
474 return G4VDiscreteProcess::PostStepDoIt( aTr 476 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
475 477
476 } 478 }
477 479
478 480
479 G4double G4LowEnergyPolarizedCompton::SetPhi(G 481 G4double G4LowEnergyPolarizedCompton::SetPhi(G4double energyRate,
480 G4double sinSqrTh) 482 G4double sinSqrTh)
481 { 483 {
482 G4double rand1; 484 G4double rand1;
483 G4double rand2; 485 G4double rand2;
484 G4double phiProbability; 486 G4double phiProbability;
485 G4double phi; 487 G4double phi;
486 G4double a, b; 488 G4double a, b;
487 489
488 do 490 do
489 { 491 {
490 rand1 = G4UniformRand(); 492 rand1 = G4UniformRand();
491 rand2 = G4UniformRand(); 493 rand2 = G4UniformRand();
492 phiProbability=0.; 494 phiProbability=0.;
493 phi = twopi*rand1; 495 phi = twopi*rand1;
494 496
495 a = 2*sinSqrTh; 497 a = 2*sinSqrTh;
496 b = energyRate + 1/energyRate; 498 b = energyRate + 1/energyRate;
497 499
498 phiProbability = 1 - (a/b)*(std::cos(phi 500 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
499 501
500 502
501 503
502 } 504 }
503 while ( rand2 > phiProbability ); 505 while ( rand2 > phiProbability );
504 return phi; 506 return phi;
505 } 507 }
506 508
507 509
508 G4ThreeVector G4LowEnergyPolarizedCompton::Set 510 G4ThreeVector G4LowEnergyPolarizedCompton::SetPerpendicularVector(G4ThreeVector& a)
509 { 511 {
510 G4double dx = a.x(); 512 G4double dx = a.x();
511 G4double dy = a.y(); 513 G4double dy = a.y();
512 G4double dz = a.z(); 514 G4double dz = a.z();
513 G4double x = dx < 0.0 ? -dx : dx; 515 G4double x = dx < 0.0 ? -dx : dx;
514 G4double y = dy < 0.0 ? -dy : dy; 516 G4double y = dy < 0.0 ? -dy : dy;
515 G4double z = dz < 0.0 ? -dz : dz; 517 G4double z = dz < 0.0 ? -dz : dz;
516 if (x < y) { 518 if (x < y) {
517 return x < z ? G4ThreeVector(-dy,dx,0) : G 519 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
518 }else{ 520 }else{
519 return y < z ? G4ThreeVector(dz,0,-dx) : G 521 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
520 } 522 }
521 } 523 }
522 524
523 G4ThreeVector G4LowEnergyPolarizedCompton::Get 525 G4ThreeVector G4LowEnergyPolarizedCompton::GetRandomPolarization(G4ThreeVector& direction0)
524 { 526 {
525 G4ThreeVector d0 = direction0.unit(); 527 G4ThreeVector d0 = direction0.unit();
526 G4ThreeVector a1 = SetPerpendicularVector(d0 528 G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
527 G4ThreeVector a0 = a1.unit(); // unit vector 529 G4ThreeVector a0 = a1.unit(); // unit vector
528 530
529 G4double rand1 = G4UniformRand(); 531 G4double rand1 = G4UniformRand();
530 532
531 G4double angle = twopi*rand1; // random pola 533 G4double angle = twopi*rand1; // random polar angle
532 G4ThreeVector b0 = d0.cross(a0); // cross pr 534 G4ThreeVector b0 = d0.cross(a0); // cross product
533 535
534 G4ThreeVector c; 536 G4ThreeVector c;
535 537
536 c.setX(std::cos(angle)*(a0.x())+std::sin(ang 538 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
537 c.setY(std::cos(angle)*(a0.y())+std::sin(ang 539 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
538 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang 540 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
539 541
540 G4ThreeVector c0 = c.unit(); 542 G4ThreeVector c0 = c.unit();
541 543
542 return c0; 544 return c0;
543 545
544 } 546 }
545 547
546 548
547 G4ThreeVector G4LowEnergyPolarizedCompton::Get 549 G4ThreeVector G4LowEnergyPolarizedCompton::GetPerpendicularPolarization
548 (const G4ThreeVector& gammaDirection, const G4 550 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
549 { 551 {
550 552
551 // 553 //
552 // The polarization of a photon is always pe 554 // The polarization of a photon is always perpendicular to its momentum direction.
553 // Therefore this function removes those vec 555 // Therefore this function removes those vector component of gammaPolarization, which
554 // points in direction of gammaDirection 556 // points in direction of gammaDirection
555 // 557 //
556 // Mathematically we search the projection o 558 // Mathematically we search the projection of the vector a on the plane E, where n is the
557 // plains normal vector. 559 // plains normal vector.
558 // The basic equation can be found in each g 560 // The basic equation can be found in each geometry book (e.g. Bronstein):
559 // p = a - (a o n)/(n o n)*n 561 // p = a - (a o n)/(n o n)*n
560 562
561 return gammaPolarization - gammaPolarization 563 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
562 } 564 }
563 565
564 566
565 G4ThreeVector G4LowEnergyPolarizedCompton::Set 567 G4ThreeVector G4LowEnergyPolarizedCompton::SetNewPolarization(G4double epsilon,
566 G4double sinSqrTh, 568 G4double sinSqrTh,
567 G4double phi, 569 G4double phi,
568 G4double costheta) 570 G4double costheta)
569 { 571 {
570 G4double rand1; 572 G4double rand1;
571 G4double rand2; 573 G4double rand2;
572 G4double cosPhi = std::cos(phi); 574 G4double cosPhi = std::cos(phi);
573 G4double sinPhi = std::sin(phi); 575 G4double sinPhi = std::sin(phi);
574 G4double sinTheta = std::sqrt(sinSqrTh); 576 G4double sinTheta = std::sqrt(sinSqrTh);
575 G4double cosSqrPhi = cosPhi*cosPhi; 577 G4double cosSqrPhi = cosPhi*cosPhi;
576 // G4double cossqrth = 1.-sinSqrTh; 578 // G4double cossqrth = 1.-sinSqrTh;
577 // G4double sinsqrphi = sinPhi*sinPhi; 579 // G4double sinsqrphi = sinPhi*sinPhi;
578 G4double normalisation = std::sqrt(1. - cosS 580 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
579 581
580 582
581 // Determination of Theta 583 // Determination of Theta
582 584
583 // ---- MGP ---- Commented out the following 585 // ---- MGP ---- Commented out the following 3 lines to avoid compilation
584 // warnings (unused variables) 586 // warnings (unused variables)
585 // G4double thetaProbability; 587 // G4double thetaProbability;
586 G4double theta; 588 G4double theta;
587 // G4double a, b; 589 // G4double a, b;
588 // G4double cosTheta; 590 // G4double cosTheta;
589 591
590 /* 592 /*
591 593
592 depaola method 594 depaola method
593 595
594 do 596 do
595 { 597 {
596 rand1 = G4UniformRand(); 598 rand1 = G4UniformRand();
597 rand2 = G4UniformRand(); 599 rand2 = G4UniformRand();
598 thetaProbability=0.; 600 thetaProbability=0.;
599 theta = twopi*rand1; 601 theta = twopi*rand1;
600 a = 4*normalisation*normalisation; 602 a = 4*normalisation*normalisation;
601 b = (epsilon + 1/epsilon) - 2; 603 b = (epsilon + 1/epsilon) - 2;
602 thetaProbability = (b + a*std::cos(theta 604 thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
603 cosTheta = std::cos(theta); 605 cosTheta = std::cos(theta);
604 } 606 }
605 while ( rand2 > thetaProbability ); 607 while ( rand2 > thetaProbability );
606 608
607 G4double cosBeta = cosTheta; 609 G4double cosBeta = cosTheta;
608 610
609 */ 611 */
610 612
611 613
612 // Dan Xu method (IEEE TNS, 52, 1160 (2005)) 614 // Dan Xu method (IEEE TNS, 52, 1160 (2005))
613 615
614 rand1 = G4UniformRand(); 616 rand1 = G4UniformRand();
615 rand2 = G4UniformRand(); 617 rand2 = G4UniformRand();
616 618
617 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsi 619 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
618 { 620 {
619 if (rand2<0.5) 621 if (rand2<0.5)
620 theta = pi/2.0; 622 theta = pi/2.0;
621 else 623 else
622 theta = 3.0*pi/2.0; 624 theta = 3.0*pi/2.0;
623 } 625 }
624 else 626 else
625 { 627 {
626 if (rand2<0.5) 628 if (rand2<0.5)
627 theta = 0; 629 theta = 0;
628 else 630 else
629 theta = pi; 631 theta = pi;
630 } 632 }
631 G4double cosBeta = std::cos(theta); 633 G4double cosBeta = std::cos(theta);
632 G4double sinBeta = std::sqrt(1-cosBeta*cosBe 634 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
633 635
634 G4ThreeVector gammaPolarization1; 636 G4ThreeVector gammaPolarization1;
635 637
636 G4double xParallel = normalisation*cosBeta; 638 G4double xParallel = normalisation*cosBeta;
637 G4double yParallel = -(sinSqrTh*cosPhi*sinPh 639 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
638 G4double zParallel = -(costheta*sinTheta*cos 640 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
639 G4double xPerpendicular = 0.; 641 G4double xPerpendicular = 0.;
640 G4double yPerpendicular = (costheta)*sinBeta 642 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
641 G4double zPerpendicular = -(sinTheta*sinPhi) 643 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
642 644
643 G4double xTotal = (xParallel + xPerpendicula 645 G4double xTotal = (xParallel + xPerpendicular);
644 G4double yTotal = (yParallel + yPerpendicula 646 G4double yTotal = (yParallel + yPerpendicular);
645 G4double zTotal = (zParallel + zPerpendicula 647 G4double zTotal = (zParallel + zPerpendicular);
646 648
647 gammaPolarization1.setX(xTotal); 649 gammaPolarization1.setX(xTotal);
648 gammaPolarization1.setY(yTotal); 650 gammaPolarization1.setY(yTotal);
649 gammaPolarization1.setZ(zTotal); 651 gammaPolarization1.setZ(zTotal);
650 652
651 return gammaPolarization1; 653 return gammaPolarization1;
652 654
653 } 655 }
654 656
655 657
656 void G4LowEnergyPolarizedCompton::SystemOfRefC 658 void G4LowEnergyPolarizedCompton::SystemOfRefChange(G4ThreeVector& direction0,
657 G4ThreeVector& direction1, 659 G4ThreeVector& direction1,
658 G4ThreeVector& polarization0, 660 G4ThreeVector& polarization0,
659 G4ThreeVector& polarization1) 661 G4ThreeVector& polarization1)
660 { 662 {
661 // direction0 is the original photon directi 663 // direction0 is the original photon direction ---> z
662 // polarization0 is the original photon pola 664 // polarization0 is the original photon polarization ---> x
663 // need to specify y axis in the real refere 665 // need to specify y axis in the real reference frame ---> y
664 G4ThreeVector Axis_Z0 = direction0.unit(); 666 G4ThreeVector Axis_Z0 = direction0.unit();
665 G4ThreeVector Axis_X0 = polarization0.unit() 667 G4ThreeVector Axis_X0 = polarization0.unit();
666 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_ 668 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
667 669
668 G4double direction_x = direction1.getX(); 670 G4double direction_x = direction1.getX();
669 G4double direction_y = direction1.getY(); 671 G4double direction_y = direction1.getY();
670 G4double direction_z = direction1.getZ(); 672 G4double direction_z = direction1.getZ();
671 673
672 direction1 = (direction_x*Axis_X0 + directio 674 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
673 G4double polarization_x = polarization1.getX 675 G4double polarization_x = polarization1.getX();
674 G4double polarization_y = polarization1.getY 676 G4double polarization_y = polarization1.getY();
675 G4double polarization_z = polarization1.getZ 677 G4double polarization_z = polarization1.getZ();
676 678
677 polarization1 = (polarization_x*Axis_X0 + po 679 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
678 680
679 } 681 }
680 682
681 683
682 G4bool G4LowEnergyPolarizedCompton::IsApplicab 684 G4bool G4LowEnergyPolarizedCompton::IsApplicable(const G4ParticleDefinition& particle)
683 { 685 {
684 return ( &particle == G4Gamma::Gamma() ); 686 return ( &particle == G4Gamma::Gamma() );
685 } 687 }
686 688
687 689
688 G4double G4LowEnergyPolarizedCompton::GetMeanF 690 G4double G4LowEnergyPolarizedCompton::GetMeanFreePath(const G4Track& track,
689 G4double, 691 G4double,
690 G4ForceCondition*) 692 G4ForceCondition*)
691 { 693 {
692 const G4DynamicParticle* photon = track.GetD 694 const G4DynamicParticle* photon = track.GetDynamicParticle();
693 G4double energy = photon->GetKineticEnergy() 695 G4double energy = photon->GetKineticEnergy();
694 const G4MaterialCutsCouple* couple = track.G 696 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
695 size_t materialIndex = couple->GetIndex(); 697 size_t materialIndex = couple->GetIndex();
696 G4double meanFreePath; 698 G4double meanFreePath;
697 if (energy > highEnergyLimit) meanFreePath = 699 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
698 else if (energy < lowEnergyLimit) meanFreePa 700 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
699 else meanFreePath = meanFreePathTable->FindV 701 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
700 return meanFreePath; 702 return meanFreePath;
701 } 703 }
702 704
703 705
704 706
705 707
706 708
707 709
708 710
709 711
710 712
711 713
712 714
713 715
714 716
715 717
716 718