| 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 // | 27 // | |
28 // | G4LowEPPolarizedComptonModel-- Geant 28 // | G4LowEPPolarizedComptonModel-- Geant4 Monash University |
29 // | polarised low energy Compton scat 29 // | polarised low energy Compton scattering model. |
30 // | J. M. C. Brown, Monash Universit 30 // | J. M. C. Brown, Monash University, Australia |
31 // | 31 // | |
32 // | 32 // | |
33 // ******************************************* 33 // *********************************************************************
34 // | 34 // | |
35 // | The following is a Geant4 class to simula 35 // | The following is a Geant4 class to simulate the process of |
36 // | bound electron Compton scattering. Genera 36 // | bound electron Compton scattering. General code structure is |
37 // | based on G4LowEnergyCompton.cc and 37 // | based on G4LowEnergyCompton.cc and |
38 // | G4LivermorePolarizedComptonModel.cc. 38 // | G4LivermorePolarizedComptonModel.cc. |
39 // | Algorithms for photon energy, and ejected 39 // | Algorithms for photon energy, and ejected Compton electron |
40 // | direction taken from: 40 // | direction taken from: |
41 // | 41 // | |
42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gill 42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin, |
43 // | "A low energy bound atomic electron Compt 43 // | "A low energy bound atomic electron Compton scattering model |
44 // | for Geant4", NIMB, Vol. 338, 77-88, 2014 44 // | for Geant4", NIMB, Vol. 338, 77-88, 2014. |
45 // | 45 // | |
46 // | The author acknowledges the work of the G 46 // | The author acknowledges the work of the Geant4 collaboration |
47 // | in developing the following algorithms th 47 // | in developing the following algorithms that have been employed |
48 // | or adapeted for the present software: 48 // | or adapeted for the present software: |
49 // | 49 // | |
50 // | # sampling of photon scattering angle, 50 // | # sampling of photon scattering angle, |
51 // | # target element selection in composite 51 // | # target element selection in composite materials, |
52 // | # target shell selection in element, 52 // | # target shell selection in element, |
53 // | # and sampling of bound electron momentu 53 // | # and sampling of bound electron momentum from Compton profiles. |
54 // | 54 // | |
55 // ******************************************* 55 // *********************************************************************
56 // | 56 // | |
57 // | History: 57 // | History: |
58 // | -------- 58 // | -------- |
59 // | 59 // | |
60 // | Jan. 2015 JMCB - 1st Version based 60 // | Jan. 2015 JMCB - 1st Version based on G4LowEPPComptonModel |
61 // | Feb. 2016 JMCB - Geant4 10.2 FPE fi 61 // | Feb. 2016 JMCB - Geant4 10.2 FPE fix for bug 1676 |
62 // | Nov. 2016 JMCB - Polarisation track 62 // | Nov. 2016 JMCB - Polarisation tracking fix in collaboration |
63 // | of Dr. Merlin Reyn 63 // | of Dr. Merlin Reynaard Kole, |
64 // | University of Gene 64 // | University of Geneva |
65 // | 65 // | |
66 // ******************************************* 66 // *********************************************************************
67 67
68 #include "G4LowEPPolarizedComptonModel.hh" 68 #include "G4LowEPPolarizedComptonModel.hh"
69 #include "G4AutoLock.hh" <<
70 #include "G4PhysicalConstants.hh" 69 #include "G4PhysicalConstants.hh"
71 #include "G4SystemOfUnits.hh" 70 #include "G4SystemOfUnits.hh"
72 71
73 //******************************************** 72 //****************************************************************************
74 73
75 using namespace std; 74 using namespace std;
76 namespace { G4Mutex LowEPPolarizedComptonModel <<
77 <<
78 75
79 G4PhysicsFreeVector* G4LowEPPolarizedComptonMo << 76 G4int G4LowEPPolarizedComptonModel::maxZ = 99;
80 G4ShellData* G4LowEPPolarizedComptonMode << 77 G4LPhysicsFreeVector* G4LowEPPolarizedComptonModel::data[] = {0};
81 G4DopplerProfile* G4LowEPPolarizedComptonMode << 78 G4ShellData* G4LowEPPolarizedComptonModel::shellData = 0;
>> 79 G4DopplerProfile* G4LowEPPolarizedComptonModel::profileData = 0;
82 80
83 static const G4double ln10 = G4Log(10.); 81 static const G4double ln10 = G4Log(10.);
84 82
85 G4LowEPPolarizedComptonModel::G4LowEPPolarized 83 G4LowEPPolarizedComptonModel::G4LowEPPolarizedComptonModel(const G4ParticleDefinition*,
86 84 const G4String& nam)
87 : G4VEmModel(nam),isInitialised(false) 85 : G4VEmModel(nam),isInitialised(false)
88 { 86 {
89 verboseLevel=1 ; 87 verboseLevel=1 ;
90 // Verbosity scale: 88 // Verbosity scale:
91 // 0 = nothing 89 // 0 = nothing
92 // 1 = warning for energy non-conservation 90 // 1 = warning for energy non-conservation
93 // 2 = details of energy budget 91 // 2 = details of energy budget
94 // 3 = calculation of cross sections, file o 92 // 3 = calculation of cross sections, file openings, sampling of atoms
95 // 4 = entering in methods 93 // 4 = entering in methods
96 94
97 if( verboseLevel>1 ) { 95 if( verboseLevel>1 ) {
98 G4cout << "Low energy photon Compton model 96 G4cout << "Low energy photon Compton model is constructed " << G4endl;
99 } 97 }
100 98
101 //Mark this model as "applicable" for atomic 99 //Mark this model as "applicable" for atomic deexcitation
102 SetDeexcitationFlag(true); 100 SetDeexcitationFlag(true);
103 101
104 fParticleChange = nullptr; << 102 fParticleChange = 0;
105 fAtomDeexcitation = nullptr; << 103 fAtomDeexcitation = 0;
106 } 104 }
107 105
108 //******************************************** 106 //****************************************************************************
109 107
110 G4LowEPPolarizedComptonModel::~G4LowEPPolarize 108 G4LowEPPolarizedComptonModel::~G4LowEPPolarizedComptonModel()
111 { 109 {
112 if(IsMaster()) { 110 if(IsMaster()) {
113 delete shellData; 111 delete shellData;
114 shellData = nullptr; << 112 shellData = 0;
115 delete profileData; 113 delete profileData;
116 profileData = nullptr; << 114 profileData = 0;
117 } 115 }
118 } 116 }
119 117
120 //******************************************** 118 //****************************************************************************
121 119
122 void G4LowEPPolarizedComptonModel::Initialise( 120 void G4LowEPPolarizedComptonModel::Initialise(const G4ParticleDefinition* particle,
123 const 121 const G4DataVector& cuts)
124 { 122 {
125 if (verboseLevel > 1) { 123 if (verboseLevel > 1) {
126 G4cout << "Calling G4LowEPPolarizedCompton 124 G4cout << "Calling G4LowEPPolarizedComptonModel::Initialise()" << G4endl;
127 } 125 }
128 126
129 // Initialise element selector 127 // Initialise element selector
>> 128
130 if(IsMaster()) { 129 if(IsMaster()) {
>> 130
131 // Access to elements 131 // Access to elements
132 const char* path = G4FindDataDir("G4LEDATA << 132
>> 133 char* path = getenv("G4LEDATA");
133 134
134 G4ProductionCutsTable* theCoupleTable = 135 G4ProductionCutsTable* theCoupleTable =
135 G4ProductionCutsTable::GetProductionCuts 136 G4ProductionCutsTable::GetProductionCutsTable();
136 G4int numOfCouples = (G4int)theCoupleTable << 137 G4int numOfCouples = theCoupleTable->GetTableSize();
137 138
138 for(G4int i=0; i<numOfCouples; ++i) { 139 for(G4int i=0; i<numOfCouples; ++i) {
139 const G4Material* material = 140 const G4Material* material =
140 theCoupleTable->GetMaterialCutsCouple( 141 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
141 const G4ElementVector* theElementVector 142 const G4ElementVector* theElementVector = material->GetElementVector();
142 std::size_t nelm = material->GetNumberOf << 143 G4int nelm = material->GetNumberOfElements();
143 144
144 for (std::size_t j=0; j<nelm; ++j) { << 145 for (G4int j=0; j<nelm; ++j) {
145 G4int Z = G4lrint((*theElementVector)[ 146 G4int Z = G4lrint((*theElementVector)[j]->GetZ());
146 if(Z < 1) { Z = 1; } 147 if(Z < 1) { Z = 1; }
147 else if(Z > maxZ){ Z = maxZ; } 148 else if(Z > maxZ){ Z = maxZ; }
148 149
149 if( (!data[Z]) ) { ReadData(Z, path); 150 if( (!data[Z]) ) { ReadData(Z, path); }
150 } 151 }
151 } 152 }
152 153
153 // For Doppler broadening 154 // For Doppler broadening
154 if(!shellData) { 155 if(!shellData) {
155 shellData = new G4ShellData(); 156 shellData = new G4ShellData();
156 shellData->SetOccupancyData(); 157 shellData->SetOccupancyData();
157 G4String file = "/doppler/shell-doppler" 158 G4String file = "/doppler/shell-doppler";
158 shellData->LoadData(file); 159 shellData->LoadData(file);
159 } 160 }
160 if(!profileData) { profileData = new G4Dop 161 if(!profileData) { profileData = new G4DopplerProfile(); }
161 162
162 InitialiseElementSelectors(particle, cuts) 163 InitialiseElementSelectors(particle, cuts);
163 } 164 }
164 165
165 if (verboseLevel > 2) { 166 if (verboseLevel > 2) {
166 G4cout << "Loaded cross section files" << 167 G4cout << "Loaded cross section files" << G4endl;
167 } 168 }
168 169
169 if( verboseLevel>1 ) { 170 if( verboseLevel>1 ) {
170 G4cout << "G4LowEPPolarizedComptonModel is 171 G4cout << "G4LowEPPolarizedComptonModel is initialized " << G4endl
171 << "Energy range: " 172 << "Energy range: "
172 << LowEnergyLimit() / eV << " eV - 173 << LowEnergyLimit() / eV << " eV - "
173 << HighEnergyLimit() / GeV << " GeV 174 << HighEnergyLimit() / GeV << " GeV"
174 << G4endl; 175 << G4endl;
175 } 176 }
176 177
177 if(isInitialised) { return; } 178 if(isInitialised) { return; }
178 179
179 fParticleChange = GetParticleChangeForGamma( 180 fParticleChange = GetParticleChangeForGamma();
180 fAtomDeexcitation = G4LossTableManager::Ins 181 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
181 isInitialised = true; 182 isInitialised = true;
182 } 183 }
183 184
184 //******************************************** 185 //****************************************************************************
185 186
186 void G4LowEPPolarizedComptonModel::InitialiseL 187 void G4LowEPPolarizedComptonModel::InitialiseLocal(const G4ParticleDefinition*,
187 188 G4VEmModel* masterModel)
188 { 189 {
189 SetElementSelectors(masterModel->GetElementS 190 SetElementSelectors(masterModel->GetElementSelectors());
190 } 191 }
191 192
192 //******************************************** 193 //****************************************************************************
193 194
194 void G4LowEPPolarizedComptonModel::ReadData(st << 195 void G4LowEPPolarizedComptonModel::ReadData(size_t Z, const char* path)
195 { 196 {
196 if (verboseLevel > 1) 197 if (verboseLevel > 1)
197 { 198 {
198 G4cout << "G4LowEPPolarizedComptonModel::R 199 G4cout << "G4LowEPPolarizedComptonModel::ReadData()"
199 << G4endl; 200 << G4endl;
200 } 201 }
201 if(data[Z]) { return; } 202 if(data[Z]) { return; }
202 const char* datadir = path; 203 const char* datadir = path;
203 if(!datadir) 204 if(!datadir)
204 { 205 {
205 datadir = G4FindDataDir("G4LEDATA"); << 206 datadir = getenv("G4LEDATA");
206 if(!datadir) 207 if(!datadir)
207 { 208 {
208 G4Exception("G4LowEPPolarizedComptonMode 209 G4Exception("G4LowEPPolarizedComptonModel::ReadData()",
209 "em0006",FatalException, 210 "em0006",FatalException,
210 "Environment variable G4LEDA 211 "Environment variable G4LEDATA not defined");
211 return; 212 return;
212 } 213 }
213 } 214 }
214 215
215 data[Z] = new G4PhysicsFreeVector(); << 216 data[Z] = new G4LPhysicsFreeVector();
>> 217
>> 218 // Activation of spline interpolation
>> 219 data[Z]->SetSpline(false);
216 220
217 std::ostringstream ost; 221 std::ostringstream ost;
218 ost << datadir << "/livermore/comp/ce-cs-" < 222 ost << datadir << "/livermore/comp/ce-cs-" << Z <<".dat";
219 std::ifstream fin(ost.str().c_str()); 223 std::ifstream fin(ost.str().c_str());
220 224
221 if( !fin.is_open()) 225 if( !fin.is_open())
222 { 226 {
223 G4ExceptionDescription ed; 227 G4ExceptionDescription ed;
224 ed << "G4LowEPPolarizedComptonModel data 228 ed << "G4LowEPPolarizedComptonModel data file <" << ost.str().c_str()
225 << "> is not opened!" << G4endl; 229 << "> is not opened!" << G4endl;
226 G4Exception("G4LowEPPolarizedComptonModel: 230 G4Exception("G4LowEPPolarizedComptonModel::ReadData()",
227 "em0003",FatalException, 231 "em0003",FatalException,
228 ed,"G4LEDATA version should be 232 ed,"G4LEDATA version should be G4EMLOW6.34 or later");
229 return; 233 return;
230 } else { 234 } else {
231 if(verboseLevel > 3) { 235 if(verboseLevel > 3) {
232 G4cout << "File " << ost.str() 236 G4cout << "File " << ost.str()
233 << " is opened by G4LowEPPolari 237 << " is opened by G4LowEPPolarizedComptonModel" << G4endl;
234 } 238 }
235 data[Z]->Retrieve(fin, true); 239 data[Z]->Retrieve(fin, true);
236 data[Z]->ScaleVector(MeV, MeV*barn); 240 data[Z]->ScaleVector(MeV, MeV*barn);
237 } 241 }
238 fin.close(); 242 fin.close();
239 } 243 }
240 244
241 //******************************************** 245 //****************************************************************************
242 246
>> 247
243 G4double 248 G4double
244 G4LowEPPolarizedComptonModel::ComputeCrossSect 249 G4LowEPPolarizedComptonModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
245 250 G4double GammaEnergy,
246 251 G4double Z, G4double,
247 252 G4double, G4double)
248 { 253 {
249 if (verboseLevel > 3) { 254 if (verboseLevel > 3) {
250 G4cout << "G4LowEPPolarizedComptonModel::C 255 G4cout << "G4LowEPPolarizedComptonModel::ComputeCrossSectionPerAtom()"
251 << G4endl; 256 << G4endl;
252 } 257 }
253 G4double cs = 0.0; 258 G4double cs = 0.0;
254 259
255 if (GammaEnergy < LowEnergyLimit()) { return 260 if (GammaEnergy < LowEnergyLimit()) { return 0.0; }
256 261
257 G4int intZ = G4lrint(Z); 262 G4int intZ = G4lrint(Z);
258 if(intZ < 1 || intZ > maxZ) { return cs; } 263 if(intZ < 1 || intZ > maxZ) { return cs; }
259 264
260 G4PhysicsFreeVector* pv = data[intZ]; << 265 G4LPhysicsFreeVector* pv = data[intZ];
261 266
262 // if element was not initialised 267 // if element was not initialised
263 // do initialisation safely for MT mode 268 // do initialisation safely for MT mode
264 if(!pv) 269 if(!pv)
265 { 270 {
266 InitialiseForElement(0, intZ); 271 InitialiseForElement(0, intZ);
267 pv = data[intZ]; 272 pv = data[intZ];
268 if(!pv) { return cs; } 273 if(!pv) { return cs; }
269 } 274 }
270 275
271 G4int n = G4int(pv->GetVectorLength() - 1); << 276 G4int n = pv->GetVectorLength() - 1;
272 G4double e1 = pv->Energy(0); 277 G4double e1 = pv->Energy(0);
273 G4double e2 = pv->Energy(n); 278 G4double e2 = pv->Energy(n);
274 279
275 if(GammaEnergy <= e1) { cs = GammaEnerg 280 if(GammaEnergy <= e1) { cs = GammaEnergy/(e1*e1)*pv->Value(e1); }
276 else if(GammaEnergy <= e2) { cs = pv->Value( 281 else if(GammaEnergy <= e2) { cs = pv->Value(GammaEnergy)/GammaEnergy; }
277 else if(GammaEnergy > e2) { cs = pv->Value( 282 else if(GammaEnergy > e2) { cs = pv->Value(e2)/GammaEnergy; }
278 283
279 return cs; 284 return cs;
280 } 285 }
281 286
282 //******************************************** 287 //****************************************************************************
283 288
284 void G4LowEPPolarizedComptonModel::SampleSecon 289 void G4LowEPPolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
285 const G4MaterialCutsCouple* c << 290 const G4MaterialCutsCouple* couple,
286 const G4DynamicParticle* aDyn << 291 const G4DynamicParticle* aDynamicGamma,
287 G4double, G4double) << 292 G4double, G4double)
288 { 293 {
289 294
290 //Determine number of digits (in decimal bas << 295 //Determine number of digits (in decimal base) that G4double can accurately represent
291 G4double g4d_order = G4double(numeric_limits << 296 G4double g4d_order = G4double(numeric_limits<G4double>::digits10);
292 G4double g4d_limit = std::pow(10.,-g4d_order << 297 G4double g4d_limit = std::pow(10.,-g4d_order);
293 298
294 // The scattered gamma energy is sampled acc 299 // The scattered gamma energy is sampled according to Klein - Nishina formula.
295 // then accepted or rejected depending on th 300 // then accepted or rejected depending on the Scattering Function multiplied
296 // by factor from Klein - Nishina formula. 301 // by factor from Klein - Nishina formula.
297 // Expression of the angular distribution as 302 // Expression of the angular distribution as Klein Nishina
298 // angular and energy distribution and Scatt 303 // angular and energy distribution and Scattering fuctions is taken from
299 // D. E. Cullen "A simple model of photon tr 304 // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
300 // Phys. Res. B 101 (1995). Method of sampli 305 // Phys. Res. B 101 (1995). Method of sampling with form factors is different
301 // data are interpolated while in the articl 306 // data are interpolated while in the article they are fitted.
302 // Reference to the article is from J. Stepa 307 // Reference to the article is from J. Stepanek New Photon, Positron
303 // and Electron Interaction Data for GEANT i 308 // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
304 // TeV (draft). 309 // TeV (draft).
305 // The random number techniques of Butcher & 310 // The random number techniques of Butcher & Messel are used
306 // (Nucl Phys 20(1960),15). 311 // (Nucl Phys 20(1960),15).
307 312
>> 313
308 G4double photonEnergy0 = aDynamicGamma->GetK 314 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy()/MeV;
309 315
310 if (verboseLevel > 3) { 316 if (verboseLevel > 3) {
311 G4cout << "G4LowEPPolarizedComptonModel::S 317 G4cout << "G4LowEPPolarizedComptonModel::SampleSecondaries() E(MeV)= "
312 << photonEnergy0/MeV << " in " << c 318 << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
313 << G4endl; 319 << G4endl;
314 } 320 }
315 // do nothing below the threshold 321 // do nothing below the threshold
316 // should never get here because the XS is z 322 // should never get here because the XS is zero below the limit
317 if (photonEnergy0 < LowEnergyLimit()) 323 if (photonEnergy0 < LowEnergyLimit())
318 return ; 324 return ;
319 325
320 G4double e0m = photonEnergy0 / electron_mass 326 G4double e0m = photonEnergy0 / electron_mass_c2 ;
321 G4ParticleMomentum photonDirection0 = aDynam 327 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
>> 328
322 329
323 // Polarisation: check orientation of photon 330 // Polarisation: check orientation of photon propagation direction and polarisation
324 // Fix if needed 331 // Fix if needed
325 332
326 G4ThreeVector photonPolarization0 = aDynamic 333 G4ThreeVector photonPolarization0 = aDynamicGamma->GetPolarization();
>> 334
327 // Check if polarisation vector is perpendic 335 // Check if polarisation vector is perpendicular and fix if not
328 336
329 if (!(photonPolarization0.isOrthogonal(photo 337 if (!(photonPolarization0.isOrthogonal(photonDirection0, 1e-6))||(photonPolarization0.mag()==0))
330 { 338 {
331 photonPolarization0 = GetRandomPolarizat 339 photonPolarization0 = GetRandomPolarization(photonDirection0);
332 } 340 }
>> 341
333 else 342 else
334 { 343 {
335 if ((photonPolarization0.howOrthogonal(p 344 if ((photonPolarization0.howOrthogonal(photonDirection0) !=0) && (photonPolarization0.howOrthogonal(photonDirection0) > g4d_limit))
336 { 345 {
337 photonPolarization0 = GetPerpendicul 346 photonPolarization0 = GetPerpendicularPolarization(photonDirection0,photonPolarization0);
338 } 347 }
339 } << 348 }
340 349
341 // Select randomly one element in the curren 350 // Select randomly one element in the current material
>> 351
342 const G4ParticleDefinition* particle = aDyn 352 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
343 const G4Element* elm = SelectRandomAtom(coup 353 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
344 G4int Z = (G4int)elm->GetZ(); 354 G4int Z = (G4int)elm->GetZ();
345 355
346 G4double LowEPPCepsilon0 = 1. / (1. + 2. * e 356 G4double LowEPPCepsilon0 = 1. / (1. + 2. * e0m);
347 G4double LowEPPCepsilon0Sq = LowEPPCepsilon0 357 G4double LowEPPCepsilon0Sq = LowEPPCepsilon0 * LowEPPCepsilon0;
348 G4double alpha1 = -std::log(LowEPPCepsilon0) 358 G4double alpha1 = -std::log(LowEPPCepsilon0);
349 G4double alpha2 = 0.5 * (1. - LowEPPCepsilon 359 G4double alpha2 = 0.5 * (1. - LowEPPCepsilon0Sq);
350 360
351 G4double wlPhoton = h_Planck*c_light/photonE 361 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
352 362
353 // Sample the energy of the scattered photon 363 // Sample the energy of the scattered photon
354 G4double LowEPPCepsilon; 364 G4double LowEPPCepsilon;
355 G4double LowEPPCepsilonSq; 365 G4double LowEPPCepsilonSq;
356 G4double oneCosT; 366 G4double oneCosT;
357 G4double sinT2; 367 G4double sinT2;
358 G4double gReject; 368 G4double gReject;
359 369
360 if (verboseLevel > 3) { 370 if (verboseLevel > 3) {
361 G4cout << "Started loop to sample gamma en 371 G4cout << "Started loop to sample gamma energy" << G4endl;
362 } 372 }
363 373
364 do 374 do
365 { 375 {
366 if ( alpha1/(alpha1+alpha2) > G4UniformR 376 if ( alpha1/(alpha1+alpha2) > G4UniformRand())
367 { 377 {
368 LowEPPCepsilon = G4Exp(-alpha1 * G4U 378 LowEPPCepsilon = G4Exp(-alpha1 * G4UniformRand());
369 LowEPPCepsilonSq = LowEPPCepsilon * 379 LowEPPCepsilonSq = LowEPPCepsilon * LowEPPCepsilon;
370 } 380 }
371 else 381 else
372 { 382 {
373 LowEPPCepsilonSq = LowEPPCepsilon0Sq 383 LowEPPCepsilonSq = LowEPPCepsilon0Sq + (1. - LowEPPCepsilon0Sq) * G4UniformRand();
374 LowEPPCepsilon = std::sqrt(LowEPPCep 384 LowEPPCepsilon = std::sqrt(LowEPPCepsilonSq);
375 } 385 }
376 386
377 oneCosT = (1. - LowEPPCepsilon) / ( LowE 387 oneCosT = (1. - LowEPPCepsilon) / ( LowEPPCepsilon * e0m);
378 sinT2 = oneCosT * (2. - oneCosT); 388 sinT2 = oneCosT * (2. - oneCosT);
379 G4double x = std::sqrt(oneCosT/2.) / (wl 389 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
380 G4double scatteringFunction = ComputeSca 390 G4double scatteringFunction = ComputeScatteringFunction(x, Z);
381 gReject = (1. - LowEPPCepsilon * sinT2 / 391 gReject = (1. - LowEPPCepsilon * sinT2 / (1. + LowEPPCepsilonSq)) * scatteringFunction;
382 392
383 } while(gReject < G4UniformRand()*Z); 393 } while(gReject < G4UniformRand()*Z);
384 394
385 G4double cosTheta = 1. - oneCosT; 395 G4double cosTheta = 1. - oneCosT;
386 G4double sinTheta = std::sqrt(sinT2); 396 G4double sinTheta = std::sqrt(sinT2);
387 G4double phi = SetPhi(LowEPPCepsilon,sinT2); 397 G4double phi = SetPhi(LowEPPCepsilon,sinT2);
388 G4double dirx = sinTheta * std::cos(phi); 398 G4double dirx = sinTheta * std::cos(phi);
389 G4double diry = sinTheta * std::sin(phi); 399 G4double diry = sinTheta * std::sin(phi);
390 G4double dirz = cosTheta ; 400 G4double dirz = cosTheta ;
391 401
392 // Set outgoing photon polarization 402 // Set outgoing photon polarization
393 403
394 G4ThreeVector photonPolarization1 = SetNewPo 404 G4ThreeVector photonPolarization1 = SetNewPolarization(LowEPPCepsilon,
395 sinT2, << 405 sinT2,
396 phi, << 406 phi,
397 cosTheta); << 407 cosTheta);
398 408
399 // Scatter photon energy and Compton electro 409 // Scatter photon energy and Compton electron direction - Method based on:
400 // J. M. C. Brown, M. R. Dimmock, J. E. Gill 410 // J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin'
401 // "A low energy bound atomic electron Compt 411 // "A low energy bound atomic electron Compton scattering model for Geant4"
402 // NIMB, Vol. 338, 77-88, 2014. 412 // NIMB, Vol. 338, 77-88, 2014.
403 413
404 // Set constants and initialize scattering p 414 // Set constants and initialize scattering parameters
>> 415
405 const G4double vel_c = c_light / (m/s); 416 const G4double vel_c = c_light / (m/s);
406 const G4double momentum_au_to_nat = halfpi* 417 const G4double momentum_au_to_nat = halfpi* hbar_Planck / Bohr_radius / (kg*m/s);
407 const G4double e_mass_kg = electron_mass_c2 418 const G4double e_mass_kg = electron_mass_c2 / c_squared / kg ;
408 419
409 const G4int maxDopplerIterations = 1000; 420 const G4int maxDopplerIterations = 1000;
410 G4double bindingE = 0.; 421 G4double bindingE = 0.;
411 G4double pEIncident = photonEnergy0 ; 422 G4double pEIncident = photonEnergy0 ;
412 G4double pERecoil = -1.; 423 G4double pERecoil = -1.;
413 G4double eERecoil = -1.; 424 G4double eERecoil = -1.;
414 G4double e_alpha =0.; 425 G4double e_alpha =0.;
415 G4double e_beta = 0.; 426 G4double e_beta = 0.;
416 427
417 G4double CE_emission_flag = 0.; 428 G4double CE_emission_flag = 0.;
418 G4double ePAU = -1; 429 G4double ePAU = -1;
419 G4int shellIdx = 0; 430 G4int shellIdx = 0;
420 G4double u_temp = 0; 431 G4double u_temp = 0;
421 G4double cosPhiE =0; 432 G4double cosPhiE =0;
422 G4double sinThetaE =0; 433 G4double sinThetaE =0;
423 G4double cosThetaE =0; 434 G4double cosThetaE =0;
424 G4int iteration = 0; 435 G4int iteration = 0;
425 436
426 if (verboseLevel > 3) { 437 if (verboseLevel > 3) {
427 G4cout << "Started loop to sample photon e 438 G4cout << "Started loop to sample photon energy and electron direction" << G4endl;
428 } 439 }
429 440
430 do{ 441 do{
431 // *************************************** <<
432 // | Determine scatter photon energy <<
433 // *************************************** <<
434 do <<
435 { <<
436 iteration++; <<
437 442
438 // ***************************************** <<
439 // | Sample bound electron information <<
440 // ***************************************** <<
441 <<
442 // Select shell based on shell occupancy <<
443 shellIdx = shellData->SelectRandomShell(Z); <<
444 bindingE = shellData->BindingEnergy(Z,shellI <<
445 <<
446 // Randomly sample bound electron momentum ( <<
447 ePAU = profileData->RandomSelectMomentum(Z,s <<
448 <<
449 // Convert to SI units <<
450 G4double ePSI = ePAU * momentum_au_to_nat; <<
451 <<
452 //Calculate bound electron velocity and norm <<
453 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / <<
454 <<
455 // Sample incident electron direction, amorp <<
456 e_alpha = pi*G4UniformRand(); <<
457 e_beta = twopi*G4UniformRand(); <<
458 <<
459 // Total energy of system <<
460 G4double eEIncident = electron_mass_c2 / sqr <<
461 G4double systemE = eEIncident + pEIncident; <<
462 <<
463 G4double gamma_temp = 1.0 / sqrt( 1 - (u_tem <<
464 G4double numerator = gamma_temp*electron_mas <<
465 G4double subdenom1 = u_temp*cosTheta*std::c <<
466 G4double subdenom2 = u_temp*sinTheta*std::si <<
467 G4double denominator = (1.0 - cosTheta) + ( <<
468 pERecoil = (numerator/denominator); <<
469 eERecoil = systemE - pERecoil; <<
470 CE_emission_flag = pEIncident - pERecoil; <<
471 } while ( (iteration <= maxDopplerIterat <<
472 <<
473 // End of recalculation of photon energy w <<
474 // *************************************** <<
475 // | Determine ejected Compton electro <<
476 // *************************************** <<
477 <<
478 // Calculate velocity of ejected Compton e <<
479 <<
480 G4double a_temp = eERecoil / electron_mass <<
481 G4double u_p_temp = sqrt(1 - (1 / (a_temp* <<
482 <<
483 // Coefficients and terms from simulatenou <<
484 G4double sinAlpha = std::sin(e_alpha); <<
485 G4double cosAlpha = std::cos(e_alpha); <<
486 G4double sinBeta = std::sin(e_beta); <<
487 G4double cosBeta = std::cos(e_beta); <<
488 <<
489 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_ <<
490 G4double gamma_p = 1.0 / sqrt(1 - (u_p_tem <<
491 <<
492 G4double var_A = pERecoil*u_p_temp*sinThet <<
493 G4double var_B = u_p_temp* (pERecoil*cosTh <<
494 G4double var_C = (pERecoil-pEIncident) - ( <<
495 <<
496 G4double var_D1 = gamma*electron_mass_c2*p <<
497 G4double var_D2 = (1 - (u_temp*cosTheta*co <<
498 G4double var_D3 = ((electron_mass_c2*elect <<
499 G4double var_D = var_D1*var_D2 + var_D3; <<
500 <<
501 G4double var_E1 = ((gamma*gamma_p)*(electr <<
502 G4double var_E2 = gamma_p*electron_mass_c2 <<
503 G4double var_E = var_E1 - var_E2; <<
504 <<
505 G4double var_F1 = ((gamma*gamma_p)*(electr <<
506 G4double var_F2 = (gamma_p*electron_mass_c <<
507 G4double var_F = var_F1 - var_F2; <<
508 <<
509 G4double var_G = (gamma*gamma_p)*(electron <<
510 <<
511 // Two equations form a quadratic form of <<
512 // Coefficents and solution to quadratic <<
513 G4double var_W1 = (var_F*var_B - var_E*var <<
514 G4double var_W2 = (var_G*var_G)*(var_A*var <<
515 G4double var_W = var_W1 + var_W2; <<
516 <<
517 G4double var_Y = 2.0*(((var_A*var_D-var_F* <<
518 <<
519 G4double var_Z1 = (var_A*var_D - var_F*var <<
520 G4double var_Z2 = (var_G*var_G)*(var_C*var <<
521 G4double var_Z = var_Z1 + var_Z2; <<
522 G4double diff1 = var_Y*var_Y; <<
523 G4double diff2 = 4*var_W*var_Z; <<
524 G4double diff = diff1 - diff2; <<
525 <<
526 // Check if diff is less than zero, if so <<
527 //Confirm that diff less than zero is due <<
528 //than 10^(-g4d_order), then set diff to z <<
529 443
530 if ((diff < 0.0) && (abs(diff / diff1) < g << 444 // ******************************************
531 { << 445 // | Determine scatter photon energy |
532 diff = 0.0; << 446 // ******************************************
533 } << 447
>> 448 do
>> 449 {
>> 450 iteration++;
>> 451
>> 452
>> 453 // ********************************************
>> 454 // | Sample bound electron information |
>> 455 // ********************************************
>> 456
>> 457 // Select shell based on shell occupancy
>> 458
>> 459 shellIdx = shellData->SelectRandomShell(Z);
>> 460 bindingE = shellData->BindingEnergy(Z,shellIdx)/MeV;
>> 461
>> 462
>> 463 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
>> 464 ePAU = profileData->RandomSelectMomentum(Z,shellIdx);
>> 465
>> 466 // Convert to SI units
>> 467 G4double ePSI = ePAU * momentum_au_to_nat;
>> 468
>> 469 //Calculate bound electron velocity and normalise to natural units
>> 470 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)) )/vel_c;
>> 471
>> 472 // Sample incident electron direction, amorphous material, to scattering photon scattering plane
>> 473
>> 474 e_alpha = pi*G4UniformRand();
>> 475 e_beta = twopi*G4UniformRand();
>> 476
>> 477 // Total energy of system
>> 478
>> 479 G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
>> 480 G4double systemE = eEIncident + pEIncident;
>> 481
>> 482
>> 483 G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp));
>> 484 G4double numerator = gamma_temp*electron_mass_c2*(1 - u_temp * std::cos(e_alpha));
>> 485 G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha);
>> 486 G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta);
>> 487 G4double denominator = (1.0 - cosTheta) + (gamma_temp*electron_mass_c2*(1 - subdenom1 - subdenom2) / pEIncident);
>> 488 pERecoil = (numerator/denominator);
>> 489 eERecoil = systemE - pERecoil;
>> 490 CE_emission_flag = pEIncident - pERecoil;
>> 491 } while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE));
>> 492
>> 493 // End of recalculation of photon energy with Doppler broadening
>> 494
>> 495
>> 496
>> 497 // *******************************************************
>> 498 // | Determine ejected Compton electron direction |
>> 499 // *******************************************************
>> 500
>> 501 // Calculate velocity of ejected Compton electron
>> 502
>> 503 G4double a_temp = eERecoil / electron_mass_c2;
>> 504 G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp)));
>> 505
>> 506 // Coefficients and terms from simulatenous equations
>> 507
>> 508 G4double sinAlpha = std::sin(e_alpha);
>> 509 G4double cosAlpha = std::cos(e_alpha);
>> 510 G4double sinBeta = std::sin(e_beta);
>> 511 G4double cosBeta = std::cos(e_beta);
>> 512
>> 513 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp));
>> 514 G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp));
534 515
535 // Plus and minus of quadratic << 516 G4double var_A = pERecoil*u_p_temp*sinTheta;
536 G4double X_p = (-var_Y + sqrt (diff))/(2*v << 517 G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident);
537 G4double X_m = (-var_Y - sqrt (diff))/(2*v << 518 G4double var_C = (pERecoil-pEIncident) - ( (pERecoil*pEIncident) / (gamma_p*electron_mass_c2))*(1 - cosTheta);
538 <<
539 // Floating point precision protection <<
540 // Check if X_p and X_m are greater than o <<
541 // Issue due to propagation of FPE and onl <<
542 if(X_p >1){X_p=1;} if(X_p<-1){X_p=-1;} <<
543 if(X_m >1){X_m=1;} if(X_m<-1){X_m=-1;} <<
544 // End of FP protection <<
545 519
546 G4double ThetaE = 0.; << 520 G4double var_D1 = gamma*electron_mass_c2*pERecoil;
>> 521 G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha));
>> 522 G4double var_D3 = ((electron_mass_c2*electron_mass_c2)*(gamma*gamma_p - 1)) - (gamma_p*electron_mass_c2*pERecoil);
>> 523 G4double var_D = var_D1*var_D2 + var_D3;
547 524
548 // Randomly sample one of the two possible << 525 G4double var_E1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosAlpha);
549 G4double sol_select = G4UniformRand(); << 526 G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta;
>> 527 G4double var_E = var_E1 - var_E2;
550 528
551 if (sol_select < 0.5) << 529 G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha);
>> 530 G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta);
>> 531 G4double var_F = var_F1 - var_F2;
>> 532
>> 533 G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha;
>> 534
>> 535 // Two equations form a quadratic form of Wx^2 + Yx + Z = 0
>> 536 // Coefficents and solution to quadratic
>> 537
>> 538 G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A);
>> 539 G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B);
>> 540 G4double var_W = var_W1 + var_W2;
>> 541
>> 542 G4double var_Y = 2.0*(((var_A*var_D-var_F*var_C)*(var_F*var_B-var_E*var_A)) - ((var_G*var_G)*var_B*var_C));
>> 543
>> 544 G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C);
>> 545 G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A);
>> 546 G4double var_Z = var_Z1 + var_Z2;
>> 547 G4double diff1 = var_Y*var_Y;
>> 548 G4double diff2 = 4*var_W*var_Z;
>> 549 G4double diff = diff1 - diff2;
>> 550
>> 551
>> 552 // Check if diff is less than zero, if so ensure it is due to FPE
>> 553
>> 554 //Confirm that diff less than zero is due FPE, i.e if abs of diff / diff1 and diff/ diff2 is less
>> 555 //than 10^(-g4d_order), then set diff to zero
>> 556
>> 557 if ((diff < 0.0) && (abs(diff / diff1) < g4d_limit) && (abs(diff / diff2) < g4d_limit) )
>> 558 {
>> 559 diff = 0.0;
>> 560 }
>> 561
>> 562 // Plus and minus of quadratic
>> 563 G4double X_p = (-var_Y + sqrt (diff))/(2*var_W);
>> 564 G4double X_m = (-var_Y - sqrt (diff))/(2*var_W);
>> 565
>> 566
>> 567 // Floating point precision protection
>> 568 // Check if X_p and X_m are greater than or less than 1 or -1, if so clean up FPE
>> 569 // Issue due to propagation of FPE and only impacts 8th sig fig onwards
>> 570
>> 571 if(X_p >1){X_p=1;} if(X_p<-1){X_p=-1;}
>> 572 if(X_m >1){X_m=1;} if(X_m<-1){X_m=-1;}
>> 573
>> 574 // End of FP protection
>> 575
>> 576 G4double ThetaE = 0.;
>> 577
>> 578
>> 579 // Randomly sample one of the two possible solutions and determin theta angle of ejected Compton electron
>> 580 G4double sol_select = G4UniformRand();
>> 581
>> 582 if (sol_select < 0.5)
552 { 583 {
553 ThetaE = std::acos(X_p); << 584 ThetaE = std::acos(X_p);
554 } 585 }
555 if (sol_select > 0.5) << 586 if (sol_select > 0.5)
556 { 587 {
557 ThetaE = std::acos(X_m); << 588 ThetaE = std::acos(X_m);
558 } 589 }
559 590
560 cosThetaE = std::cos(ThetaE); << 591 cosThetaE = std::cos(ThetaE);
561 sinThetaE = std::sin(ThetaE); << 592 sinThetaE = std::sin(ThetaE);
562 G4double Theta = std::acos(cosTheta); << 593 G4double Theta = std::acos(cosTheta);
563 << 594
564 //Calculate electron Phi << 595 //Calculate electron Phi
565 G4double iSinThetaE = std::sqrt(1+std::tan << 596 G4double iSinThetaE = std::sqrt(1+std::tan((pi/2.0)-ThetaE)*std::tan((pi/2.0)-ThetaE));
566 G4double iSinTheta = std::sqrt(1+std::tan( << 597 G4double iSinTheta = std::sqrt(1+std::tan((pi/2.0)-Theta)*std::tan((pi/2.0)-Theta));
567 G4double ivar_A = iSinTheta/ (pERecoil*u_p << 598 G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp);
568 // Trigs << 599 // Trigs
569 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_ << 600 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE);
570 // End of calculation of ejection Compton << 601
571 //Fix for floating point errors << 602 // End of calculation of ejection Compton electron direction
572 } while ( (iteration <= maxDopplerIterations << 603
>> 604 //Fix for floating point errors
>> 605
>> 606 } while ( (iteration <= maxDopplerIterations) && (abs(cosPhiE) > 1));
573 607
574 // Revert to original if maximum number of i << 608 // Revert to original if maximum number of iterations threshold has been reached
575 if (iteration >= maxDopplerIterations) 609 if (iteration >= maxDopplerIterations)
576 { 610 {
577 pERecoil = photonEnergy0 ; 611 pERecoil = photonEnergy0 ;
578 bindingE = 0.; 612 bindingE = 0.;
579 dirx=0.0; 613 dirx=0.0;
580 diry=0.0; 614 diry=0.0;
581 dirz=1.0; 615 dirz=1.0;
582 } 616 }
583 617
584 // Set "scattered" photon direction and ener 618 // Set "scattered" photon direction and energy
>> 619
585 G4ThreeVector photonDirection1(dirx,diry,dir 620 G4ThreeVector photonDirection1(dirx,diry,dirz);
586 SystemOfRefChange(photonDirection0,photonDir 621 SystemOfRefChange(photonDirection0,photonDirection1,
587 photonPolarization0,photonPolarization 622 photonPolarization0,photonPolarization1);
588 623
>> 624
589 if (pERecoil > 0.) 625 if (pERecoil > 0.)
590 { 626 {
591 fParticleChange->SetProposedKineticEnerg << 627 fParticleChange->SetProposedKineticEnergy(pERecoil) ;
592 fParticleChange->ProposeMomentumDirectio << 628 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
593 fParticleChange->ProposePolarization(pho << 629 fParticleChange->ProposePolarization(photonPolarization1);
594 << 630
595 // Set ejected Compton electron directio << 631 // Set ejected Compton electron direction and energy
596 G4double PhiE = std::acos(cosPhiE); << 632 G4double PhiE = std::acos(cosPhiE);
597 G4double eDirX = sinThetaE * std::cos(ph << 633 G4double eDirX = sinThetaE * std::cos(phi+PhiE);
598 G4double eDirY = sinThetaE * std::sin(ph << 634 G4double eDirY = sinThetaE * std::sin(phi+PhiE);
599 G4double eDirZ = cosThetaE; << 635 G4double eDirZ = cosThetaE;
600 << 636
601 G4double eKineticEnergy = pEIncident - p << 637 G4double eKineticEnergy = pEIncident - pERecoil - bindingE;
602 << 638
603 G4ThreeVector eDirection(eDirX,eDirY,eDi << 639 G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
604 SystemOfRefChangeElect(photonDirection0, << 640 SystemOfRefChangeElect(photonDirection0,eDirection,
605 photonPolarization0); << 641 photonPolarization0);
606 << 642
607 G4DynamicParticle* dp = new G4DynamicPar << 643 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
608 eDirection,eKineticEnergy) ; << 644 eDirection,eKineticEnergy) ;
609 fvect->push_back(dp); << 645 fvect->push_back(dp);
>> 646
610 } 647 }
611 else 648 else
612 { 649 {
613 fParticleChange->SetProposedKineticEnerg 650 fParticleChange->SetProposedKineticEnergy(0.);
614 fParticleChange->ProposeTrackStatus(fSto 651 fParticleChange->ProposeTrackStatus(fStopAndKill);
615 } 652 }
616 653
617 // sample deexcitation 654 // sample deexcitation
618 // 655 //
>> 656
619 if (verboseLevel > 3) { 657 if (verboseLevel > 3) {
620 G4cout << "Started atomic de-excitation " 658 G4cout << "Started atomic de-excitation " << fAtomDeexcitation << G4endl;
621 } 659 }
622 660
623 if(fAtomDeexcitation && iteration < maxDoppl 661 if(fAtomDeexcitation && iteration < maxDopplerIterations) {
624 G4int index = couple->GetIndex(); 662 G4int index = couple->GetIndex();
625 if(fAtomDeexcitation->CheckDeexcitationAct 663 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
626 std::size_t nbefore = fvect->size(); << 664 size_t nbefore = fvect->size();
627 G4AtomicShellEnumerator as = G4AtomicShe 665 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
628 const G4AtomicShell* shell = fAtomDeexci 666 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
629 fAtomDeexcitation->GenerateParticles(fve 667 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
630 std::size_t nafter = fvect->size(); << 668 size_t nafter = fvect->size();
631 if(nafter > nbefore) { 669 if(nafter > nbefore) {
632 for (std::size_t i=nbefore; i<nafter; << 670 for (size_t i=nbefore; i<nafter; ++i) {
633 //Check if there is enough residual 671 //Check if there is enough residual energy
634 if (bindingE >= ((*fvect)[i])->GetKi 672 if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
635 { << 673 {
636 //Ok, this is a valid secondary: keep << 674 //Ok, this is a valid secondary: keep it
637 bindingE -= ((*fvect)[i])->GetKineticE << 675 bindingE -= ((*fvect)[i])->GetKineticEnergy();
638 } << 676 }
639 else 677 else
640 { << 678 {
641 //Invalid secondary: not enough energy << 679 //Invalid secondary: not enough energy to create it!
642 //Keep its energy in the local deposit << 680 //Keep its energy in the local deposit
643 delete (*fvect)[i]; << 681 delete (*fvect)[i];
644 (*fvect)[i]=nullptr; << 682 (*fvect)[i]=0;
645 } << 683 }
646 } 684 }
647 } 685 }
648 } 686 }
649 } 687 }
650 688
651 //This should never happen 689 //This should never happen
652 if(bindingE < 0.0) 690 if(bindingE < 0.0)
653 G4Exception("G4LowEPPolarizedComptonModel: << 691 G4Exception("G4LowEPPolarizedComptonModel::SampleSecondaries()",
654 "em2051",FatalException,"Negative local en << 692 "em2051",FatalException,"Negative local energy deposit");
655 693
656 fParticleChange->ProposeLocalEnergyDeposit(b 694 fParticleChange->ProposeLocalEnergyDeposit(bindingE);
>> 695
657 } 696 }
658 697
659 //******************************************** 698 //****************************************************************************
660 699
661 G4double 700 G4double
662 G4LowEPPolarizedComptonModel::ComputeScatterin 701 G4LowEPPolarizedComptonModel::ComputeScatteringFunction(G4double x, G4int Z)
663 { 702 {
664 G4double value = Z; 703 G4double value = Z;
665 if (x <= ScatFuncFitParam[Z][2]) { 704 if (x <= ScatFuncFitParam[Z][2]) {
666 705
667 G4double lgq = G4Log(x)/ln10; 706 G4double lgq = G4Log(x)/ln10;
668 707
669 if (lgq < ScatFuncFitParam[Z][1]) { 708 if (lgq < ScatFuncFitParam[Z][1]) {
670 value = ScatFuncFitParam[Z][3] + lgq*Sca 709 value = ScatFuncFitParam[Z][3] + lgq*ScatFuncFitParam[Z][4];
671 } else { 710 } else {
672 value = ScatFuncFitParam[Z][5] + lgq*Sca 711 value = ScatFuncFitParam[Z][5] + lgq*ScatFuncFitParam[Z][6] +
673 lgq*lgq*ScatFuncFitParam[Z][7] + lgq*l 712 lgq*lgq*ScatFuncFitParam[Z][7] + lgq*lgq*lgq*ScatFuncFitParam[Z][8];
674 } 713 }
675 value = G4Exp(value*ln10); 714 value = G4Exp(value*ln10);
676 } 715 }
677 return value; 716 return value;
678 } 717 }
679 718
680 719
681 //******************************************** 720 //****************************************************************************
682 721
>> 722 #include "G4AutoLock.hh"
>> 723 namespace { G4Mutex LowEPPolarizedComptonModelMutex = G4MUTEX_INITIALIZER; }
>> 724
683 void 725 void
684 G4LowEPPolarizedComptonModel::InitialiseForEle 726 G4LowEPPolarizedComptonModel::InitialiseForElement(const G4ParticleDefinition*,
685 G4int Z) << 727 G4int Z)
686 { 728 {
687 G4AutoLock l(&LowEPPolarizedComptonModelMute 729 G4AutoLock l(&LowEPPolarizedComptonModelMutex);
688 if(!data[Z]) { ReadData(Z); } 730 if(!data[Z]) { ReadData(Z); }
689 l.unlock(); 731 l.unlock();
690 } 732 }
691 733
692 //******************************************** 734 //****************************************************************************
693 735
694 //Fitting data to compute scattering function 736 //Fitting data to compute scattering function
>> 737
695 const G4double G4LowEPPolarizedComptonModel::S 738 const G4double G4LowEPPolarizedComptonModel::ScatFuncFitParam[101][9] = {
696 { 0, 0., 0., 0., 0., << 739 { 0, 0., 0., 0., 0., 0., 0., 0., 0.},
697 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -1 << 740 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -143.374, 50.787, -5.951, 0.2304418},
698 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -5 << 741 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -53.649, 13.892, -0.948, 0.006996759},
699 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -6 << 742 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -62.090, 21.462, -2.453, 0.093416},
700 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -1 << 743 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -127.906, 46.491, -5.614, 0.2262103},
701 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -1 << 744 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -131.153, 47.132, -5.619, 0.2233819},
702 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -1 << 745 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -128.143, 45.379, -5.325, 0.2083009},
703 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -1 << 746 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -131.048, 46.314, -5.421, 0.2114925},
704 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -1 << 747 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -128.225, 44.818, -5.183, 0.1997155},
705 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -1 << 748 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -122.112, 42.103, -4.796, 0.1819099},
706 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -1 << 749 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -110.143, 37.225, -4.143, 0.1532094},
707 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -4 << 750 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -41.137, 12.313, -1.152, 0.03384553},
708 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -5 << 751 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -53.549, 17.420, -1.840, 0.06431849},
709 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -6 << 752 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -66.243, 22.297, -2.460, 0.09045854},
710 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -7 << 753 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -78.271, 26.757, -3.008, 0.1128195},
711 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -8 << 754 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -85.069, 29.164, -3.291, 0.1239113},
712 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -9 << 755 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -93.640, 32.274, -3.665, 0.1388633},
713 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -9 << 756 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -98.534, 33.958, -3.857, 0.1461557},
714 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -1 << 757 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -100.077, 34.379, -3.891, 0.1468902},
715 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -5 << 758 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -53.819, 17.528, -1.851, 0.0648722},
716 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -5 << 759 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -52.221, 17.169, -1.832, 0.06502094},
717 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -5 << 760 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -55.365, 18.276, -1.961, 0.07002778},
718 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -5 << 761 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -57.412, 18.957, -2.036, 0.07278856},
719 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -5 << 762 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -58.478, 19.270, -2.065, 0.07362722},
720 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -6 << 763 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -62.192, 20.358, -2.167, 0.07666583},
721 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -5 << 764 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -58.007, 18.924, -2.003, 0.0704305},
722 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -6 << 765 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -61.176, 20.067, -2.141, 0.0760269},
723 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -6 << 766 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -61.909, 20.271, -2.159, 0.07653559},
724 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -6 << 767 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -62.402, 20.391, -2.167, 0.07664243},
725 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -6 << 768 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -67.305, 21.954, -2.331, 0.0823267},
726 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -6 << 769 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -62.064, 20.136, -2.122, 0.07437589},
727 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -6 << 770 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -61.068, 19.808, -2.086, 0.07307488},
728 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -6 << 771 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -63.126, 20.553, -2.175, 0.07660222},
729 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -6 << 772 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -65.674, 21.445, -2.278, 0.08054694},
730 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -6 << 773 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -69.457, 22.811, -2.442, 0.08709536},
731 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -7 << 774 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -72.283, 23.808, -2.558, 0.09156808},
732 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -7 << 775 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -74.696, 24.641, -2.653, 0.09516597},
733 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -4 << 776 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -46.235, 14.519, -1.458, 0.04837659},
734 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -4 << 777 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -41.784, 13.065, -1.300, 0.04267703},
735 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -4 << 778 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -44.609, 14.114, -1.429, 0.0479348},
736 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -4 << 779 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -47.142, 15.019, -1.536, 0.0521347},
737 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -5 << 780 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -53.114, 17.052, -1.766, 0.06079426},
738 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -5 << 781 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -54.590, 17.550, -1.822, 0.06290335},
739 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -5 << 782 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -51.272, 16.423, -1.694, 0.05806108},
740 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -5 << 783 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -58.314, 18.839, -1.969, 0.0684608},
741 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -5 << 784 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -59.674, 19.295, -2.020, 0.07037294},
742 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -7 << 785 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -72.228, 23.693, -2.532, 0.09017969},
743 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -6 << 786 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -60.890, 19.647, -2.053, 0.07138694},
744 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -5 << 787 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -56.152, 18.002, -1.863, 0.06410123},
745 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -5 << 788 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -56.208, 18.052, -1.872, 0.06456884},
746 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -5 << 789 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -58.478, 18.887, -1.973, 0.06860356},
747 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -6 << 790 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -60.708, 19.676, -2.066, 0.07225841},
748 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -6 << 791 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -63.341, 20.632, -2.180, 0.0767412},
749 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -6 << 792 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -66.339, 21.716, -2.310, 0.08191981},
750 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -6 << 793 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -67.649, 22.151, -2.357, 0.08357825},
751 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -4 << 794 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -45.302, 14.219, -1.423, 0.04712317},
752 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -3 << 795 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -39.825, 12.363, -1.214, 0.03931009},
753 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -3 << 796 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -38.952, 11.982, -1.160, 0.03681554},
754 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -4 << 797 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -41.959, 13.118, -1.302, 0.04271291},
755 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -4 << 798 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -40.528, 12.555, -1.230, 0.03971407},
756 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -3 << 799 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -39.986, 12.329, -1.200, 0.03843737},
757 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -4 << 800 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -42.756, 13.362, -1.327, 0.0436124},
758 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -4 << 801 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -42.682, 13.314, -1.319, 0.04322509},
759 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -4 << 802 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -42.399, 13.185, -1.301, 0.04243861},
760 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -4 << 803 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -43.226, 13.475, -1.335, 0.04377341},
761 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -4 << 804 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -43.232, 13.456, -1.330, 0.04347536},
762 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -4 << 805 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -42.047, 12.990, -1.270, 0.04095499},
763 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -4 << 806 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -42.405, 13.106, -1.283, 0.04146024},
764 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -4 << 807 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -41.974, 12.926, -1.259, 0.040435},
765 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -4 << 808 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -42.132, 12.967, -1.262, 0.04048908},
766 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -4 << 809 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -42.582, 13.122, -1.280, 0.04119599},
767 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -4 << 810 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -42.586, 13.115, -1.278, 0.04107587},
768 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -4 << 811 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -43.708, 13.509, -1.324, 0.04286491},
769 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -4 << 812 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -44.838, 13.902, -1.369, 0.04457132},
770 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -4 << 813 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -45.545, 14.137, -1.395, 0.04553459},
771 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -4 << 814 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -46.426, 14.431, -1.428, 0.04678218},
772 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -4 << 815 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -47.513, 14.806, -1.471, 0.04842566},
773 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -4 << 816 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -48.225, 15.042, -1.497, 0.04938364},
774 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -5 << 817 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -53.256, 16.739, -1.687, 0.05645173},
775 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -5 << 818 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -53.900, 16.946, -1.709, 0.05723134},
776 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -4 << 819 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -49.801, 15.536, -1.547, 0.05103522},
777 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -4 << 820 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -49.651, 15.512, -1.548, 0.05123203},
778 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -5 << 821 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -51.021, 16.018, -1.609, 0.05364831},
779 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -5 << 822 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -52.693, 16.612, -1.679, 0.05638698},
780 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -5 << 823 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -54.415, 17.238, -1.754, 0.05935566},
781 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -5 << 824 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -56.083, 17.834, -1.824, 0.06206068},
782 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -5 << 825 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -57.860, 18.463, -1.898, 0.0649633},
783 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -3 << 826 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -39.973, 12.164, -1.162, 0.0364598},
784 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -3 << 827 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -34.591, 10.338, -0.956, 0.0287409},
785 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -3 << 828 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -36.002, 10.867, -1.021, 0.03136835},
786 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -3 << 829 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -37.660, 11.475, -1.095, 0.03435334},
787 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -3 << 830 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -37.268, 11.301, -1.071, 0.0330539},
788 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -3 << 831 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -37.695, 11.438, -1.085, 0.03376669},
789 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -3 << 832 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -39.010, 11.927, -1.146, 0.03630848},
790 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -3 << 833 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -37.192, 11.229, -1.057, 0.0325621},
791 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -3 << 834 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -37.589, 11.363, -1.072, 0.03312393},
792 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -3 << 835 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -39.573, 12.095, -1.162, 0.03680527},
793 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -4 << 836 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -40.007, 12.242, -1.178, 0.03737377},
794 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -3 << 837 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -39.509, 12.041, -1.152, 0.03629023},
795 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -3 << 838 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -39.939, 12.183, -1.168, 0.03690464},
796 {100, 6.412, 1.00000E+10, -12.900, 1.993, -3 << 839 {100, 6.412, 1.00000E+10, -12.900, 1.993, -39.973, 12.180, -1.166, 0.036773}
797 }; << 840 };
798 841
799 //******************************************** 842 //****************************************************************************
800 843
801 //Supporting functions for photon polarisation 844 //Supporting functions for photon polarisation effects
>> 845
802 G4double G4LowEPPolarizedComptonModel::SetPhi( 846 G4double G4LowEPPolarizedComptonModel::SetPhi(G4double energyRate,
803 G4double sinT2) << 847 G4double sinT2)
804 { 848 {
805 G4double rand1; 849 G4double rand1;
806 G4double rand2; 850 G4double rand2;
807 G4double phiProbability; 851 G4double phiProbability;
808 G4double phi; 852 G4double phi;
809 G4double a, b; 853 G4double a, b;
810 854
811 do 855 do
812 { 856 {
813 rand1 = G4UniformRand(); 857 rand1 = G4UniformRand();
814 rand2 = G4UniformRand(); 858 rand2 = G4UniformRand();
815 phiProbability=0.; 859 phiProbability=0.;
816 phi = twopi*rand1; 860 phi = twopi*rand1;
817 861
818 a = 2*sinT2; 862 a = 2*sinT2;
819 b = energyRate + 1/energyRate; 863 b = energyRate + 1/energyRate;
820 864
821 phiProbability = 1 - (a/b)*(std::cos(phi 865 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
>> 866
>> 867
>> 868
822 } 869 }
823 while ( rand2 > phiProbability ); 870 while ( rand2 > phiProbability );
824 return phi; 871 return phi;
825 } 872 }
826 873
827 //******************************************** 874 //****************************************************************************
828 875
829 G4ThreeVector G4LowEPPolarizedComptonModel::Se 876 G4ThreeVector G4LowEPPolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a)
830 { 877 {
831 G4double dx = a.x(); 878 G4double dx = a.x();
832 G4double dy = a.y(); 879 G4double dy = a.y();
833 G4double dz = a.z(); 880 G4double dz = a.z();
834 G4double x = dx < 0.0 ? -dx : dx; 881 G4double x = dx < 0.0 ? -dx : dx;
835 G4double y = dy < 0.0 ? -dy : dy; 882 G4double y = dy < 0.0 ? -dy : dy;
836 G4double z = dz < 0.0 ? -dz : dz; 883 G4double z = dz < 0.0 ? -dz : dz;
837 if (x < y) { 884 if (x < y) {
838 return x < z ? G4ThreeVector(-dy,dx,0) : G 885 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
839 }else{ 886 }else{
840 return y < z ? G4ThreeVector(dz,0,-dx) : G 887 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
841 } 888 }
842 } 889 }
843 890
844 //******************************************** 891 //****************************************************************************
845 892
846 G4ThreeVector G4LowEPPolarizedComptonModel::Ge 893 G4ThreeVector G4LowEPPolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0)
847 { 894 {
848 G4ThreeVector d0 = direction0.unit(); 895 G4ThreeVector d0 = direction0.unit();
849 G4ThreeVector a1 = SetPerpendicularVector(d0 896 G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
850 G4ThreeVector a0 = a1.unit(); // unit vector 897 G4ThreeVector a0 = a1.unit(); // unit vector
851 898
852 G4double rand1 = G4UniformRand(); 899 G4double rand1 = G4UniformRand();
853 900
854 G4double angle = twopi*rand1; // random pola 901 G4double angle = twopi*rand1; // random polar angle
855 G4ThreeVector b0 = d0.cross(a0); // cross pr 902 G4ThreeVector b0 = d0.cross(a0); // cross product
856 903
857 G4ThreeVector c; 904 G4ThreeVector c;
858 905
859 c.setX(std::cos(angle)*(a0.x())+std::sin(ang 906 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
860 c.setY(std::cos(angle)*(a0.y())+std::sin(ang 907 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
861 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang 908 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
862 909
863 G4ThreeVector c0 = c.unit(); 910 G4ThreeVector c0 = c.unit();
864 911
865 return c0; 912 return c0;
>> 913
866 } 914 }
867 915
868 //******************************************** 916 //****************************************************************************
869 917
870 G4ThreeVector G4LowEPPolarizedComptonModel::Ge 918 G4ThreeVector G4LowEPPolarizedComptonModel::GetPerpendicularPolarization
871 (const G4ThreeVector& photonDirection, const G 919 (const G4ThreeVector& photonDirection, const G4ThreeVector& photonPolarization) const
872 { 920 {
>> 921
873 // 922 //
874 // The polarization of a photon is always pe 923 // The polarization of a photon is always perpendicular to its momentum direction.
875 // Therefore this function removes those vec 924 // Therefore this function removes those vector component of photonPolarization, which
876 // points in direction of photonDirection 925 // points in direction of photonDirection
877 // 926 //
878 // Mathematically we search the projection o 927 // Mathematically we search the projection of the vector a on the plane E, where n is the
879 // plains normal vector. 928 // plains normal vector.
880 // The basic equation can be found in each g 929 // The basic equation can be found in each geometry book (e.g. Bronstein):
881 // p = a - (a o n)/(n o n)*n 930 // p = a - (a o n)/(n o n)*n
882 931
883 return photonPolarization - photonPolarizati 932 return photonPolarization - photonPolarization.dot(photonDirection)/photonDirection.dot(photonDirection) * photonDirection;
884 } 933 }
885 934
886 //******************************************** 935 //****************************************************************************
887 936
888 G4ThreeVector G4LowEPPolarizedComptonModel::Se 937 G4ThreeVector G4LowEPPolarizedComptonModel::SetNewPolarization(G4double LowEPPCepsilon,
889 G4double sinT2, << 938 G4double sinT2,
890 G4double phi, << 939 G4double phi,
891 G4double costheta) << 940 G4double costheta)
892 { 941 {
893 G4double rand1; 942 G4double rand1;
894 G4double rand2; 943 G4double rand2;
895 G4double cosPhi = std::cos(phi); 944 G4double cosPhi = std::cos(phi);
896 G4double sinPhi = std::sin(phi); 945 G4double sinPhi = std::sin(phi);
897 G4double sinTheta = std::sqrt(sinT2); 946 G4double sinTheta = std::sqrt(sinT2);
898 G4double cosP2 = cosPhi*cosPhi; 947 G4double cosP2 = cosPhi*cosPhi;
899 G4double normalisation = std::sqrt(1. - cosP 948 G4double normalisation = std::sqrt(1. - cosP2*sinT2);
900 949
>> 950
901 // Method based on: 951 // Method based on:
902 // D. Xu, Z. He and F. Zhang 952 // D. Xu, Z. He and F. Zhang
903 // "Detection of Gamma Ray Polarization Usin 953 // "Detection of Gamma Ray Polarization Using a 3-D Position Sensitive CdZnTe Detector"
904 // IEEE TNS, Vol. 52(4), 1160-1164, 2005. 954 // IEEE TNS, Vol. 52(4), 1160-1164, 2005.
905 955
906 // Determination of Theta 956 // Determination of Theta
907 957
908 G4double theta; 958 G4double theta;
>> 959
909 rand1 = G4UniformRand(); 960 rand1 = G4UniformRand();
910 rand2 = G4UniformRand(); 961 rand2 = G4UniformRand();
911 962
912 if (rand1<(LowEPPCepsilon+1.0/LowEPPCepsilon 963 if (rand1<(LowEPPCepsilon+1.0/LowEPPCepsilon-2)/(2.0*(LowEPPCepsilon+1.0/LowEPPCepsilon)-4.0*sinT2*cosP2))
913 { 964 {
914 if (rand2<0.5) 965 if (rand2<0.5)
915 theta = pi/2.0; 966 theta = pi/2.0;
916 else 967 else
917 theta = 3.0*pi/2.0; 968 theta = 3.0*pi/2.0;
918 } 969 }
919 else 970 else
920 { 971 {
921 if (rand2<0.5) 972 if (rand2<0.5)
922 theta = 0; 973 theta = 0;
923 else 974 else
924 theta = pi; 975 theta = pi;
925 } 976 }
926 G4double cosBeta = std::cos(theta); 977 G4double cosBeta = std::cos(theta);
927 G4double sinBeta = std::sqrt(1-cosBeta*cosBe 978 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
928 979
929 G4ThreeVector photonPolarization1; 980 G4ThreeVector photonPolarization1;
930 981
931 G4double xParallel = normalisation*cosBeta; 982 G4double xParallel = normalisation*cosBeta;
932 G4double yParallel = -(sinT2*cosPhi*sinPhi)* 983 G4double yParallel = -(sinT2*cosPhi*sinPhi)*cosBeta/normalisation;
933 G4double zParallel = -(costheta*sinTheta*cos 984 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
934 G4double xPerpendicular = 0.; 985 G4double xPerpendicular = 0.;
935 G4double yPerpendicular = (costheta)*sinBeta 986 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
936 G4double zPerpendicular = -(sinTheta*sinPhi) 987 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
937 988
938 G4double xTotal = (xParallel + xPerpendicula 989 G4double xTotal = (xParallel + xPerpendicular);
939 G4double yTotal = (yParallel + yPerpendicula 990 G4double yTotal = (yParallel + yPerpendicular);
940 G4double zTotal = (zParallel + zPerpendicula 991 G4double zTotal = (zParallel + zPerpendicular);
941 992
942 photonPolarization1.setX(xTotal); 993 photonPolarization1.setX(xTotal);
943 photonPolarization1.setY(yTotal); 994 photonPolarization1.setY(yTotal);
944 photonPolarization1.setZ(zTotal); 995 photonPolarization1.setZ(zTotal);
945 996
946 return photonPolarization1; 997 return photonPolarization1;
947 998
948 } 999 }
949 <<
950 //******************************************** <<
951 void G4LowEPPolarizedComptonModel::SystemOfRef 1000 void G4LowEPPolarizedComptonModel::SystemOfRefChange(G4ThreeVector& direction0,
952 G4ThreeVector& direction1, << 1001 G4ThreeVector& direction1,
953 G4ThreeVector& polarization0, << 1002 G4ThreeVector& polarization0,
954 G4ThreeVector& polarization1) << 1003 G4ThreeVector& polarization1)
955 { 1004 {
956 // direction0 is the original photon directi 1005 // direction0 is the original photon direction ---> z
957 // polarization0 is the original photon pola 1006 // polarization0 is the original photon polarization ---> x
958 // need to specify y axis in the real refere 1007 // need to specify y axis in the real reference frame ---> y
959 G4ThreeVector Axis_Z0 = direction0.unit(); 1008 G4ThreeVector Axis_Z0 = direction0.unit();
960 G4ThreeVector Axis_X0 = polarization0.unit() 1009 G4ThreeVector Axis_X0 = polarization0.unit();
961 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_ 1010 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
962 1011
963 G4double direction_x = direction1.getX(); 1012 G4double direction_x = direction1.getX();
964 G4double direction_y = direction1.getY(); 1013 G4double direction_y = direction1.getY();
965 G4double direction_z = direction1.getZ(); 1014 G4double direction_z = direction1.getZ();
966 1015
967 direction1 = (direction_x*Axis_X0 + directio 1016 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
968 G4double polarization_x = polarization1.getX 1017 G4double polarization_x = polarization1.getX();
969 G4double polarization_y = polarization1.getY 1018 G4double polarization_y = polarization1.getY();
970 G4double polarization_z = polarization1.getZ 1019 G4double polarization_z = polarization1.getZ();
971 1020
972 polarization1 = (polarization_x*Axis_X0 + po 1021 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
973 1022
974 } 1023 }
975 1024
976 //******************************************** <<
977 void G4LowEPPolarizedComptonModel::SystemOfRef 1025 void G4LowEPPolarizedComptonModel::SystemOfRefChangeElect(G4ThreeVector& pdirection,
978 G4ThreeVector& edirection, << 1026 G4ThreeVector& edirection,
979 G4ThreeVector& ppolarization) << 1027 G4ThreeVector& ppolarization)
980 { 1028 {
981 // direction0 is the original photon directi 1029 // direction0 is the original photon direction ---> z
982 // polarization0 is the original photon pola 1030 // polarization0 is the original photon polarization ---> x
983 // need to specify y axis in the real refere 1031 // need to specify y axis in the real reference frame ---> y
984 G4ThreeVector Axis_Z0 = pdirection.unit(); 1032 G4ThreeVector Axis_Z0 = pdirection.unit();
985 G4ThreeVector Axis_X0 = ppolarization.unit() 1033 G4ThreeVector Axis_X0 = ppolarization.unit();
986 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_ 1034 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
987 1035
988 G4double direction_x = edirection.getX(); 1036 G4double direction_x = edirection.getX();
989 G4double direction_y = edirection.getY(); 1037 G4double direction_y = edirection.getY();
990 G4double direction_z = edirection.getZ(); 1038 G4double direction_z = edirection.getZ();
991 1039
992 edirection = (direction_x*Axis_X0 + directio 1040 edirection = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
>> 1041
993 } 1042 }
994 1043
995 1044
996 1045