| 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 // G4SPSAngDistribution class implementation << 26 ///////////////////////////////////////////////////////////////////////////////
>> 27 //
>> 28 // MODULE: G4SPSAngDistribution.cc
>> 29 //
>> 30 // Version: 1.0
>> 31 // Date: 5/02/04
>> 32 // Author: Fan Lei
>> 33 // Organisation: QinetiQ ltd.
>> 34 // Customer: ESA/ESTEC
>> 35 //
>> 36 ///////////////////////////////////////////////////////////////////////////////
>> 37 //
>> 38 //
>> 39 // CHANGE HISTORY
>> 40 // --------------
>> 41 //
>> 42 //
>> 43 // Version 1.0, 05/02/2004, Fan Lei, Created.
>> 44 // Based on the G4GeneralParticleSource class in Geant4 v6.0
>> 45 //
>> 46 ///////////////////////////////////////////////////////////////////////////////
27 // 47 //
28 // Author: Fan Lei, QinetiQ ltd. - 05/02/2004 <<
29 // Customer: ESA/ESTEC <<
30 // Revisions: Andrea Dotti, SLAC <<
31 // ------------------------------------------- <<
32 48
33 #include "G4SPSAngDistribution.hh" 49 #include "G4SPSAngDistribution.hh"
34 50
35 #include "Randomize.hh" 51 #include "Randomize.hh"
36 #include "G4PhysicalConstants.hh" 52 #include "G4PhysicalConstants.hh"
37 53
38 G4SPSAngDistribution::G4SPSAngDistribution() << 54
>> 55 G4SPSAngDistribution::G4SPSAngDistribution()
>> 56 : Theta(0.), Phi(0.), posDist(0),angRndm(0)
39 { 57 {
40 // Angular distribution Variables 58 // Angular distribution Variables
41 G4ThreeVector zero; 59 G4ThreeVector zero;
42 particle_momentum_direction = G4ParticleMome 60 particle_momentum_direction = G4ParticleMomentum(0,0,-1);
43 61
44 AngDistType = "planar"; 62 AngDistType = "planar";
45 AngRef1 = CLHEP::HepXHat; 63 AngRef1 = CLHEP::HepXHat;
46 AngRef2 = CLHEP::HepYHat; 64 AngRef2 = CLHEP::HepYHat;
47 AngRef3 = CLHEP::HepZHat; 65 AngRef3 = CLHEP::HepZHat;
48 MinTheta = 0.; 66 MinTheta = 0.;
49 MaxTheta = pi; 67 MaxTheta = pi;
50 MinPhi = 0.; 68 MinPhi = 0.;
51 MaxPhi = twopi; 69 MaxPhi = twopi;
52 DR = 0.; 70 DR = 0.;
53 DX = 0.; 71 DX = 0.;
54 DY = 0.; 72 DY = 0.;
55 FocusPoint = G4ThreeVector(0., 0., 0.); 73 FocusPoint = G4ThreeVector(0., 0., 0.);
56 UserDistType = "NULL"; 74 UserDistType = "NULL";
57 UserWRTSurface = true; 75 UserWRTSurface = true;
58 UserAngRef = false; 76 UserAngRef = false;
59 IPDFThetaExist = false; 77 IPDFThetaExist = false;
60 IPDFPhiExist = false; 78 IPDFPhiExist = false;
61 verbosityLevel = 0; << 79 verbosityLevel = 0 ;
62 << 80
63 G4MUTEXINIT(mutex); 81 G4MUTEXINIT(mutex);
64 } 82 }
65 83
66 G4SPSAngDistribution::~G4SPSAngDistribution() 84 G4SPSAngDistribution::~G4SPSAngDistribution()
67 { 85 {
68 G4MUTEXDESTROY(mutex); 86 G4MUTEXDESTROY(mutex);
69 } 87 }
70 88
71 void G4SPSAngDistribution::SetAngDistType(cons << 89 //
>> 90 void G4SPSAngDistribution::SetAngDistType(G4String atype)
72 { 91 {
73 G4AutoLock l(&mutex); << 92 G4AutoLock l(&mutex);
74 if(atype != "iso" && atype != "cos" && atype 93 if(atype != "iso" && atype != "cos" && atype != "user" && atype != "planar"
75 && atype != "beam1d" && atype != "beam2d" 94 && atype != "beam1d" && atype != "beam2d" && atype != "focused")
76 { << 95 G4cout << "Error, distribution must be iso, cos, planar, beam1d, beam2d, focused or user" << G4endl;
77 G4cout << "Error, distribution must be iso <<
78 << G4endl; <<
79 } <<
80 else 96 else
81 { <<
82 AngDistType = atype; 97 AngDistType = atype;
83 } << 98 if (AngDistType == "cos") MaxTheta = pi/2. ;
84 if (AngDistType == "cos") { MaxTheta = pi/2 << 99 if (AngDistType == "user") {
85 if (AngDistType == "user") << 100 UDefThetaH = IPDFThetaH = ZeroPhysVector ;
86 { << 101 IPDFThetaExist = false ;
87 UDefThetaH = IPDFThetaH = ZeroPhysVector; << 102 UDefPhiH = IPDFPhiH = ZeroPhysVector ;
88 IPDFThetaExist = false; << 103 IPDFPhiExist = false ;
89 UDefPhiH = IPDFPhiH = ZeroPhysVector; <<
90 IPDFPhiExist = false; <<
91 } 104 }
92 } 105 }
93 106
94 void G4SPSAngDistribution::DefineAngRefAxes(co << 107 void G4SPSAngDistribution::DefineAngRefAxes(G4String refname, G4ThreeVector ref)
95 co <<
96 { 108 {
97 G4AutoLock l(&mutex); << 109 G4AutoLock l(&mutex);
98 if (refname == "angref1") << 110 if(refname == "angref1")
99 AngRef1 = ref.unit(); // x' 111 AngRef1 = ref.unit(); // x'
100 else if (refname == "angref2") << 112 else if(refname == "angref2")
101 AngRef2 = ref.unit(); // vector in x'y' pl 113 AngRef2 = ref.unit(); // vector in x'y' plane
102 114
103 // User defines x' (AngRef1) and a vector in 115 // User defines x' (AngRef1) and a vector in the x'y'
104 // plane (AngRef2). Then, AngRef1 x AngRef2 116 // plane (AngRef2). Then, AngRef1 x AngRef2 = AngRef3
105 // the z' vector. Then, AngRef3 x AngRef1 = 117 // the z' vector. Then, AngRef3 x AngRef1 = AngRef2
106 // which will now be y'. 118 // which will now be y'.
107 119
108 AngRef3 = AngRef1.cross(AngRef2); // z' 120 AngRef3 = AngRef1.cross(AngRef2); // z'
109 AngRef2 = AngRef3.cross(AngRef1); // y' 121 AngRef2 = AngRef3.cross(AngRef1); // y'
110 UserAngRef = true ; 122 UserAngRef = true ;
111 if(verbosityLevel == 2) 123 if(verbosityLevel == 2)
112 { << 124 {
113 G4cout << "Angular distribution rotation a << 125 G4cout << "Angular distribution rotation axes " << AngRef1 << " " << AngRef2 << " " << AngRef3 << G4endl;
114 << " " << AngRef2 << " " << AngRef3 << 126 }
115 } <<
116 } 127 }
117 128
118 void G4SPSAngDistribution::SetMinTheta(G4doubl 129 void G4SPSAngDistribution::SetMinTheta(G4double mint)
119 { 130 {
120 G4AutoLock l(&mutex); << 131 G4AutoLock l(&mutex);
121 MinTheta = mint; 132 MinTheta = mint;
122 } 133 }
123 134
124 void G4SPSAngDistribution::SetMinPhi(G4double 135 void G4SPSAngDistribution::SetMinPhi(G4double minp)
125 { 136 {
126 G4AutoLock l(&mutex); << 137 G4AutoLock l(&mutex);
127 MinPhi = minp; 138 MinPhi = minp;
128 } 139 }
129 140
130 void G4SPSAngDistribution::SetMaxTheta(G4doubl 141 void G4SPSAngDistribution::SetMaxTheta(G4double maxt)
131 { 142 {
132 G4AutoLock l(&mutex); << 143 G4AutoLock l(&mutex);
133 MaxTheta = maxt; 144 MaxTheta = maxt;
134 } 145 }
135 146
136 void G4SPSAngDistribution::SetMaxPhi(G4double 147 void G4SPSAngDistribution::SetMaxPhi(G4double maxp)
137 { 148 {
138 G4AutoLock l(&mutex); << 149 G4AutoLock l(&mutex);
139 MaxPhi = maxp; 150 MaxPhi = maxp;
140 } 151 }
141 152
142 void G4SPSAngDistribution::SetBeamSigmaInAngR( 153 void G4SPSAngDistribution::SetBeamSigmaInAngR(G4double r)
143 { 154 {
144 G4AutoLock l(&mutex); << 155 G4AutoLock l(&mutex);
145 DR = r; 156 DR = r;
146 } 157 }
147 158
148 void G4SPSAngDistribution::SetBeamSigmaInAngX( 159 void G4SPSAngDistribution::SetBeamSigmaInAngX(G4double r)
149 { 160 {
150 G4AutoLock l(&mutex); << 161 G4AutoLock l(&mutex);
151 DX = r; 162 DX = r;
152 } 163 }
153 164
154 void G4SPSAngDistribution::SetBeamSigmaInAngY( 165 void G4SPSAngDistribution::SetBeamSigmaInAngY(G4double r)
155 { 166 {
156 G4AutoLock l(&mutex); << 167 G4AutoLock l(&mutex);
157 DY = r; 168 DY = r;
158 } 169 }
159 170
160 void G4SPSAngDistribution:: << 171 void G4SPSAngDistribution::SetParticleMomentumDirection(G4ParticleMomentum aMomentumDirection)
161 SetParticleMomentumDirection(const G4ParticleM <<
162 { 172 {
163 G4AutoLock l(&mutex); << 173 G4AutoLock l(&mutex);
164 particle_momentum_direction = aMomentumDirec << 174 particle_momentum_direction = aMomentumDirection.unit();
>> 175
165 } 176 }
166 177
167 void G4SPSAngDistribution::SetPosDistribution( 178 void G4SPSAngDistribution::SetPosDistribution(G4SPSPosDistribution* a)
168 { 179 {
169 G4AutoLock l(&mutex); << 180 G4AutoLock l(&mutex);
170 posDist = a; << 181 posDist = a;
171 } 182 }
172 183
173 void G4SPSAngDistribution::SetBiasRndm(G4SPSRa 184 void G4SPSAngDistribution::SetBiasRndm(G4SPSRandomGenerator* a)
174 { 185 {
175 G4AutoLock l(&mutex); << 186 G4AutoLock l(&mutex);
176 angRndm = a; << 187 angRndm = a;
177 } 188 }
178 189
179 void G4SPSAngDistribution::SetVerbosity(G4int 190 void G4SPSAngDistribution::SetVerbosity(G4int a)
180 { 191 {
181 G4AutoLock l(&mutex); << 192 G4AutoLock l(&mutex);
182 verbosityLevel = a; << 193 verbosityLevel = a;
183 } 194 }
184 195
185 void G4SPSAngDistribution::UserDefAngTheta(con << 196
>> 197 void G4SPSAngDistribution::UserDefAngTheta(G4ThreeVector input)
186 { 198 {
187 G4AutoLock l(&mutex); << 199 G4AutoLock l(&mutex);
188 if(UserDistType == "NULL") UserDistType = "t 200 if(UserDistType == "NULL") UserDistType = "theta";
189 if(UserDistType == "phi") UserDistType = "bo 201 if(UserDistType == "phi") UserDistType = "both";
190 G4double thi, val; 202 G4double thi, val;
191 thi = input.x(); 203 thi = input.x();
192 val = input.y(); 204 val = input.y();
193 if(verbosityLevel >= 1) G4cout << "In UserDe << 205 if(verbosityLevel >= 1)
>> 206 G4cout << "In UserDefAngTheta" << G4endl;
194 UDefThetaH.InsertValues(thi, val); 207 UDefThetaH.InsertValues(thi, val);
195 } 208 }
196 209
197 G4String G4SPSAngDistribution::GetDistType() << 210 G4String G4SPSAngDistribution::GetDistType() { G4AutoLock l(&mutex); return AngDistType;}
198 { << 211 G4double G4SPSAngDistribution::GetMinTheta() { G4AutoLock l(&mutex); return MinTheta; }
199 G4AutoLock l(&mutex); << 212 G4double G4SPSAngDistribution::GetMaxTheta() { G4AutoLock l(&mutex); return MaxTheta; }
200 return AngDistType; << 213 G4double G4SPSAngDistribution::GetMinPhi() { G4AutoLock l(&mutex); return MinPhi; }
201 } << 214 G4double G4SPSAngDistribution::GetMaxPhi() { G4AutoLock l(&mutex); return MaxPhi; }
202 << 215 G4ThreeVector G4SPSAngDistribution::GetDirection() { G4AutoLock l(&mutex); return particle_momentum_direction; }
203 G4double G4SPSAngDistribution::GetMinTheta() <<
204 { <<
205 G4AutoLock l(&mutex); <<
206 return MinTheta; <<
207 } <<
208 <<
209 G4double G4SPSAngDistribution::GetMaxTheta() <<
210 { <<
211 G4AutoLock l(&mutex); <<
212 return MaxTheta; <<
213 } <<
214 216
215 G4double G4SPSAngDistribution::GetMinPhi() << 217 void G4SPSAngDistribution::UserDefAngPhi(G4ThreeVector input)
216 { 218 {
217 G4AutoLock l(&mutex); << 219 G4AutoLock l(&mutex);
218 return MinPhi; <<
219 } <<
220 <<
221 G4double G4SPSAngDistribution::GetMaxPhi() <<
222 { <<
223 G4AutoLock l(&mutex); <<
224 return MaxPhi; <<
225 } <<
226 <<
227 G4ThreeVector G4SPSAngDistribution::GetDirecti <<
228 { <<
229 G4AutoLock l(&mutex); <<
230 return particle_momentum_direction; <<
231 } <<
232 <<
233 void G4SPSAngDistribution::UserDefAngPhi(const <<
234 { <<
235 G4AutoLock l(&mutex); <<
236 if(UserDistType == "NULL") UserDistType = "p 220 if(UserDistType == "NULL") UserDistType = "phi";
237 if(UserDistType == "theta") UserDistType = " 221 if(UserDistType == "theta") UserDistType = "both";
238 G4double phhi, val; 222 G4double phhi, val;
239 phhi = input.x(); 223 phhi = input.x();
240 val = input.y(); 224 val = input.y();
241 if(verbosityLevel >= 1) G4cout << "In UserDe << 225 if(verbosityLevel >= 1)
>> 226 G4cout << "In UserDefAngPhi" << G4endl;
242 UDefPhiH.InsertValues(phhi, val); 227 UDefPhiH.InsertValues(phhi, val);
243 } 228 }
244 229
245 void G4SPSAngDistribution::SetFocusPoint(const << 230 void G4SPSAngDistribution::SetFocusPoint(G4ThreeVector input)
246 { 231 {
247 G4AutoLock l(&mutex); << 232 G4AutoLock l(&mutex);
248 FocusPoint = input; 233 FocusPoint = input;
249 } 234 }
250 235
251 void G4SPSAngDistribution::SetUserWRTSurface(G 236 void G4SPSAngDistribution::SetUserWRTSurface(G4bool wrtSurf)
252 { 237 {
253 G4AutoLock l(&mutex); << 238 G4AutoLock l(&mutex);
254 << 239 // This is only applied in user mode?
255 // if UserWRTSurface = true then the user wa 240 // if UserWRTSurface = true then the user wants momenta with respect
256 // to the surface normals. 241 // to the surface normals.
257 // When doing this theta has to be 0-90 only 242 // When doing this theta has to be 0-90 only otherwise there will be
258 // errors, which currently are flagged anywh 243 // errors, which currently are flagged anywhere.
259 // <<
260 UserWRTSurface = wrtSurf; 244 UserWRTSurface = wrtSurf;
261 } 245 }
262 246
263 void G4SPSAngDistribution::SetUseUserAngAxis(G 247 void G4SPSAngDistribution::SetUseUserAngAxis(G4bool userang)
264 { 248 {
265 G4AutoLock l(&mutex); << 249 G4AutoLock l(&mutex);
266 <<
267 // if UserAngRef = true the angular distrib 250 // if UserAngRef = true the angular distribution is defined wrt
268 // the user defined coordinates << 251 // the user defined co-ordinates
269 // <<
270 UserAngRef = userang; 252 UserAngRef = userang;
271 } 253 }
272 254
273 void G4SPSAngDistribution::GenerateBeamFlux(G4 255 void G4SPSAngDistribution::GenerateBeamFlux(G4ParticleMomentum& mom)
274 { 256 {
275 G4double theta, phi; 257 G4double theta, phi;
276 G4double px, py, pz; 258 G4double px, py, pz;
277 if (AngDistType == "beam1d") 259 if (AngDistType == "beam1d")
278 { <<
279 theta = G4RandGauss::shoot(0.0,DR); <<
280 phi = twopi * G4UniformRand(); <<
281 } <<
282 else <<
283 { <<
284 px = G4RandGauss::shoot(0.0,DX); <<
285 py = G4RandGauss::shoot(0.0,DY); <<
286 theta = std::sqrt (px*px + py*py); <<
287 if (theta != 0.) <<
288 { 260 {
289 phi = std::acos(px/theta); << 261 theta = G4RandGauss::shoot(0.0,DR);
290 if ( py < 0.) phi = -phi; << 262 phi = twopi * G4UniformRand();
291 } 263 }
292 else << 264 else
293 { 265 {
294 phi = 0.0; << 266 px = G4RandGauss::shoot(0.0,DX);
>> 267 py = G4RandGauss::shoot(0.0,DY);
>> 268 theta = std::sqrt (px*px + py*py);
>> 269 if (theta != 0.) {
>> 270 phi = std::acos(px/theta);
>> 271 if ( py < 0.) phi = -phi;
>> 272 } else {
>> 273 phi = 0.0;
>> 274 }
295 } 275 }
296 } <<
297 px = -std::sin(theta) * std::cos(phi); 276 px = -std::sin(theta) * std::cos(phi);
298 py = -std::sin(theta) * std::sin(phi); 277 py = -std::sin(theta) * std::sin(phi);
299 pz = -std::cos(theta); 278 pz = -std::cos(theta);
300 G4double finx, finy, finz; << 279 G4double finx, finy, finz ;
301 finx=px, finy=py, finz=pz; << 280 finx = px, finy =py, finz =pz;
302 if (UserAngRef) << 281 if (UserAngRef){
303 { <<
304 // Apply Angular Rotation Matrix 282 // Apply Angular Rotation Matrix
305 // x * AngRef1, y * AngRef2 and z * AngRef 283 // x * AngRef1, y * AngRef2 and z * AngRef3
306 finx = (px * AngRef1.x()) + (py * AngRef2. 284 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
307 finy = (px * AngRef1.y()) + (py * AngRef2. 285 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
308 finz = (px * AngRef1.z()) + (py * AngRef2. 286 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
309 G4double ResMag = std::sqrt((finx*finx) + 287 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
310 finx = finx/ResMag; 288 finx = finx/ResMag;
311 finy = finy/ResMag; 289 finy = finy/ResMag;
312 finz = finz/ResMag; 290 finz = finz/ResMag;
313 } 291 }
314 mom.setX(finx); 292 mom.setX(finx);
315 mom.setY(finy); 293 mom.setY(finy);
316 mom.setZ(finz); 294 mom.setZ(finz);
317 295
318 // particle_momentum_direction now holds uni << 296 // particle_momentum_direction now holds unit momentum vector.
319 <<
320 if(verbosityLevel >= 1) 297 if(verbosityLevel >= 1)
321 { <<
322 G4cout << "Generating beam vector: " << mo 298 G4cout << "Generating beam vector: " << mom << G4endl;
323 } <<
324 } 299 }
325 300
326 void G4SPSAngDistribution::GenerateFocusedFlux 301 void G4SPSAngDistribution::GenerateFocusedFlux(G4ParticleMomentum& mom)
327 { 302 {
328 mom = (FocusPoint - posDist->GetParticlePos( 303 mom = (FocusPoint - posDist->GetParticlePos()).unit();
329 << 304 //
330 // particle_momentum_direction now holds uni 305 // particle_momentum_direction now holds unit momentum vector.
331 <<
332 if(verbosityLevel >= 1) 306 if(verbosityLevel >= 1)
333 { <<
334 G4cout << "Generating focused vector: " << 307 G4cout << "Generating focused vector: " << mom << G4endl;
335 } <<
336 } 308 }
337 309
338 void G4SPSAngDistribution::GenerateIsotropicFl 310 void G4SPSAngDistribution::GenerateIsotropicFlux(G4ParticleMomentum& mom)
339 { 311 {
340 // generates isotropic flux. 312 // generates isotropic flux.
341 // No vectors are needed. 313 // No vectors are needed.
342 <<
343 G4double rndm, rndm2; 314 G4double rndm, rndm2;
344 G4double px, py, pz; 315 G4double px, py, pz;
345 316
>> 317 //
346 G4double sintheta, sinphi,costheta,cosphi; 318 G4double sintheta, sinphi,costheta,cosphi;
347 rndm = angRndm->GenRandTheta(); 319 rndm = angRndm->GenRandTheta();
348 costheta = std::cos(MinTheta) - rndm * (std: << 320 costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) - std::cos(MaxTheta));
349 - std: <<
350 sintheta = std::sqrt(1. - costheta*costheta) 321 sintheta = std::sqrt(1. - costheta*costheta);
351 322
352 rndm2 = angRndm->GenRandPhi(); 323 rndm2 = angRndm->GenRandPhi();
353 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 324 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
354 sinphi = std::sin(Phi); 325 sinphi = std::sin(Phi);
355 cosphi = std::cos(Phi); 326 cosphi = std::cos(Phi);
356 327
357 px = -sintheta * cosphi; 328 px = -sintheta * cosphi;
358 py = -sintheta * sinphi; 329 py = -sintheta * sinphi;
359 pz = -costheta; 330 pz = -costheta;
360 331
361 // For volume and point source use mother or << 332 // for volume and ponit source use mother or user defined co-ordinates
362 // for plane and surface source user surface << 333 // for plane and surface source user surface-normal or userdefined co-ordinates
363 // coordinates <<
364 // 334 //
365 G4double finx, finy, finz; 335 G4double finx, finy, finz;
366 if (posDist->GetSourcePosType() == "Point" << 336 if (posDist->GetSourcePosType() == "Point" || posDist->GetSourcePosType() == "Volume") {
367 || posDist->GetSourcePosType() == "Volume") << 337 if (UserAngRef){
368 { <<
369 if (UserAngRef) <<
370 { <<
371 // Apply Rotation Matrix 338 // Apply Rotation Matrix
372 // x * AngRef1, y * AngRef2 and z * AngR 339 // x * AngRef1, y * AngRef2 and z * AngRef3
373 finx = (px * AngRef1.x()) + (py * AngRef 340 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
374 finy = (px * AngRef1.y()) + (py * AngRef 341 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
375 finz = (px * AngRef1.z()) + (py * AngRef 342 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
376 } << 343 } else {
377 else <<
378 { <<
379 finx = px; 344 finx = px;
380 finy = py; 345 finy = py;
381 finz = pz; 346 finz = pz;
382 } 347 }
383 } << 348 } else { // for plane and surface source
384 else << 349 if (UserAngRef){
385 { // for plane and surface source <<
386 if (UserAngRef) <<
387 { <<
388 // Apply Rotation Matrix 350 // Apply Rotation Matrix
389 // x * AngRef1, y * AngRef2 and z * AngR 351 // x * AngRef1, y * AngRef2 and z * AngRef3
390 finx = (px * AngRef1.x()) + (py * AngRef 352 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
391 finy = (px * AngRef1.y()) + (py * AngRef 353 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
392 finz = (px * AngRef1.z()) + (py * AngRef 354 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
393 } << 355 } else {
394 else << 356 finx = (px*posDist->GetSideRefVec1().x()) + (py*posDist->GetSideRefVec2().x()) + (pz*posDist->GetSideRefVec3().x());
395 { << 357 finy = (px*posDist->GetSideRefVec1().y()) + (py*posDist->GetSideRefVec2().y()) + (pz*posDist->GetSideRefVec3().y());
396 finx = (px*posDist->GetSideRefVec1().x() << 358 finz = (px*posDist->GetSideRefVec1().z()) + (py*posDist->GetSideRefVec2().z()) + (pz*posDist->GetSideRefVec3().z());
397 + (py*posDist->GetSideRefVec2().x() <<
398 + (pz*posDist->GetSideRefVec3().x() <<
399 finy = (px*posDist->GetSideRefVec1().y() <<
400 + (py*posDist->GetSideRefVec2().y() <<
401 + (pz*posDist->GetSideRefVec3().y() <<
402 finz = (px*posDist->GetSideRefVec1().z() <<
403 + (py*posDist->GetSideRefVec2().z() <<
404 + (pz*posDist->GetSideRefVec3().z() <<
405 } 359 }
406 } 360 }
407 G4double ResMag = std::sqrt((finx*finx) + (f 361 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
408 finx = finx/ResMag; 362 finx = finx/ResMag;
409 finy = finy/ResMag; 363 finy = finy/ResMag;
410 finz = finz/ResMag; 364 finz = finz/ResMag;
411 365
412 mom.setX(finx); 366 mom.setX(finx);
413 mom.setY(finy); 367 mom.setY(finy);
414 mom.setZ(finz); 368 mom.setZ(finz);
415 369
416 // particle_momentum_direction now holds uni 370 // particle_momentum_direction now holds unit momentum vector.
417 <<
418 if(verbosityLevel >= 1) 371 if(verbosityLevel >= 1)
419 { <<
420 G4cout << "Generating isotropic vector: " 372 G4cout << "Generating isotropic vector: " << mom << G4endl;
421 } <<
422 } 373 }
423 374
424 void G4SPSAngDistribution::GenerateCosineLawFl 375 void G4SPSAngDistribution::GenerateCosineLawFlux(G4ParticleMomentum& mom)
425 { 376 {
426 // Method to generate flux distributed with 377 // Method to generate flux distributed with a cosine law
427 <<
428 G4double px, py, pz; 378 G4double px, py, pz;
429 G4double rndm, rndm2; 379 G4double rndm, rndm2;
430 << 380 //
431 G4double sintheta, sinphi,costheta,cosphi; 381 G4double sintheta, sinphi,costheta,cosphi;
432 rndm = angRndm->GenRandTheta(); 382 rndm = angRndm->GenRandTheta();
433 sintheta = std::sqrt( rndm * (std::sin(MaxTh << 383 sintheta = std::sqrt( rndm * (std::sin(MaxTheta)*std::sin(MaxTheta) - std::sin(MinTheta)*std::sin(MinTheta) )
434 - std::sin(MinTh << 384 +std::sin(MinTheta)*std::sin(MinTheta) );
435 + std::sin(MinTheta)*std <<
436 costheta = std::sqrt(1. -sintheta*sintheta); 385 costheta = std::sqrt(1. -sintheta*sintheta);
437 386
438 rndm2 = angRndm->GenRandPhi(); 387 rndm2 = angRndm->GenRandPhi();
439 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 388 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
440 sinphi = std::sin(Phi); 389 sinphi = std::sin(Phi);
441 cosphi = std::cos(Phi); 390 cosphi = std::cos(Phi);
442 391
443 px = -sintheta * cosphi; 392 px = -sintheta * cosphi;
444 py = -sintheta * sinphi; 393 py = -sintheta * sinphi;
445 pz = -costheta; 394 pz = -costheta;
446 395
447 // for volume and point source use mother or << 396 // for volume and ponit source use mother or user defined co-ordinates
448 // for plane and surface source user surface << 397 // for plane and surface source user surface-normal or userdefined co-ordinates
449 // coordinates <<
450 // 398 //
451 G4double finx, finy, finz; 399 G4double finx, finy, finz;
452 if (posDist->GetSourcePosType() == "Point" << 400 if (posDist->GetSourcePosType() == "Point" || posDist->GetSourcePosType() == "Volume") {
453 || posDist->GetSourcePosType() == "Volume") << 401 if (UserAngRef){
454 { <<
455 if (UserAngRef) <<
456 { <<
457 // Apply Rotation Matrix 402 // Apply Rotation Matrix
458 finx = (px * AngRef1.x()) + (py * AngRef 403 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
459 finy = (px * AngRef1.y()) + (py * AngRef 404 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
460 finz = (px * AngRef1.z()) + (py * AngRef 405 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
461 } << 406 } else {
462 else <<
463 { <<
464 finx = px; 407 finx = px;
465 finy = py; 408 finy = py;
466 finz = pz; 409 finz = pz;
467 } 410 }
468 } << 411 } else { // for plane and surface source
469 else << 412 if (UserAngRef){
470 { // for plane and surface source <<
471 if (UserAngRef) <<
472 { <<
473 // Apply Rotation Matrix 413 // Apply Rotation Matrix
474 finx = (px * AngRef1.x()) + (py * AngRef 414 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
475 finy = (px * AngRef1.y()) + (py * AngRef 415 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
476 finz = (px * AngRef1.z()) + (py * AngRef 416 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
477 } << 417 } else {
478 else << 418 finx = (px*posDist->GetSideRefVec1().x()) + (py*posDist->GetSideRefVec2().x()) + (pz*posDist->GetSideRefVec3().x());
479 { << 419 finy = (px*posDist->GetSideRefVec1().y()) + (py*posDist->GetSideRefVec2().y()) + (pz*posDist->GetSideRefVec3().y());
480 finx = (px*posDist->GetSideRefVec1().x() << 420 finz = (px*posDist->GetSideRefVec1().z()) + (py*posDist->GetSideRefVec2().z()) + (pz*posDist->GetSideRefVec3().z());
481 + (py*posDist->GetSideRefVec2().x() <<
482 + (pz*posDist->GetSideRefVec3().x() <<
483 finy = (px*posDist->GetSideRefVec1().y() <<
484 + (py*posDist->GetSideRefVec2().y() <<
485 + (pz*posDist->GetSideRefVec3().y() <<
486 finz = (px*posDist->GetSideRefVec1().z() <<
487 + (py*posDist->GetSideRefVec2().z() <<
488 + (pz*posDist->GetSideRefVec3().z() <<
489 } 421 }
490 } 422 }
491 G4double ResMag = std::sqrt((finx*finx) + (f 423 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
492 finx = finx/ResMag; 424 finx = finx/ResMag;
493 finy = finy/ResMag; 425 finy = finy/ResMag;
494 finz = finz/ResMag; 426 finz = finz/ResMag;
495 427
496 mom.setX(finx); 428 mom.setX(finx);
497 mom.setY(finy); 429 mom.setY(finy);
498 mom.setZ(finz); 430 mom.setZ(finz);
499 431
500 // particle_momentum_direction now contains 432 // particle_momentum_direction now contains unit momentum vector.
501 <<
502 if(verbosityLevel >= 1) 433 if(verbosityLevel >= 1)
503 { << 434 {
504 G4cout << "Resultant cosine-law unit momen << 435 G4cout << "Resultant cosine-law unit momentum vector " << mom << G4endl;
505 } << 436 }
506 } 437 }
507 438
508 void G4SPSAngDistribution::GeneratePlanarFlux( 439 void G4SPSAngDistribution::GeneratePlanarFlux(G4ParticleMomentum& mom)
509 { 440 {
510 // particle_momentum_direction now contains 441 // particle_momentum_direction now contains unit momentum vector.
511 // nothing need be done here as the m-direct 442 // nothing need be done here as the m-directions have been set directly
512 // under this option 443 // under this option
513 <<
514 if(verbosityLevel >= 1) 444 if(verbosityLevel >= 1)
515 { << 445 {
516 G4cout << "Resultant Planar wave momentum << 446 G4cout << "Resultant Planar wave momentum vector " << mom << G4endl;
517 } << 447 }
518 } 448 }
519 449
520 void G4SPSAngDistribution::GenerateUserDefFlux 450 void G4SPSAngDistribution::GenerateUserDefFlux(G4ParticleMomentum& mom)
521 { 451 {
522 G4double rndm, px, py, pz, pmag; 452 G4double rndm, px, py, pz, pmag;
523 453
524 if(UserDistType == "NULL") 454 if(UserDistType == "NULL")
525 { <<
526 G4cout << "Error: UserDistType undefined" 455 G4cout << "Error: UserDistType undefined" << G4endl;
527 } << 456 else if(UserDistType == "theta") {
528 else if(UserDistType == "theta") <<
529 { <<
530 Theta = 10.; 457 Theta = 10.;
531 while(Theta > MaxTheta || Theta < MinTheta 458 while(Theta > MaxTheta || Theta < MinTheta)
532 { <<
533 Theta = GenerateUserDefTheta(); 459 Theta = GenerateUserDefTheta();
534 } <<
535 Phi = 10.; 460 Phi = 10.;
536 while(Phi > MaxPhi || Phi < MinPhi) << 461 while(Phi > MaxPhi || Phi < MinPhi) {
537 { <<
538 rndm = angRndm->GenRandPhi(); 462 rndm = angRndm->GenRandPhi();
539 Phi = twopi * rndm; 463 Phi = twopi * rndm;
540 } 464 }
541 } 465 }
542 else if(UserDistType == "phi") << 466 else if(UserDistType == "phi") {
543 { <<
544 Theta = 10.; 467 Theta = 10.;
545 while(Theta > MaxTheta || Theta < MinTheta 468 while(Theta > MaxTheta || Theta < MinTheta)
546 { << 469 {
547 rndm = angRndm->GenRandTheta(); << 470 rndm = angRndm->GenRandTheta();
548 Theta = std::acos(1. - (2. * rndm)); << 471 Theta = std::acos(1. - (2. * rndm));
549 } << 472 }
550 Phi = 10.; 473 Phi = 10.;
551 while(Phi > MaxPhi || Phi < MinPhi) 474 while(Phi > MaxPhi || Phi < MinPhi)
552 { <<
553 Phi = GenerateUserDefPhi(); 475 Phi = GenerateUserDefPhi();
554 } <<
555 } 476 }
556 else if(UserDistType == "both") 477 else if(UserDistType == "both")
557 { <<
558 Theta = 10.; <<
559 while(Theta > MaxTheta || Theta < MinTheta <<
560 { <<
561 Theta = GenerateUserDefTheta(); <<
562 } <<
563 Phi = 10.; <<
564 while(Phi > MaxPhi || Phi < MinPhi) <<
565 { 478 {
566 Phi = GenerateUserDefPhi(); << 479 Theta = 10.;
>> 480 while(Theta > MaxTheta || Theta < MinTheta)
>> 481 Theta = GenerateUserDefTheta();
>> 482 Phi = 10.;
>> 483 while(Phi > MaxPhi || Phi < MinPhi)
>> 484 Phi = GenerateUserDefPhi();
567 } 485 }
568 } <<
569 px = -std::sin(Theta) * std::cos(Phi); 486 px = -std::sin(Theta) * std::cos(Phi);
570 py = -std::sin(Theta) * std::sin(Phi); 487 py = -std::sin(Theta) * std::sin(Phi);
571 pz = -std::cos(Theta); 488 pz = -std::cos(Theta);
572 489
573 pmag = std::sqrt((px*px) + (py*py) + (pz*pz) 490 pmag = std::sqrt((px*px) + (py*py) + (pz*pz));
574 491
575 if(!UserWRTSurface) << 492 if(!UserWRTSurface) {
576 { <<
577 G4double finx, finy, finz; 493 G4double finx, finy, finz;
578 if (UserAngRef) << 494 if (UserAngRef) {
579 { <<
580 // Apply Rotation Matrix 495 // Apply Rotation Matrix
581 // x * AngRef1, y * AngRef2 and z * AngR 496 // x * AngRef1, y * AngRef2 and z * AngRef3
582 finx = (px * AngRef1.x()) + (py * AngRef 497 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
583 finy = (px * AngRef1.y()) + (py * AngRef 498 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
584 finz = (px * AngRef1.z()) + (py * AngRef 499 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
585 } << 500 } else { // use mother co-ordinates
586 else // use mother coordinates <<
587 { <<
588 finx = px; 501 finx = px;
589 finy = py; 502 finy = py;
590 finz = pz; 503 finz = pz;
591 } 504 }
592 G4double ResMag = std::sqrt((finx*finx) + 505 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
593 finx = finx/ResMag; 506 finx = finx/ResMag;
594 finy = finy/ResMag; 507 finy = finy/ResMag;
595 finz = finz/ResMag; 508 finz = finz/ResMag;
596 509
597 mom.setX(finx); 510 mom.setX(finx);
598 mom.setY(finy); 511 mom.setY(finy);
599 mom.setZ(finz); 512 mom.setZ(finz);
600 } 513 }
601 else // UserWRTSurface = true << 514 else { // UserWRTSurface = true
602 { <<
603 G4double pxh = px/pmag; 515 G4double pxh = px/pmag;
604 G4double pyh = py/pmag; 516 G4double pyh = py/pmag;
605 G4double pzh = pz/pmag; 517 G4double pzh = pz/pmag;
606 if(verbosityLevel > 1) << 518 if(verbosityLevel > 1) {
607 { << 519 G4cout <<"SideRefVecs " <<posDist->GetSideRefVec1()<<posDist->GetSideRefVec2()<<posDist->GetSideRefVec3()<<G4endl;
608 G4cout << "SideRefVecs " << posDist->Get << 520 G4cout <<"Raw Unit vector "<<pxh<<","<<pyh<<","<<pzh<<G4endl;
609 << posDist->GetSideRefVec2() << p << 521 }
610 << G4endl; << 522 G4double resultx = (pxh*posDist->GetSideRefVec1().x()) + (pyh*posDist->GetSideRefVec2().x()) +
611 G4cout << "Raw Unit vector " << pxh << 523 (pzh*posDist->GetSideRefVec3().x());
612 << "," << pyh << "," << pzh << G4 <<
613 } <<
614 G4double resultx = (pxh*posDist->GetSideRe <<
615 + (pyh*posDist->GetSideRe <<
616 + (pzh*posDist->GetSideRe <<
617 524
618 G4double resulty = (pxh*posDist->GetSideRe << 525 G4double resulty = (pxh*posDist->GetSideRefVec1().y()) + (pyh*posDist->GetSideRefVec2().y()) +
619 + (pyh*posDist->GetSideRe << 526 (pzh*posDist->GetSideRefVec3().y());
620 + (pzh*posDist->GetSideRe <<
621 527
622 G4double resultz = (pxh*posDist->GetSideRe << 528 G4double resultz = (pxh*posDist->GetSideRefVec1().z()) + (pyh*posDist->GetSideRefVec2().z()) +
623 + (pyh*posDist->GetSideRe << 529 (pzh*posDist->GetSideRefVec3().z());
624 + (pzh*posDist->GetSideRe <<
625 530
626 G4double ResMag = std::sqrt((resultx*resul << 531 G4double ResMag = std::sqrt((resultx*resultx) + (resulty*resulty) + (resultz*resultz));
627 + (resulty*resul <<
628 + (resultz*resul <<
629 resultx = resultx/ResMag; 532 resultx = resultx/ResMag;
630 resulty = resulty/ResMag; 533 resulty = resulty/ResMag;
631 resultz = resultz/ResMag; 534 resultz = resultz/ResMag;
632 535
633 mom.setX(resultx); 536 mom.setX(resultx);
634 mom.setY(resulty); 537 mom.setY(resulty);
635 mom.setZ(resultz); 538 mom.setZ(resultz);
636 } 539 }
637 540
638 // particle_momentum_direction now contains 541 // particle_momentum_direction now contains unit momentum vector.
639 <<
640 if(verbosityLevel > 0 ) 542 if(verbosityLevel > 0 )
641 { << 543 {
642 G4cout << "Final User Defined momentum vec << 544 G4cout << "Final User Defined momentum vector " << particle_momentum_direction << G4endl;
643 << particle_momentum_direction << G << 545 }
644 } <<
645 } 546 }
646 547
647 G4double G4SPSAngDistribution::GenerateUserDef 548 G4double G4SPSAngDistribution::GenerateUserDefTheta()
648 { 549 {
649 // Create cumulative histogram if not alread << 550 // Create cumulative histogram if not already done so. Then use RandFlat
650 // Then use RandFlat::shoot to generate the << 551 //::shoot to generate the output Theta value.
651 <<
652 if(UserDistType == "NULL" || UserDistType == 552 if(UserDistType == "NULL" || UserDistType == "phi")
653 { <<
654 // No user defined theta distribution <<
655 G4cout << "Error ***********************" <<
656 G4cout << "UserDistType = " << UserDistTyp <<
657 return (0.); <<
658 } <<
659 <<
660 // UserDistType = theta or both and so a the <<
661 // is defined. This should be integrated if <<
662 G4AutoLock l(&mutex); <<
663 if(!IPDFThetaExist) <<
664 { <<
665 // IPDF has not been created, so create it <<
666 // <<
667 G4double bins[1024],vals[1024], sum; <<
668 G4int ii; <<
669 auto maxbin = G4int(UDefThetaH.GetVectorL <<
670 bins[0] = UDefThetaH.GetLowEdgeEnergy(std: <<
671 vals[0] = UDefThetaH(std::size_t(0)); <<
672 sum = vals[0]; <<
673 for(ii=1; ii<maxbin; ++ii) <<
674 { 553 {
675 bins[ii] = UDefThetaH.GetLowEdgeEnergy(s << 554 // No user defined theta distribution
676 vals[ii] = UDefThetaH(std::size_t(ii)) + << 555 G4cout << "Error ***********************" << G4endl;
677 sum = sum + UDefThetaH(std::size_t(ii)); << 556 G4cout << "UserDistType = " << UserDistType << G4endl;
>> 557 return (0.);
678 } 558 }
679 for(ii=0; ii<maxbin; ++ii) << 559 else
680 { 560 {
681 vals[ii] = vals[ii]/sum; << 561 // UserDistType = theta or both and so a theta distribution
682 IPDFThetaH.InsertValues(bins[ii], vals[i << 562 // is defined. This should be integrated if not already done.
>> 563 G4AutoLock l(&mutex);
>> 564 if(IPDFThetaExist == false)
>> 565 {
>> 566 // IPDF has not been created, so create it
>> 567 G4double bins[1024],vals[1024], sum;
>> 568 G4int ii;
>> 569 G4int maxbin = G4int(UDefThetaH.GetVectorLength());
>> 570 bins[0] = UDefThetaH.GetLowEdgeEnergy(size_t(0));
>> 571 vals[0] = UDefThetaH(size_t(0));
>> 572 sum = vals[0];
>> 573 for(ii=1;ii<maxbin;ii++)
>> 574 {
>> 575 bins[ii] = UDefThetaH.GetLowEdgeEnergy(size_t(ii));
>> 576 vals[ii] = UDefThetaH(size_t(ii)) + vals[ii-1];
>> 577 sum = sum + UDefThetaH(size_t(ii));
>> 578 }
>> 579 for(ii=0;ii<maxbin;ii++)
>> 580 {
>> 581 vals[ii] = vals[ii]/sum;
>> 582 IPDFThetaH.InsertValues(bins[ii], vals[ii]);
>> 583 }
>> 584 // Make IPDFThetaExist = true
>> 585 IPDFThetaExist = true;
>> 586 }
>> 587 l.unlock();
>> 588 // IPDF has been create so carry on
>> 589 G4double rndm = G4UniformRand();
>> 590 return(IPDFThetaH.GetEnergy(rndm));
683 } 591 }
684 IPDFThetaExist = true; <<
685 } <<
686 l.unlock(); <<
687 <<
688 // IPDF has been created so carry on <<
689 // <<
690 G4double rndm = G4UniformRand(); <<
691 return IPDFThetaH.GetEnergy(rndm); <<
692 } 592 }
693 593
694 G4double G4SPSAngDistribution::GenerateUserDef 594 G4double G4SPSAngDistribution::GenerateUserDefPhi()
695 { 595 {
696 // Create cumulative histogram if not alread << 596 // Create cumulative histogram if not already done so. Then use RandFlat
697 // Then use RandFlat::shoot to generate the << 597 //::shoot to generate the output Theta value.
698 598
699 if(UserDistType == "NULL" || UserDistType == 599 if(UserDistType == "NULL" || UserDistType == "theta")
700 { <<
701 // No user defined phi distribution <<
702 G4cout << "Error ***********************" <<
703 G4cout << "UserDistType = " << UserDistTyp <<
704 return(0.); <<
705 } <<
706 <<
707 // UserDistType = phi or both and so a phi d <<
708 // is defined. This should be integrated if <<
709 G4AutoLock l(&mutex); <<
710 if(!IPDFPhiExist) <<
711 { <<
712 // IPDF has not been created, so create it <<
713 // <<
714 G4double bins[1024],vals[1024], sum; <<
715 G4int ii; <<
716 auto maxbin = G4int(UDefPhiH.GetVectorLen <<
717 bins[0] = UDefPhiH.GetLowEdgeEnergy(std::s <<
718 vals[0] = UDefPhiH(std::size_t(0)); <<
719 sum = vals[0]; <<
720 for(ii=1; ii<maxbin; ++ii) <<
721 { 600 {
722 bins[ii] = UDefPhiH.GetLowEdgeEnergy(std << 601 // No user defined phi distribution
723 vals[ii] = UDefPhiH(std::size_t(ii)) + v << 602 G4cout << "Error ***********************" << G4endl;
724 sum = sum + UDefPhiH(std::size_t(ii)); << 603 G4cout << "UserDistType = " << UserDistType << G4endl;
>> 604 return(0.);
725 } 605 }
726 for(ii=0; ii<maxbin; ++ii) << 606 else
727 { 607 {
728 vals[ii] = vals[ii]/sum; << 608 // UserDistType = phi or both and so a phi distribution
729 IPDFPhiH.InsertValues(bins[ii], vals[ii] << 609 // is defined. This should be integrated if not already done.
>> 610 G4AutoLock l(&mutex);
>> 611 if(IPDFPhiExist == false)
>> 612 {
>> 613 // IPDF has not been created, so create it
>> 614 G4double bins[1024],vals[1024], sum;
>> 615 G4int ii;
>> 616 G4int maxbin = G4int(UDefPhiH.GetVectorLength());
>> 617 bins[0] = UDefPhiH.GetLowEdgeEnergy(size_t(0));
>> 618 vals[0] = UDefPhiH(size_t(0));
>> 619 sum = vals[0];
>> 620 for(ii=1;ii<maxbin;ii++)
>> 621 {
>> 622 bins[ii] = UDefPhiH.GetLowEdgeEnergy(size_t(ii));
>> 623 vals[ii] = UDefPhiH(size_t(ii)) + vals[ii-1];
>> 624 sum = sum + UDefPhiH(size_t(ii));
>> 625 }
>> 626 for(ii=0;ii<maxbin;ii++)
>> 627 {
>> 628 vals[ii] = vals[ii]/sum;
>> 629 IPDFPhiH.InsertValues(bins[ii], vals[ii]);
>> 630 }
>> 631 // Make IPDFPhiExist = true
>> 632 IPDFPhiExist = true;
>> 633 }
>> 634 l.unlock();
>> 635 // IPDF has been create so carry on
>> 636 G4double rndm = G4UniformRand();
>> 637 return(IPDFPhiH.GetEnergy(rndm));
730 } 638 }
731 IPDFPhiExist = true; <<
732 } <<
733 l.unlock(); <<
734 <<
735 // IPDF has been create so carry on <<
736 // <<
737 G4double rndm = G4UniformRand(); <<
738 return IPDFPhiH.GetEnergy(rndm); <<
739 } 639 }
740 << 640 //
741 void G4SPSAngDistribution::ReSetHist(const G4S << 641 void G4SPSAngDistribution::ReSetHist(G4String atype)
742 { 642 {
743 G4AutoLock l(&mutex); << 643 G4AutoLock l(&mutex);
744 if (atype == "theta") << 644 if (atype == "theta") {
745 { <<
746 UDefThetaH = IPDFThetaH = ZeroPhysVector ; 645 UDefThetaH = IPDFThetaH = ZeroPhysVector ;
747 IPDFThetaExist = false ; << 646 IPDFThetaExist = false ;}
748 } << 647 else if (atype == "phi"){
749 else if (atype == "phi") <<
750 { <<
751 UDefPhiH = IPDFPhiH = ZeroPhysVector ; 648 UDefPhiH = IPDFPhiH = ZeroPhysVector ;
752 IPDFPhiExist = false ; << 649 IPDFPhiExist = false ;}
753 } << 650 else {
754 else <<
755 { <<
756 G4cout << "Error, histtype not accepted " 651 G4cout << "Error, histtype not accepted " << G4endl;
757 } 652 }
758 } 653 }
759 654
>> 655
760 G4ParticleMomentum G4SPSAngDistribution::Gener 656 G4ParticleMomentum G4SPSAngDistribution::GenerateOne()
761 { 657 {
762 // Local copy for thread safety << 658 //Local copy for thread safety
763 // <<
764 G4ParticleMomentum localM = particle_momentu 659 G4ParticleMomentum localM = particle_momentum_direction;
765 <<
766 // Angular stuff 660 // Angular stuff
767 // <<
768 if(AngDistType == "iso") 661 if(AngDistType == "iso")
769 GenerateIsotropicFlux(localM); 662 GenerateIsotropicFlux(localM);
770 else if(AngDistType == "cos") 663 else if(AngDistType == "cos")
771 GenerateCosineLawFlux(localM); 664 GenerateCosineLawFlux(localM);
772 else if(AngDistType == "planar") 665 else if(AngDistType == "planar")
773 GeneratePlanarFlux(localM); 666 GeneratePlanarFlux(localM);
774 else if(AngDistType == "beam1d" || AngDistTy 667 else if(AngDistType == "beam1d" || AngDistType == "beam2d" )
775 GenerateBeamFlux(localM); 668 GenerateBeamFlux(localM);
776 else if(AngDistType == "user") 669 else if(AngDistType == "user")
777 GenerateUserDefFlux(localM); 670 GenerateUserDefFlux(localM);
778 else if(AngDistType == "focused") 671 else if(AngDistType == "focused")
779 GenerateFocusedFlux(localM); 672 GenerateFocusedFlux(localM);
780 else 673 else
781 G4cout << "Error: AngDistType has unusual 674 G4cout << "Error: AngDistType has unusual value" << G4endl;
782 return localM; 675 return localM;
783 } 676 }
>> 677
>> 678
>> 679
>> 680
>> 681
>> 682
>> 683
>> 684
>> 685
784 686