| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4SPSPosDistribution class implementation << 23 ///////////////////////////////////////////////////////////////////////////////
>> 24 //
>> 25 // MODULE: G4SPSPosDistribution.cc
>> 26 //
>> 27 // Version: 1.0
>> 28 // Date: 5/02/04
>> 29 // Author: Fan Lei
>> 30 // Organisation: QinetiQ ltd.
>> 31 // Customer: ESA/ESTEC
>> 32 //
>> 33 ///////////////////////////////////////////////////////////////////////////////
>> 34 //
>> 35 // CHANGE HISTORY
>> 36 // --------------
>> 37 //
>> 38 //
>> 39 // Version 1.0, 05/02/2004, Fan Lei, Created.
>> 40 // Based on the G4GeneralParticleSource class in Geant4 v6.0
>> 41 //
>> 42 ///////////////////////////////////////////////////////////////////////////////
27 // 43 //
28 // Author: Fan Lei, QinetiQ ltd. - 05/02/2004 <<
29 // Customer: ESA/ESTEC <<
30 // Revisions: Andrea Dotti, John Allison, Mako <<
31 // ------------------------------------------- <<
32 <<
33 #include "G4SPSPosDistribution.hh" <<
34 <<
35 #include "G4PhysicalConstants.hh" <<
36 #include "Randomize.hh" 44 #include "Randomize.hh"
>> 45 //#include <math.h>
37 #include "G4TransportationManager.hh" 46 #include "G4TransportationManager.hh"
38 #include "G4VPhysicalVolume.hh" 47 #include "G4VPhysicalVolume.hh"
39 #include "G4PhysicalVolumeStore.hh" 48 #include "G4PhysicalVolumeStore.hh"
40 #include "G4AutoLock.hh" <<
41 #include "G4AutoDelete.hh" <<
42 49
43 G4SPSPosDistribution::thread_data_t::thread_da << 50 #include "G4SPSPosDistribution.hh"
44 { <<
45 CSideRefVec1 = G4ThreeVector(CLHEP::HepXHat) <<
46 CSideRefVec2 = G4ThreeVector(CLHEP::HepYHat) <<
47 CSideRefVec3 = G4ThreeVector(CLHEP::HepZHat) <<
48 CParticlePos = G4ThreeVector(0,0,0); <<
49 } <<
50 51
51 G4SPSPosDistribution::G4SPSPosDistribution() 52 G4SPSPosDistribution::G4SPSPosDistribution()
52 { 53 {
>> 54
>> 55 // Initialise all variables
>> 56 // Position distribution Variables
>> 57
53 SourcePosType = "Point"; 58 SourcePosType = "Point";
54 Shape = "NULL"; 59 Shape = "NULL";
55 CentreCoords = G4ThreeVector(0,0,0); <<
56 Rotx = CLHEP::HepXHat; <<
57 Roty = CLHEP::HepYHat; <<
58 Rotz = CLHEP::HepZHat; <<
59 halfx = 0.; 60 halfx = 0.;
60 halfy = 0.; 61 halfy = 0.;
61 halfz = 0.; 62 halfz = 0.;
62 Radius = 0.; 63 Radius = 0.;
63 Radius0 = 0.; 64 Radius0 = 0.;
64 SR = 0.; 65 SR = 0.;
65 SX = 0.; 66 SX = 0.;
66 SY = 0.; 67 SY = 0.;
67 ParAlpha = 0.; 68 ParAlpha = 0.;
68 ParTheta = 0.; 69 ParTheta = 0.;
69 ParPhi = 0.; 70 ParPhi = 0.;
>> 71 CentreCoords = G4ThreeVector(0., 0., 0.);
>> 72 Rotx = HepXHat;
>> 73 Roty = HepYHat;
>> 74 Rotz = HepZHat;
>> 75 Confine = false; //If true confines source distribution to VolName
70 VolName = "NULL"; 76 VolName = "NULL";
>> 77 SideRefVec1 = HepXHat; // x-axis
>> 78 SideRefVec2 = HepYHat; // y-axis
>> 79 SideRefVec3 = HepZHat; // z-axis
71 verbosityLevel = 0 ; 80 verbosityLevel = 0 ;
72 G4MUTEXINIT(a_mutex); << 81 gNavigator = G4TransportationManager::GetTransportationManager()
>> 82 ->GetNavigatorForTracking();
73 } 83 }
74 84
75 G4SPSPosDistribution::~G4SPSPosDistribution() 85 G4SPSPosDistribution::~G4SPSPosDistribution()
76 { 86 {
77 G4MUTEXDESTROY(a_mutex); <<
78 } 87 }
79 88
80 void G4SPSPosDistribution::SetPosDisType(const << 89 void G4SPSPosDistribution::SetPosDisType(G4String PosType)
81 { 90 {
82 SourcePosType = PosType; 91 SourcePosType = PosType;
83 } 92 }
84 93
85 void G4SPSPosDistribution::SetPosDisShape(cons << 94 void G4SPSPosDistribution::SetPosDisShape(G4String shapeType)
86 { 95 {
87 Shape = shapeType; 96 Shape = shapeType;
88 } 97 }
89 98
90 void G4SPSPosDistribution::SetCentreCoords(con << 99 void G4SPSPosDistribution::SetCentreCoords(G4ThreeVector coordsOfCentre)
91 { 100 {
92 CentreCoords = coordsOfCentre; 101 CentreCoords = coordsOfCentre;
93 } 102 }
94 103
95 void G4SPSPosDistribution::SetPosRot1(const G4 << 104 void G4SPSPosDistribution::SetPosRot1(G4ThreeVector posrot1)
96 { 105 {
97 // This should be x' 106 // This should be x'
98 <<
99 Rotx = posrot1; 107 Rotx = posrot1;
100 if(verbosityLevel == 2) 108 if(verbosityLevel == 2)
101 { << 109 {
102 G4cout << "Vector x' " << Rotx << G4endl; << 110 G4cout << "Vector x' " << Rotx << G4endl;
103 } << 111 }
104 GenerateRotationMatrices(); 112 GenerateRotationMatrices();
105 } 113 }
106 114
107 void G4SPSPosDistribution::SetPosRot2(const G4 << 115 void G4SPSPosDistribution::SetPosRot2(G4ThreeVector posrot2)
108 { 116 {
109 // This is a vector in the plane x'y' but ne << 117 // This is a vector in the plane x'y' but need not
110 << 118 // be y'
111 Roty = posrot2; 119 Roty = posrot2;
112 if(verbosityLevel == 2) 120 if(verbosityLevel == 2)
113 { << 121 {
114 G4cout << "The vector in the x'-y' plane " << 122 G4cout << "The vector in the x'-y' plane " << Roty << G4endl;
115 } << 123 }
116 GenerateRotationMatrices(); 124 GenerateRotationMatrices();
117 } 125 }
118 126
119 void G4SPSPosDistribution::SetHalfX(G4double x 127 void G4SPSPosDistribution::SetHalfX(G4double xhalf)
120 { 128 {
121 halfx = xhalf; 129 halfx = xhalf;
122 } 130 }
123 131
124 void G4SPSPosDistribution::SetHalfY(G4double y 132 void G4SPSPosDistribution::SetHalfY(G4double yhalf)
125 { 133 {
126 halfy = yhalf; 134 halfy = yhalf;
127 } 135 }
128 136
129 void G4SPSPosDistribution::SetHalfZ(G4double z 137 void G4SPSPosDistribution::SetHalfZ(G4double zhalf)
130 { 138 {
131 halfz = zhalf; 139 halfz = zhalf;
132 } 140 }
133 141
134 void G4SPSPosDistribution::SetRadius(G4double << 142 void G4SPSPosDistribution::SetRadius(G4double rad)
135 { 143 {
136 Radius = rds; << 144 Radius = rad;
137 } 145 }
138 146
139 void G4SPSPosDistribution::SetRadius0(G4double << 147 void G4SPSPosDistribution::SetRadius0(G4double rad)
140 { 148 {
141 Radius0 = rds; << 149 Radius0 = rad;
142 } 150 }
143 151
144 void G4SPSPosDistribution::SetBeamSigmaInR(G4d 152 void G4SPSPosDistribution::SetBeamSigmaInR(G4double r)
145 { 153 {
146 SX = SY = r; <<
147 SR = r; 154 SR = r;
>> 155 SX = SY = r/std::sqrt(2.);
148 } 156 }
149 157
150 void G4SPSPosDistribution::SetBeamSigmaInX(G4d 158 void G4SPSPosDistribution::SetBeamSigmaInX(G4double r)
151 { 159 {
152 SX = r; 160 SX = r;
153 } 161 }
154 162
155 void G4SPSPosDistribution::SetBeamSigmaInY(G4d 163 void G4SPSPosDistribution::SetBeamSigmaInY(G4double r)
156 { 164 {
157 SY = r; 165 SY = r;
158 } 166 }
159 167
160 void G4SPSPosDistribution::SetParAlpha(G4doubl 168 void G4SPSPosDistribution::SetParAlpha(G4double paralp)
161 { 169 {
162 ParAlpha = paralp; 170 ParAlpha = paralp;
163 } 171 }
164 172
165 void G4SPSPosDistribution::SetParTheta(G4doubl 173 void G4SPSPosDistribution::SetParTheta(G4double parthe)
166 { 174 {
167 ParTheta = parthe; 175 ParTheta = parthe;
168 } 176 }
169 177
170 void G4SPSPosDistribution::SetParPhi(G4double 178 void G4SPSPosDistribution::SetParPhi(G4double parphi)
171 { 179 {
172 ParPhi = parphi; 180 ParPhi = parphi;
173 } 181 }
174 182
175 const G4String& G4SPSPosDistribution::GetPosDi <<
176 { <<
177 return SourcePosType; <<
178 } <<
179 <<
180 const G4String& G4SPSPosDistribution::GetPosDi <<
181 { <<
182 return Shape; <<
183 } <<
184 <<
185 const G4ThreeVector& G4SPSPosDistribution::Get <<
186 { <<
187 return CentreCoords; <<
188 } <<
189 <<
190 G4double G4SPSPosDistribution::GetHalfX() cons <<
191 { <<
192 return halfx; <<
193 } <<
194 <<
195 G4double G4SPSPosDistribution::GetHalfY() cons <<
196 { <<
197 return halfy; <<
198 } <<
199 <<
200 G4double G4SPSPosDistribution::GetHalfZ() cons <<
201 { <<
202 return halfz; <<
203 } <<
204 <<
205 G4double G4SPSPosDistribution::GetRadius() con <<
206 { <<
207 return Radius; <<
208 } <<
209 <<
210 void G4SPSPosDistribution::SetBiasRndm (G4SPSR <<
211 { <<
212 G4AutoLock l(&a_mutex); <<
213 PosRndm = a; <<
214 } <<
215 <<
216 void G4SPSPosDistribution::SetVerbosity(G4int <<
217 { <<
218 verbosityLevel = a; <<
219 } <<
220 <<
221 const G4String& G4SPSPosDistribution::GetSourc <<
222 { <<
223 return SourcePosType; <<
224 } <<
225 <<
226 const G4ThreeVector& G4SPSPosDistribution::Get <<
227 { <<
228 return ThreadData.Get().CParticlePos; <<
229 } <<
230 <<
231 const G4ThreeVector& G4SPSPosDistribution::Get <<
232 { <<
233 return ThreadData.Get().CSideRefVec1; <<
234 } <<
235 <<
236 const G4ThreeVector& G4SPSPosDistribution::Get <<
237 { <<
238 return ThreadData.Get().CSideRefVec2; <<
239 } <<
240 <<
241 const G4ThreeVector& G4SPSPosDistribution::Get <<
242 { <<
243 return ThreadData.Get().CSideRefVec3; <<
244 } <<
245 <<
246 void G4SPSPosDistribution::GenerateRotationMat 183 void G4SPSPosDistribution::GenerateRotationMatrices()
247 { 184 {
248 // This takes in 2 vectors, x' and one in th 185 // This takes in 2 vectors, x' and one in the plane x'-y',
249 // and from these takes a cross product to c 186 // and from these takes a cross product to calculate z'.
250 // Then a cross product is taken between x' << 187 // Then a cross product is taken between x' and z' to give
251 << 188 // y'.
252 Rotx = Rotx.unit(); // x' 189 Rotx = Rotx.unit(); // x'
253 Roty = Roty.unit(); // vector in x'y' plane 190 Roty = Roty.unit(); // vector in x'y' plane
254 Rotz = Rotx.cross(Roty); // z' 191 Rotz = Rotx.cross(Roty); // z'
255 Rotz = Rotz.unit(); << 192 Roty = Rotz.cross(Rotx); // y'
256 Roty =Rotz.cross(Rotx); // y' <<
257 Roty =Roty.unit(); <<
258 if(verbosityLevel == 2) 193 if(verbosityLevel == 2)
259 { << 194 {
260 G4cout << "The new axes, x', y', z' " << 195 G4cout << "The new axes, x', y', z' " << Rotx << " " << Roty << " " << Rotz << G4endl;
261 << Rotx << " " << Roty << " " << Ro << 196 }
262 } <<
263 } 197 }
264 198
265 void G4SPSPosDistribution::ConfineSourceToVolu << 199 void G4SPSPosDistribution::ConfineSourceToVolume(G4String Vname)
266 { 200 {
267 VolName = Vname; 201 VolName = Vname;
268 if(verbosityLevel == 2) { G4cout << VolName << 202 if(verbosityLevel == 2)
269 << 203 G4cout << VolName << G4endl;
270 if(VolName=="NULL") << 204 G4VPhysicalVolume *tempPV = NULL;
271 { << 205 G4PhysicalVolumeStore *PVStore = 0;
272 if(verbosityLevel >= 1) << 206 G4String theRequiredVolumeName = VolName;
273 { G4cout << "Volume confinement is set off << 207 PVStore = G4PhysicalVolumeStore::GetInstance();
274 Confine = false; << 208 G4int i = 0;
275 return; << 209 G4bool found = false;
>> 210 if(verbosityLevel == 2)
>> 211 G4cout << PVStore->size() << G4endl;
>> 212 while (!found && i<G4int(PVStore->size())) {
>> 213 tempPV = (*PVStore)[i];
>> 214 found = tempPV->GetName() == theRequiredVolumeName;
>> 215 if(verbosityLevel == 2)
>> 216 G4cout << i << " " << " " << tempPV->GetName() << " " << theRequiredVolumeName << " " << found << G4endl;
>> 217 if (!found)
>> 218 {i++;}
276 } 219 }
277 << 220 // found = true then the volume exists else it doesnt.
278 G4VPhysicalVolume* tempPV = nullptr; << 221 if(found == true)
279 G4PhysicalVolumeStore* PVStore = G4PhysicalV <<
280 if(verbosityLevel == 2) { G4cout << PVStore- <<
281 <<
282 tempPV = PVStore->GetVolume(VolName); <<
283 <<
284 // the volume exists else it doesn't <<
285 // <<
286 if (tempPV != nullptr) <<
287 { <<
288 if(verbosityLevel >= 1) <<
289 { 222 {
290 G4cout << "Volume " << VolName << " exis << 223 if(verbosityLevel >= 1)
>> 224 G4cout << "Volume " << VolName << " exists" << G4endl;
>> 225 Confine = true;
291 } 226 }
292 Confine = true; <<
293 } <<
294 else 227 else
295 { << 228 {
296 G4cout << " **** Error: Volume <" << VolNa << 229 G4cout << " **** Error: Volume does not exist **** " << G4endl;
297 << "> does not exist **** " << G4en << 230 G4cout << " Ignoring confine condition" << G4endl;
298 G4cout << " Ignoring confine condition" << << 231 Confine = false;
299 Confine = false; << 232 VolName = "NULL";
300 VolName = "NULL"; << 233 }
301 } << 234
302 } 235 }
303 236
304 void G4SPSPosDistribution::GeneratePointSource << 237 void G4SPSPosDistribution::GeneratePointSource()
305 { 238 {
306 // Generates Points given the point source << 239 // Generates Points given the point source.
307 <<
308 if(SourcePosType == "Point") 240 if(SourcePosType == "Point")
309 { << 241 particle_position = CentreCoords;
310 pos = CentreCoords; <<
311 } <<
312 else 242 else
313 { <<
314 if(verbosityLevel >= 1) 243 if(verbosityLevel >= 1)
315 { << 244 G4cout << "Error SourcePosType is not set to Point" << G4endl;
316 G4cerr << "Error SourcePosType is not se <<
317 } <<
318 } <<
319 } 245 }
320 246
321 void G4SPSPosDistribution::GeneratePointsInBea << 247 void G4SPSPosDistribution::GeneratePointsInBeam()
322 { 248 {
323 G4double x, y, z; 249 G4double x, y, z;
324 250
325 G4ThreeVector RandPos; 251 G4ThreeVector RandPos;
326 G4double tempx, tempy, tempz; 252 G4double tempx, tempy, tempz;
327 z = 0.; 253 z = 0.;
328 254
329 // Private Method to create points in a plan 255 // Private Method to create points in a plane
330 // <<
331 if(Shape == "Circle") 256 if(Shape == "Circle")
332 { << 257 {
333 x = Radius + 100.; << 258 x = Radius + 100.;
334 y = Radius + 100.; << 259 y = Radius + 100.;
335 while(std::sqrt((x*x) + (y*y)) > Radius) << 260 while(std::sqrt((x*x) + (y*y)) > Radius)
336 { << 261 {
337 x = PosRndm->GenRandX(); << 262 x = posRndm->GenRandX();
338 y = PosRndm->GenRandY(); << 263 y = posRndm->GenRandY();
339 << 264
340 x = (x*2.*Radius) - Radius; << 265 x = (x*2.*Radius) - Radius;
341 y = (y*2.*Radius) - Radius; << 266 y = (y*2.*Radius) - Radius;
342 } << 267 }
343 x += G4RandGauss::shoot(0.0,SX) ; << 268 x += G4RandGauss::shoot(0.0,SX) ;
344 y += G4RandGauss::shoot(0.0,SY) ; << 269 y += G4RandGauss::shoot(0.0,SY) ;
345 } << 270 }
346 else 271 else
347 { << 272 {
348 // All other cases default to Rectangle ca << 273 // all other cases default to Rectangle case
349 // << 274 x = posRndm->GenRandX();
350 x = PosRndm->GenRandX(); << 275 y = posRndm->GenRandY();
351 y = PosRndm->GenRandY(); << 276 x = (x*2.*halfx) - halfx;
352 x = (x*2.*halfx) - halfx; << 277 y = (y*2.*halfy) - halfy;
353 y = (y*2.*halfy) - halfy; << 278 x += G4RandGauss::shoot(0.0,SX);
354 x += G4RandGauss::shoot(0.0,SX); << 279 y += G4RandGauss::shoot(0.0,SY);
355 y += G4RandGauss::shoot(0.0,SY); << 280 }
356 } <<
357 <<
358 // Apply Rotation Matrix 281 // Apply Rotation Matrix
359 // x * Rotx, y * Roty and z * Rotz 282 // x * Rotx, y * Roty and z * Rotz
360 // <<
361 if(verbosityLevel >= 2) 283 if(verbosityLevel >= 2)
362 { << 284 {
363 G4cout << "Raw position " << x << "," << y << 285 G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
364 } << 286 }
365 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z 287 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
366 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z 288 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
367 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z 289 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
368 290
369 RandPos.setX(tempx); 291 RandPos.setX(tempx);
370 RandPos.setY(tempy); 292 RandPos.setY(tempy);
371 RandPos.setZ(tempz); 293 RandPos.setZ(tempz);
372 294
373 // Translate 295 // Translate
374 // << 296 particle_position = CentreCoords + RandPos;
375 pos = CentreCoords + RandPos; <<
376 if(verbosityLevel >= 1) 297 if(verbosityLevel >= 1)
377 { <<
378 if(verbosityLevel >= 2) <<
379 { 298 {
380 G4cout << "Rotated Position " << RandPos << 299 if(verbosityLevel >= 2)
>> 300 {
>> 301 G4cout << "Rotated Position " << RandPos << G4endl;
>> 302 }
>> 303 G4cout << "Rotated and Translated position " << particle_position << G4endl;
381 } 304 }
382 G4cout << "Rotated and Translated position <<
383 } <<
384 } 305 }
385 306
386 void G4SPSPosDistribution::GeneratePointsInPla << 307 void G4SPSPosDistribution::GeneratePointsInPlane()
387 { 308 {
388 G4double x, y, z; 309 G4double x, y, z;
389 G4double expression; 310 G4double expression;
390 G4ThreeVector RandPos; 311 G4ThreeVector RandPos;
391 G4double tempx, tempy, tempz; 312 G4double tempx, tempy, tempz;
392 x = y = z = 0.; 313 x = y = z = 0.;
393 thread_data_t& td = ThreadData.Get(); <<
394 314
395 if(SourcePosType != "Plane" && verbosityLeve 315 if(SourcePosType != "Plane" && verbosityLevel >= 1)
396 { << 316 G4cout << "Error: SourcePosType is not Plane" << G4endl;
397 G4cerr << "Error: SourcePosType is not Pla <<
398 } <<
399 317
400 // Private Method to create points in a plan 318 // Private Method to create points in a plane
401 // <<
402 if(Shape == "Circle") 319 if(Shape == "Circle")
403 { <<
404 x = Radius + 100.; <<
405 y = Radius + 100.; <<
406 while(std::sqrt((x*x) + (y*y)) > Radius) <<
407 { 320 {
408 x = PosRndm->GenRandX(); << 321 x = Radius + 100.;
409 y = PosRndm->GenRandY(); << 322 y = Radius + 100.;
410 << 323 while(std::sqrt((x*x) + (y*y)) > Radius)
411 x = (x*2.*Radius) - Radius; << 324 {
412 y = (y*2.*Radius) - Radius; << 325 x = posRndm->GenRandX();
>> 326 y = posRndm->GenRandY();
>> 327
>> 328 x = (x*2.*Radius) - Radius;
>> 329 y = (y*2.*Radius) - Radius;
>> 330 }
413 } 331 }
414 } <<
415 else if(Shape == "Annulus") 332 else if(Shape == "Annulus")
416 { <<
417 x = Radius + 100.; <<
418 y = Radius + 100.; <<
419 while(std::sqrt((x*x) + (y*y)) > Radius <<
420 || std::sqrt((x*x) + (y*y)) < Radius0 ) <<
421 { 333 {
422 x = PosRndm->GenRandX(); << 334 x = Radius + 100.;
423 y = PosRndm->GenRandY(); << 335 y = Radius + 100.;
424 << 336 while(std::sqrt((x*x) + (y*y)) > Radius || std::sqrt((x*x) + (y*y)) < Radius0 )
425 x = (x*2.*Radius) - Radius; << 337 {
426 y = (y*2.*Radius) - Radius; << 338 x = posRndm->GenRandX();
>> 339 y = posRndm->GenRandY();
>> 340
>> 341 x = (x*2.*Radius) - Radius;
>> 342 y = (y*2.*Radius) - Radius;
>> 343 }
427 } 344 }
428 } <<
429 else if(Shape == "Ellipse") 345 else if(Shape == "Ellipse")
430 { <<
431 expression = 20.; <<
432 while(expression > 1.) <<
433 { 346 {
434 x = PosRndm->GenRandX(); << 347 expression = 20.;
435 y = PosRndm->GenRandY(); << 348 while(expression > 1.)
>> 349 {
>> 350 x = posRndm->GenRandX();
>> 351 y = posRndm->GenRandY();
436 352
437 x = (x*2.*halfx) - halfx; << 353 x = (x*2.*halfx) - halfx;
438 y = (y*2.*halfy) - halfy; << 354 y = (y*2.*halfy) - halfy;
439 355
440 expression = ((x*x)/(halfx*halfx)) + ((y << 356 expression = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy));
>> 357 }
441 } 358 }
442 } <<
443 else if(Shape == "Square") 359 else if(Shape == "Square")
444 { << 360 {
445 x = PosRndm->GenRandX(); << 361 x = posRndm->GenRandX();
446 y = PosRndm->GenRandY(); << 362 y = posRndm->GenRandY();
447 x = (x*2.*halfx) - halfx; << 363 x = (x*2.*halfx) - halfx;
448 y = (y*2.*halfy) - halfy; << 364 y = (y*2.*halfy) - halfy;
449 } << 365 }
450 else if(Shape == "Rectangle") 366 else if(Shape == "Rectangle")
451 { << 367 {
452 x = PosRndm->GenRandX(); << 368 x = posRndm->GenRandX();
453 y = PosRndm->GenRandY(); << 369 y = posRndm->GenRandY();
454 x = (x*2.*halfx) - halfx; << 370 x = (x*2.*halfx) - halfx;
455 y = (y*2.*halfy) - halfy; << 371 y = (y*2.*halfy) - halfy;
456 } << 372 }
457 else 373 else
458 { <<
459 G4cout << "Shape not one of the plane type 374 G4cout << "Shape not one of the plane types" << G4endl;
460 } <<
461 375
462 // Apply Rotation Matrix 376 // Apply Rotation Matrix
463 // x * Rotx, y * Roty and z * Rotz 377 // x * Rotx, y * Roty and z * Rotz
464 // <<
465 if(verbosityLevel == 2) 378 if(verbosityLevel == 2)
466 { << 379 {
467 G4cout << "Raw position " << x << "," << y << 380 G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
468 } << 381 }
469 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z 382 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
470 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z 383 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
471 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z 384 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
472 385
473 RandPos.setX(tempx); 386 RandPos.setX(tempx);
474 RandPos.setY(tempy); 387 RandPos.setY(tempy);
475 RandPos.setZ(tempz); 388 RandPos.setZ(tempz);
476 389
477 // Translate 390 // Translate
478 // << 391 particle_position = CentreCoords + RandPos;
479 pos = CentreCoords + RandPos; <<
480 if(verbosityLevel >= 1) 392 if(verbosityLevel >= 1)
481 { <<
482 if(verbosityLevel == 2) <<
483 { 393 {
484 G4cout << "Rotated Position " << RandPos << 394 if(verbosityLevel == 2)
>> 395 {
>> 396 G4cout << "Rotated Position " << RandPos << G4endl;
>> 397 }
>> 398 G4cout << "Rotated and Translated position " << particle_position << G4endl;
485 } 399 }
486 G4cout << "Rotated and Translated position <<
487 } <<
488 400
489 // For Cosine-Law make SideRefVecs = to Rota 401 // For Cosine-Law make SideRefVecs = to Rotation matrix vectors
490 // Note: these need to be thread-local, use << 402 SideRefVec1 = Rotx;
491 // << 403 SideRefVec2 = Roty;
492 td.CSideRefVec1 = Rotx; << 404 SideRefVec3 = Rotz;
493 td.CSideRefVec2 = Roty; <<
494 td.CSideRefVec3 = Rotz; <<
495 <<
496 // If rotation matrix z' point to origin the 405 // If rotation matrix z' point to origin then invert the matrix
497 // So that SideRefVecs point away << 406 // So that SideRefVecs point away.
498 // <<
499 if((CentreCoords.x() > 0. && Rotz.x() < 0.) 407 if((CentreCoords.x() > 0. && Rotz.x() < 0.)
500 || (CentreCoords.x() < 0. && Rotz.x() > 0 408 || (CentreCoords.x() < 0. && Rotz.x() > 0.)
501 || (CentreCoords.y() > 0. && Rotz.y() < 0 409 || (CentreCoords.y() > 0. && Rotz.y() < 0.)
502 || (CentreCoords.y() < 0. && Rotz.y() > 0 410 || (CentreCoords.y() < 0. && Rotz.y() > 0.)
503 || (CentreCoords.z() > 0. && Rotz.z() < 0 411 || (CentreCoords.z() > 0. && Rotz.z() < 0.)
504 || (CentreCoords.z() < 0. && Rotz.z() > 0 412 || (CentreCoords.z() < 0. && Rotz.z() > 0.))
505 { << 413 {
506 // Invert y and z << 414 // Invert y and z.
507 // << 415 SideRefVec2 = -SideRefVec2;
508 td.CSideRefVec2 = - td.CSideRefVec2; << 416 SideRefVec3 = -SideRefVec3;
509 td.CSideRefVec3 = - td.CSideRefVec3; << 417 }
510 } <<
511 if(verbosityLevel == 2) 418 if(verbosityLevel == 2)
512 { << 419 {
513 G4cout << "Reference vectors for cosine-la << 420 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
514 << td.CSideRefVec1<< " " << td.CSid << 421 }
515 << " " << td.CSideRefVec3 << G4endl <<
516 } <<
517 } 422 }
518 423
519 void G4SPSPosDistribution::GeneratePointsOnSur << 424 void G4SPSPosDistribution::GeneratePointsOnSurface()
520 { 425 {
521 // Private method to create points on a surf << 426 //Private method to create points on a surface
522 // <<
523 G4double theta, phi; 427 G4double theta, phi;
524 G4double x, y, z; 428 G4double x, y, z;
525 x = y = z = 0.; 429 x = y = z = 0.;
526 G4double expression; <<
527 G4ThreeVector RandPos; 430 G4ThreeVector RandPos;
>> 431 // G4double tempx, tempy, tempz;
528 432
529 thread_data_t& td = ThreadData.Get(); <<
530 if(SourcePosType != "Surface" && verbosityLe 433 if(SourcePosType != "Surface" && verbosityLevel >= 1)
531 { <<
532 G4cout << "Error SourcePosType not Surface 434 G4cout << "Error SourcePosType not Surface" << G4endl;
533 } <<
534 435
535 if(Shape == "Sphere") 436 if(Shape == "Sphere")
536 { << 437 {
537 theta = PosRndm->GenRandPosTheta(); << 438 G4double tantheta;
538 phi = PosRndm->GenRandPosPhi(); << 439 theta = posRndm->GenRandPosTheta();
539 theta = std::acos(1. - 2.*theta); // theta << 440 phi = posRndm->GenRandPosPhi();
540 phi = phi * 2. * pi; << 441 theta = std::acos(1. - 2.*theta); // theta isotropic
>> 442 phi = phi * 2. * pi;
>> 443 tantheta = std::tan(theta);
541 444
542 x = Radius * std::sin(theta) * std::cos(ph << 445 x = Radius * std::sin(theta) * std::cos(phi);
543 y = Radius * std::sin(theta) * std::sin(ph << 446 y = Radius * std::sin(theta) * std::sin(phi);
544 z = Radius * std::cos(theta); << 447 z = Radius * std::cos(theta);
545 448
546 RandPos.setX(x); << 449 RandPos.setX(x);
547 RandPos.setY(y); << 450 RandPos.setY(y);
548 RandPos.setZ(z); << 451 RandPos.setZ(z);
549 << 452
550 // Cosine-law (not a good idea to use this << 453 // Cosine-law (not a good idea to use this here)
551 // << 454 G4ThreeVector zdash(x,y,z);
552 G4ThreeVector zdash(x,y,z); << 455 zdash = zdash.unit();
553 zdash = zdash.unit(); << 456 G4ThreeVector xdash = Rotz.cross(zdash);
554 G4ThreeVector xdash = Rotz.cross(zdash); << 457 G4ThreeVector ydash = xdash.cross(zdash);
555 G4ThreeVector ydash = xdash.cross(zdash); << 458 SideRefVec1 = xdash.unit();
556 td.CSideRefVec1 = xdash.unit(); << 459 SideRefVec2 = ydash.unit();
557 td.CSideRefVec2 = ydash.unit(); << 460 SideRefVec3 = zdash.unit();
558 td.CSideRefVec3 = zdash.unit(); << 461 }
559 } <<
560 else if(Shape == "Ellipsoid") 462 else if(Shape == "Ellipsoid")
561 { << 463 {
562 G4double minphi, maxphi, middlephi; << 464 G4double theta, phi, minphi, maxphi, middlephi;
563 G4double answer, constant; << 465 G4double answer, constant;
564 466
565 constant = pi/(halfx*halfx) + pi/(halfy*ha << 467 constant = pi/(halfx*halfx) + pi/(halfy*halfy) +
>> 468 twopi/(halfz*halfz);
566 469
567 // Simplified approach << 470 // simplified approach
568 // << 471 theta = posRndm->GenRandPosTheta();
569 theta = PosRndm->GenRandPosTheta(); << 472 phi = posRndm->GenRandPosPhi();
570 phi = PosRndm->GenRandPosPhi(); <<
571 473
572 theta = std::acos(1. - 2.*theta); << 474 theta = std::acos(1. - 2.*theta);
573 minphi = 0.; << 475 minphi = 0.;
574 maxphi = twopi; << 476 maxphi = twopi;
575 while(maxphi-minphi > 0.) << 477 while(maxphi-minphi > 0.)
576 { << 478 {
577 middlephi = (maxphi+minphi)/2.; << 479 middlephi = (maxphi+minphi)/2.;
578 answer = (1./(halfx*halfx))*(middlephi/2 << 480 answer = (1./(halfx*halfx))*(middlephi/2. + std::sin(2*middlephi)/4.)
579 + (1./(halfy*halfy))*(middlephi/2 << 481 + (1./(halfy*halfy))*(middlephi/2. - std::sin(2*middlephi)/4.)
580 + middlephi/(halfz*halfz); << 482 + middlephi/(halfz*halfz);
581 answer = answer/constant; << 483 answer = answer/constant;
582 if(answer > phi) maxphi = middlephi; << 484 if(answer > phi) maxphi = middlephi;
583 if(answer < phi) minphi = middlephi; << 485 if(answer < phi) minphi = middlephi;
584 if(std::fabs(answer-phi) <= 0.00001) << 486 if(std::fabs(answer-phi) <= 0.00001)
585 { << 487 {
586 minphi = maxphi +1; << 488 minphi = maxphi +1;
587 phi = middlephi; << 489 phi = middlephi;
588 } << 490 }
589 } << 491 }
590 << 492
591 // x,y and z form a unit vector. Put this << 493 x = std::sin(theta)*std::cos(phi);
592 // << 494 y = std::sin(theta)*std::sin(phi);
593 x = std::sin(theta)*std::cos(phi); << 495 z = std::cos(theta);
594 y = std::sin(theta)*std::sin(phi); << 496 // x,y and z form a unit vector. Put this onto the ellipse.
595 z = std::cos(theta); << 497 G4double lhs;
596 << 498 // solve for x
597 G4double lhs; << 499 G4double numYinX = y/x;
598 << 500 G4double numZinX = z/x;
599 // Solve for x << 501 G4double tempxvar;
600 // << 502 tempxvar= 1./(halfx*halfx)+(numYinX*numYinX)/(halfy*halfy)
601 G4double numYinX = y/x; << 503 + (numZinX*numZinX)/(halfz*halfz);
602 G4double numZinX = z/x; << 504
603 G4double tempxvar; << 505 tempxvar = 1./tempxvar;
604 tempxvar = 1./(halfx*halfx)+(numYinX*numYi << 506 G4double coordx = std::sqrt(tempxvar);
605 + (numZinX*numZinX)/(halfz*halfz) <<
606 tempxvar = 1./tempxvar; <<
607 G4double coordx = std::sqrt(tempxvar); <<
608 507
609 // Solve for y << 508 //solve for y
610 // << 509 G4double numXinY = x/y;
611 G4double numXinY = x/y; << 510 G4double numZinY = z/y;
612 G4double numZinY = z/y; << 511 G4double tempyvar;
613 G4double tempyvar; << 512 tempyvar=(numXinY*numXinY)/(halfx*halfx)+1./(halfy*halfy)
614 tempyvar = (numXinY*numXinY)/(halfx*halfx) << 513 +(numZinY*numZinY)/(halfz*halfz);
615 + (numZinY*numZinY)/(halfz*halfz) << 514 tempyvar = 1./tempyvar;
616 tempyvar = 1./tempyvar; << 515 G4double coordy = std::sqrt(tempyvar);
617 G4double coordy = std::sqrt(tempyvar); <<
618 <<
619 // Solve for z <<
620 // <<
621 G4double numXinZ = x/z; <<
622 G4double numYinZ = y/z; <<
623 G4double tempzvar; <<
624 tempzvar = (numXinZ*numXinZ)/(halfx*halfx) <<
625 + (numYinZ*numYinZ)/(halfy*halfy) <<
626 tempzvar = 1./tempzvar; <<
627 G4double coordz = std::sqrt(tempzvar); <<
628 <<
629 lhs = std::sqrt((coordx*coordx)/(halfx*hal <<
630 (coordy*coordy)/(halfy*hal <<
631 (coordz*coordz)/(halfz*hal <<
632 516
633 if(std::fabs(lhs-1.) > 0.001 && verbosityL << 517 //solve for z
634 { << 518 G4double numXinZ = x/z;
635 G4cout << "Error: theta, phi not really << 519 G4double numYinZ = y/z;
636 } << 520 G4double tempzvar;
637 << 521 tempzvar=(numXinZ*numXinZ)/(halfx*halfx)
638 // coordx, coordy and coordz are all posit << 522 +(numYinZ*numYinZ)/(halfy*halfy)+1./(halfz*halfz);
639 // << 523 tempzvar = 1./tempzvar;
640 G4double TestRandVar = G4UniformRand(); << 524 G4double coordz = std::sqrt(tempzvar);
641 if(TestRandVar > 0.5) << 525
642 { << 526 lhs = std::sqrt((coordx*coordx)/(halfx*halfx) +
643 coordx = -coordx; << 527 (coordy*coordy)/(halfy*halfy) +
644 } << 528 (coordz*coordz)/(halfz*halfz));
645 TestRandVar = G4UniformRand(); <<
646 if(TestRandVar > 0.5) <<
647 { <<
648 coordy = -coordy; <<
649 } <<
650 TestRandVar = G4UniformRand(); <<
651 if(TestRandVar > 0.5) <<
652 { <<
653 coordz = -coordz; <<
654 } <<
655 <<
656 RandPos.setX(coordx); <<
657 RandPos.setY(coordy); <<
658 RandPos.setZ(coordz); <<
659 <<
660 // Cosine-law (not a good idea to use this <<
661 // <<
662 G4ThreeVector zdash(coordx,coordy,coordz); <<
663 zdash = zdash.unit(); <<
664 G4ThreeVector xdash = Rotz.cross(zdash); <<
665 G4ThreeVector ydash = xdash.cross(zdash); <<
666 td.CSideRefVec1 = xdash.unit(); <<
667 td.CSideRefVec2 = ydash.unit(); <<
668 td.CSideRefVec3 = zdash.unit(); <<
669 } <<
670 else if(Shape == "Cylinder") <<
671 { <<
672 G4double AreaTop, AreaBot, AreaLat; <<
673 G4double AreaTotal, prob1, prob2, prob3; <<
674 G4double testrand; <<
675 <<
676 // User giver Radius and z-half length <<
677 // Calculate surface areas, maybe move thi <<
678 // a different method <<
679 <<
680 AreaTop = pi * Radius * Radius; <<
681 AreaBot = AreaTop; <<
682 AreaLat = 2. * pi * Radius * 2. * halfz; <<
683 AreaTotal = AreaTop + AreaBot + AreaLat; <<
684 529
685 prob1 = AreaTop / AreaTotal; << 530 if(std::fabs(lhs-1.) > 0.001 && verbosityLevel >= 1)
686 prob2 = AreaBot / AreaTotal; << 531 G4cout << "Error: theta, phi not really on ellipsoid" << G4endl;
687 prob3 = 1.00 - prob1 - prob2; <<
688 if(std::fabs(prob3 - (AreaLat/AreaTotal)) <<
689 { <<
690 if(verbosityLevel >= 1) <<
691 { <<
692 G4cout << AreaLat/AreaTotal << " " << <<
693 } <<
694 G4cout << "Error in prob3" << G4endl; <<
695 } <<
696 <<
697 // Decide surface to calculate point on. <<
698 532
699 testrand = G4UniformRand(); << 533 // coordx, coordy and coordz are all positive
700 if(testrand <= prob1) // Point on Top sur << 534 G4double TestRandVar = G4UniformRand();
701 { << 535 if(TestRandVar > 0.5)
702 z = halfz; << 536 {
703 x = Radius + 100.; << 537 coordx = -coordx;
704 y = Radius + 100.; << 538 }
705 while(((x*x)+(y*y)) > (Radius*Radius)) << 539 TestRandVar = G4UniformRand();
706 { << 540 if(TestRandVar > 0.5)
707 x = PosRndm->GenRandX(); << 541 {
708 y = PosRndm->GenRandY(); << 542 coordy = -coordy;
709 << 543 }
710 x = x * 2. * Radius; << 544 TestRandVar = G4UniformRand();
711 y = y * 2. * Radius; << 545 if(TestRandVar > 0.5)
712 x = x - Radius; << 546 {
713 y = y - Radius; << 547 coordz = -coordz;
714 } << 548 }
715 // Cosine law << 549
716 // << 550 RandPos.setX(coordx);
717 td.CSideRefVec1 = Rotx; << 551 RandPos.setY(coordy);
718 td.CSideRefVec2 = Roty; << 552 RandPos.setZ(coordz);
719 td.CSideRefVec3 = Rotz; <<
720 } <<
721 else if((testrand > prob1) && (testrand <= <<
722 { // Point on Bot <<
723 z = -halfz; <<
724 x = Radius + 100.; <<
725 y = Radius + 100.; <<
726 while(((x*x)+(y*y)) > (Radius*Radius)) <<
727 { <<
728 x = PosRndm->GenRandX(); <<
729 y = PosRndm->GenRandY(); <<
730 <<
731 x = x * 2. * Radius; <<
732 y = y * 2. * Radius; <<
733 x = x - Radius; <<
734 y = y - Radius; <<
735 } <<
736 // Cosine law <<
737 // <<
738 td.CSideRefVec1 = Rotx; <<
739 td.CSideRefVec2 = -Roty; <<
740 td.CSideRefVec3 = -Rotz; <<
741 } <<
742 else if(testrand > (prob1+prob2)) // Poin <<
743 { <<
744 G4double rand; <<
745 <<
746 rand = PosRndm->GenRandPosPhi(); <<
747 rand = rand * 2. * pi; <<
748 553
749 x = Radius * std::cos(rand); << 554 // Cosine-law (not a good idea to use this here)
750 y = Radius * std::sin(rand); << 555 G4ThreeVector zdash(coordx,coordy,coordz);
751 <<
752 z = PosRndm->GenRandZ(); <<
753 <<
754 z = z * 2. *halfz; <<
755 z = z - halfz; <<
756 <<
757 // Cosine law <<
758 // <<
759 G4ThreeVector zdash(x,y,0.); <<
760 zdash = zdash.unit(); 556 zdash = zdash.unit();
761 G4ThreeVector xdash = Rotz.cross(zdash); 557 G4ThreeVector xdash = Rotz.cross(zdash);
762 G4ThreeVector ydash = xdash.cross(zdash) 558 G4ThreeVector ydash = xdash.cross(zdash);
763 td.CSideRefVec1 = xdash.unit(); << 559 SideRefVec1 = xdash.unit();
764 td.CSideRefVec2 = ydash.unit(); << 560 SideRefVec2 = ydash.unit();
765 td.CSideRefVec3 = zdash.unit(); << 561 SideRefVec3 = zdash.unit();
766 } <<
767 else <<
768 { <<
769 G4cout << "Error: testrand " << testrand <<
770 } <<
771 <<
772 RandPos.setX(x); <<
773 RandPos.setY(y); <<
774 RandPos.setZ(z); <<
775 } <<
776 else if(Shape == "EllipticCylinder") <<
777 { <<
778 G4double AreaTop, AreaBot, AreaLat, AreaTo <<
779 G4double h; <<
780 G4double prob1, prob2, prob3; <<
781 G4double testrand; <<
782 <<
783 // User giver x, y and z-half length <<
784 // Calculate surface areas, maybe move thi <<
785 // a different method <<
786 <<
787 AreaTop = pi * halfx * halfy; <<
788 AreaBot = AreaTop; <<
789 <<
790 // Using circumference approximation by Ra <<
791 // AreaLat = 2*halfz * pi*( 3*(halfx + hal <<
792 // - std::sqrt((3*halfx + halfy) * <<
793 // Using circumference approximation by Ra <<
794 // <<
795 h = ((halfx - halfy)*(halfx - halfy)) / (( <<
796 AreaLat = 2*halfz * pi*(halfx + halfy) <<
797 * (1 + (3*h)/(10 + std::sqrt(4 - 3 <<
798 AreaTotal = AreaTop + AreaBot + AreaLat; <<
799 <<
800 prob1 = AreaTop / AreaTotal; <<
801 prob2 = AreaBot / AreaTotal; <<
802 prob3 = 1.00 - prob1 - prob2; <<
803 if(std::fabs(prob3 - (AreaLat/AreaTotal)) <<
804 { <<
805 if(verbosityLevel >= 1) <<
806 { <<
807 G4cout << AreaLat/AreaTotal << " " << <<
808 } <<
809 G4cout << "Error in prob3" << G4endl; <<
810 } <<
811 <<
812 // Decide surface to calculate point on <<
813 <<
814 testrand = G4UniformRand(); <<
815 if(testrand <= prob1) // Point on Top sur <<
816 { <<
817 z = halfz; <<
818 expression = 20.; <<
819 while(expression > 1.) <<
820 { <<
821 x = PosRndm->GenRandX(); <<
822 y = PosRndm->GenRandY(); <<
823 <<
824 x = (x * 2. * halfx) - halfx; <<
825 y = (y * 2. * halfy) - halfy; <<
826 <<
827 expression = ((x*x)/(halfx*halfx)) + ( <<
828 } <<
829 } <<
830 else if((testrand > prob1) && (testrand <= <<
831 { // Point on Bot <<
832 z = -halfz; <<
833 expression = 20.; <<
834 while(expression > 1.) <<
835 { <<
836 x = PosRndm->GenRandX(); <<
837 y = PosRndm->GenRandY(); <<
838 <<
839 x = (x * 2. * halfx) - halfx; <<
840 y = (y * 2. * halfy) - halfy; <<
841 <<
842 expression = ((x*x)/(halfx*halfx)) + ( <<
843 } <<
844 } 562 }
845 else if(testrand > (prob1+prob2)) // Poin << 563 else if(Shape == "Cylinder")
846 { 564 {
847 G4double rand; << 565 G4double AreaTop, AreaBot, AreaLat;
848 << 566 G4double AreaTotal, prob1, prob2, prob3;
849 rand = PosRndm->GenRandPosPhi(); << 567 G4double testrand;
850 rand = rand * 2. * pi; << 568
851 << 569 // User giver Radius and z-half length
852 x = halfx * std::cos(rand); << 570 // Calculate surface areas, maybe move this to
853 y = halfy * std::sin(rand); << 571 // a different routine.
854 << 572
855 z = PosRndm->GenRandZ(); << 573 AreaTop = pi * Radius * Radius;
>> 574 AreaBot = AreaTop;
>> 575 AreaLat = 2. * pi * Radius * 2. * halfz;
>> 576 AreaTotal = AreaTop + AreaBot + AreaLat;
>> 577
>> 578 prob1 = AreaTop / AreaTotal;
>> 579 prob2 = AreaBot / AreaTotal;
>> 580 prob3 = 1.00 - prob1 - prob2;
>> 581 if(std::fabs(prob3 - (AreaLat/AreaTotal)) >= 0.001)
>> 582 {
>> 583 if(verbosityLevel >= 1)
>> 584 G4cout << AreaLat/AreaTotal << " " << prob3<<G4endl;
>> 585 G4cout << "Error in prob3" << G4endl;
>> 586 }
>> 587
>> 588 // Decide surface to calculate point on.
>> 589
>> 590 testrand = G4UniformRand();
>> 591 if(testrand <= prob1)
>> 592 {
>> 593 //Point on Top surface
>> 594 z = halfz;
>> 595 x = Radius + 100.;
>> 596 y = Radius + 100.;
>> 597 while(((x*x)+(y*y)) > (Radius*Radius))
>> 598 {
>> 599 x = posRndm->GenRandX();
>> 600 y = posRndm->GenRandY();
>> 601
>> 602 x = x * 2. * Radius;
>> 603 y = y * 2. * Radius;
>> 604 x = x - Radius;
>> 605 y = y - Radius;
>> 606 }
>> 607 // Cosine law
>> 608 SideRefVec1 = Rotx;
>> 609 SideRefVec2 = Roty;
>> 610 SideRefVec3 = Rotz;
>> 611 }
>> 612 else if((testrand > prob1) && (testrand <= (prob1 + prob2)))
>> 613 {
>> 614 //Point on Bottom surface
>> 615 z = -halfz;
>> 616 x = Radius + 100.;
>> 617 y = Radius + 100.;
>> 618 while(((x*x)+(y*y)) > (Radius*Radius))
>> 619 {
>> 620 x = posRndm->GenRandX();
>> 621 y = posRndm->GenRandY();
>> 622
>> 623 x = x * 2. * Radius;
>> 624 y = y * 2. * Radius;
>> 625 x = x - Radius;
>> 626 y = y - Radius;
>> 627 }
>> 628 // Cosine law
>> 629 SideRefVec1 = Rotx;
>> 630 SideRefVec2 = -Roty;
>> 631 SideRefVec3 = -Rotz;
>> 632 }
>> 633 else if(testrand > (prob1+prob2))
>> 634 {
>> 635 G4double rand;
>> 636 //Point on Lateral Surface
>> 637
>> 638 rand = posRndm->GenRandPosPhi();
>> 639 rand = rand * 2. * pi;
>> 640
>> 641 x = Radius * std::cos(rand);
>> 642 y = Radius * std::sin(rand);
>> 643
>> 644 z = posRndm->GenRandZ();
>> 645
>> 646 z = z * 2. * halfz;
>> 647 z = z - halfz;
>> 648
>> 649 // Cosine law
>> 650 G4ThreeVector zdash(x,y,0.);
>> 651 zdash = zdash.unit();
>> 652 G4ThreeVector xdash = Rotz.cross(zdash);
>> 653 G4ThreeVector ydash = xdash.cross(zdash);
>> 654 SideRefVec1 = xdash.unit();
>> 655 SideRefVec2 = ydash.unit();
>> 656 SideRefVec3 = zdash.unit();
>> 657 }
>> 658 else
>> 659 G4cout << "Error: testrand " << testrand << G4endl;
>> 660
>> 661 RandPos.setX(x);
>> 662 RandPos.setY(y);
>> 663 RandPos.setZ(z);
856 664
857 z = (z * 2. * halfz) - halfz; <<
858 } <<
859 else <<
860 { <<
861 G4cout << "Error: testrand " << testrand <<
862 } 665 }
863 <<
864 RandPos.setX(x); <<
865 RandPos.setY(y); <<
866 RandPos.setZ(z); <<
867 <<
868 // Cosine law <<
869 // <<
870 G4ThreeVector zdash(x,y,z); <<
871 zdash = zdash.unit(); <<
872 G4ThreeVector xdash = Rotz.cross(zdash); <<
873 G4ThreeVector ydash = xdash.cross(zdash); <<
874 td.CSideRefVec1 = xdash.unit(); <<
875 td.CSideRefVec2 = ydash.unit(); <<
876 td.CSideRefVec3 = zdash.unit(); <<
877 } <<
878 else if(Shape == "Para") 666 else if(Shape == "Para")
879 { << 667 {
880 // Right Parallelepiped. << 668 G4double testrand;
881 // User gives x,y,z half lengths and ParAl << 669 //Right Parallelepiped.
882 // ParTheta and ParPhi << 670 // User gives x,y,z half lengths and ParAlpha
883 // +x = <1, -x >1 & <2, +y >2 & <3, -y >3 << 671 // ParTheta and ParPhi
884 // +z >4 & < 5, -z >5 &<6 << 672 // +x = <1, -x >1 & <2, +y >2 & <3, -y >3 &<4
885 << 673 // +z >4 & < 5, -z >5 &<6.
886 G4double testrand = G4UniformRand(); << 674 testrand = G4UniformRand();
887 G4double AreaX = halfy * halfz * 4.; << 675 G4double AreaX = halfy * halfz * 4.;
888 G4double AreaY = halfx * halfz * 4.; << 676 G4double AreaY = halfx * halfz * 4.;
889 G4double AreaZ = halfx * halfy * 4.; << 677 G4double AreaZ = halfx * halfy * 4.;
890 G4double AreaTotal = 2*(AreaX + AreaY + Ar << 678 G4double AreaTotal = 2*(AreaX + AreaY + AreaZ);
891 G4double Probs[6]; << 679 G4double Probs[6];
892 Probs[0] = AreaX/AreaTotal; << 680 Probs[0] = AreaX/AreaTotal;
893 Probs[1] = Probs[0] + AreaX/AreaTotal; << 681 Probs[1] = Probs[0] + AreaX/AreaTotal;
894 Probs[2] = Probs[1] + AreaY/AreaTotal; << 682 Probs[2] = Probs[1] + AreaY/AreaTotal;
895 Probs[3] = Probs[2] + AreaY/AreaTotal; << 683 Probs[3] = Probs[2] + AreaY/AreaTotal;
896 Probs[4] = Probs[3] + AreaZ/AreaTotal; << 684 Probs[4] = Probs[3] + AreaZ/AreaTotal;
897 Probs[5] = Probs[4] + AreaZ/AreaTotal; << 685 Probs[5] = Probs[4] + AreaZ/AreaTotal;
898 686
899 x = PosRndm->GenRandX(); << 687 x = posRndm->GenRandX();
900 y = PosRndm->GenRandY(); << 688 y = posRndm->GenRandY();
901 z = PosRndm->GenRandZ(); << 689 z = posRndm->GenRandZ();
902 690
903 x = x * halfx * 2.; << 691 x = x * halfx * 2.;
904 x = x - halfx; << 692 x = x - halfx;
905 y = y * halfy * 2.; << 693 y = y * halfy * 2.;
906 y = y - halfy; << 694 y = y - halfy;
907 z = z * halfz * 2.; << 695 z = z * halfz * 2.;
908 z = z - halfz; << 696 z = z - halfz;
909 << 697 // Pick a side first
910 // Pick a side first << 698 if(testrand < Probs[0])
911 // << 699 {
912 if(testrand < Probs[0]) << 700 // side is +x
913 { << 701 x = halfx + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
914 // side is +x << 702 y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
915 << 703 z = z;
916 x = halfx + z*std::tan(ParTheta)*std::co << 704 // Cosine-law
917 y = y + z*std::tan(ParTheta)*std::sin(Pa << 705 G4ThreeVector xdash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
918 // z = z; << 706 halfz*std::tan(ParTheta)*std::sin(ParPhi),
919 << 707 halfz/std::cos(ParPhi));
920 // Cosine-law << 708 G4ThreeVector ydash(halfy*std::tan(ParAlpha), -halfy, 0.0);
921 // << 709 xdash = xdash.unit();
922 G4ThreeVector xdash(halfz*std::tan(ParTh << 710 ydash = ydash.unit();
923 halfz*std::tan(ParTh << 711 G4ThreeVector zdash = xdash.cross(ydash);
924 halfz/std::cos(ParPh << 712 SideRefVec1 = xdash.unit();
925 G4ThreeVector ydash(halfy*std::tan(ParAl << 713 SideRefVec2 = ydash.unit();
926 xdash = xdash.unit(); << 714 SideRefVec3 = zdash.unit();
927 ydash = ydash.unit(); << 715 }
928 G4ThreeVector zdash = xdash.cross(ydash) << 716 else if(testrand >= Probs[0] && testrand < Probs[1])
929 td.CSideRefVec1 = xdash.unit(); << 717 {
930 td.CSideRefVec2 = ydash.unit(); << 718 // side is -x
931 td.CSideRefVec3 = zdash.unit(); << 719 x = -halfx + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
932 } << 720 y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
933 else if(testrand >= Probs[0] && testrand < << 721 z = z;
934 { << 722 // Cosine-law
935 // side is -x << 723 G4ThreeVector xdash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
936 << 724 halfz*std::tan(ParTheta)*std::sin(ParPhi),
937 x = -halfx + z*std::tan(ParTheta)*std::c << 725 halfz/std::cos(ParPhi));
938 y = y + z*std::tan(ParTheta)*std::sin(Pa << 726 G4ThreeVector ydash(halfy*std::tan(ParAlpha), halfy, 0.0);
939 // z = z; << 727 xdash = xdash.unit();
940 << 728 ydash = ydash.unit();
941 // Cosine-law << 729 G4ThreeVector zdash = xdash.cross(ydash);
942 // << 730 SideRefVec1 = xdash.unit();
943 G4ThreeVector xdash(halfz*std::tan(ParTh << 731 SideRefVec2 = ydash.unit();
944 halfz*std::tan(ParTh << 732 SideRefVec3 = zdash.unit();
945 halfz/std::cos(ParPh << 733 }
946 G4ThreeVector ydash(halfy*std::tan(ParAl << 734 else if(testrand >= Probs[1] && testrand < Probs[2])
947 xdash = xdash.unit(); << 735 {
948 ydash = ydash.unit(); << 736 // side is +y
949 G4ThreeVector zdash = xdash.cross(ydash) << 737 x = x + z*std::tan(ParTheta)*std::cos(ParPhi) + halfy*std::tan(ParAlpha);
950 td.CSideRefVec1 = xdash.unit(); << 738 y = halfy + z*std::tan(ParTheta)*std::sin(ParPhi);
951 td.CSideRefVec2 = ydash.unit(); << 739 z = z;
952 td.CSideRefVec3 = zdash.unit(); << 740 // Cosine-law
953 } << 741 G4ThreeVector ydash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
954 else if(testrand >= Probs[1] && testrand < << 742 halfz*std::tan(ParTheta)*std::sin(ParPhi),
955 { << 743 halfz/std::cos(ParPhi));
956 // side is +y << 744 ydash = ydash.unit();
957 << 745 G4ThreeVector xdash = Roty.cross(ydash);
958 x = x + z*std::tan(ParTheta)*std::cos(Pa << 746 G4ThreeVector zdash = xdash.cross(ydash);
959 y = halfy + z*std::tan(ParTheta)*std::si << 747 SideRefVec1 = xdash.unit();
960 // z = z; << 748 SideRefVec2 = -ydash.unit();
961 << 749 SideRefVec3 = -zdash.unit();
962 // Cosine-law << 750 }
963 // << 751 else if(testrand >= Probs[2] && testrand < Probs[3])
964 G4ThreeVector ydash(halfz*std::tan(ParTh << 752 {
965 halfz*std::tan(ParTh << 753 // side is -y
966 halfz/std::cos(ParPh << 754 x = x + z*std::tan(ParTheta)*std::cos(ParPhi) - halfy*std::tan(ParAlpha);
967 ydash = ydash.unit(); << 755 y = -halfy + z*std::tan(ParTheta)*std::sin(ParPhi);
968 G4ThreeVector xdash = Roty.cross(ydash); << 756 z = z;
969 G4ThreeVector zdash = xdash.cross(ydash) << 757 // Cosine-law
970 td.CSideRefVec1 = xdash.unit(); << 758 G4ThreeVector ydash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
971 td.CSideRefVec2 = -ydash.unit(); << 759 halfz*std::tan(ParTheta)*std::sin(ParPhi),
972 td.CSideRefVec3 = -zdash.unit(); << 760 halfz/std::cos(ParPhi));
973 } << 761 ydash = ydash.unit();
974 else if(testrand >= Probs[2] && testrand < << 762 G4ThreeVector xdash = Roty.cross(ydash);
975 { << 763 G4ThreeVector zdash = xdash.cross(ydash);
976 // side is -y << 764 SideRefVec1 = xdash.unit();
977 << 765 SideRefVec2 = ydash.unit();
978 x = x + z*std::tan(ParTheta)*std::cos(Pa << 766 SideRefVec3 = zdash.unit();
979 y = -halfy + z*std::tan(ParTheta)*std::s << 767 }
980 // z = z; << 768 else if(testrand >= Probs[3] && testrand < Probs[4])
981 << 769 {
982 // Cosine-law << 770 // side is +z
983 // << 771 z = halfz;
984 G4ThreeVector ydash(halfz*std::tan(ParTh << 772 y = y + halfz*std::sin(ParPhi)*std::tan(ParTheta);
985 halfz*std::tan(ParTh << 773 x = x + halfz*std::cos(ParPhi)*std::tan(ParTheta) + y*std::tan(ParAlpha);
986 halfz/std::cos(ParPh << 774 // Cosine-law
987 ydash = ydash.unit(); << 775 SideRefVec1 = Rotx;
988 G4ThreeVector xdash = Roty.cross(ydash); << 776 SideRefVec2 = Roty;
989 G4ThreeVector zdash = xdash.cross(ydash) << 777 SideRefVec3 = Rotz;
990 td.CSideRefVec1 = xdash.unit(); << 778 }
991 td.CSideRefVec2 = ydash.unit(); << 779 else if(testrand >= Probs[4] && testrand < Probs[5])
992 td.CSideRefVec3 = zdash.unit(); << 780 {
993 } << 781 // side is -z
994 else if(testrand >= Probs[3] && testrand < << 782 z = -halfz;
995 { << 783 y = y - halfz*std::sin(ParPhi)*std::tan(ParTheta);
996 // side is +z << 784 x = x - halfz*std::cos(ParPhi)*std::tan(ParTheta) + y*std::tan(ParAlpha);
997 << 785 // Cosine-law
998 z = halfz; << 786 SideRefVec1 = Rotx;
999 y = y + halfz*std::sin(ParPhi)*std::tan( << 787 SideRefVec2 = -Roty;
1000 x = x + halfz*std::cos(ParPhi)*std::tan << 788 SideRefVec3 = -Rotz;
1001 << 789 }
1002 // Cosine-law << 790 else
1003 // << 791 {
1004 td.CSideRefVec1 = Rotx; << 792 G4cout << "Error: testrand out of range" << G4endl;
1005 td.CSideRefVec2 = Roty; << 793 if(verbosityLevel >= 1)
1006 td.CSideRefVec3 = Rotz; << 794 G4cout << "testrand=" << testrand << " Probs[5]=" << Probs[5] <<G4endl;
1007 } << 795 }
1008 else if(testrand >= Probs[4] && testrand << 796
1009 { << 797 RandPos.setX(x);
1010 // side is -z << 798 RandPos.setY(y);
1011 << 799 RandPos.setZ(z);
1012 z = -halfz; <<
1013 y = y - halfz*std::sin(ParPhi)*std::tan <<
1014 x = x - halfz*std::cos(ParPhi)*std::tan <<
1015 <<
1016 // Cosine-law <<
1017 // <<
1018 td.CSideRefVec1 = Rotx; <<
1019 td.CSideRefVec2 = -Roty; <<
1020 td.CSideRefVec3 = -Rotz; <<
1021 } <<
1022 else <<
1023 { <<
1024 G4cout << "Error: testrand out of range <<
1025 if(verbosityLevel >= 1) <<
1026 { <<
1027 G4cout << "testrand=" << testrand << <<
1028 } <<
1029 } 800 }
1030 801
1031 RandPos.setX(x); <<
1032 RandPos.setY(y); <<
1033 RandPos.setZ(z); <<
1034 } <<
1035 <<
1036 // Apply Rotation Matrix 802 // Apply Rotation Matrix
1037 // x * Rotx, y * Roty and z * Rotz 803 // x * Rotx, y * Roty and z * Rotz
1038 // <<
1039 if(verbosityLevel == 2) 804 if(verbosityLevel == 2)
1040 { <<
1041 G4cout << "Raw position " << RandPos << G 805 G4cout << "Raw position " << RandPos << G4endl;
1042 } << 806
1043 x=(RandPos.x()*Rotx.x())+(RandPos.y()*Roty. 807 x=(RandPos.x()*Rotx.x())+(RandPos.y()*Roty.x())+(RandPos.z()*Rotz.x());
1044 y=(RandPos.x()*Rotx.y())+(RandPos.y()*Roty. 808 y=(RandPos.x()*Rotx.y())+(RandPos.y()*Roty.y())+(RandPos.z()*Rotz.y());
1045 z=(RandPos.x()*Rotx.z())+(RandPos.y()*Roty. 809 z=(RandPos.x()*Rotx.z())+(RandPos.y()*Roty.z())+(RandPos.z()*Rotz.z());
1046 810
1047 RandPos.setX(x); 811 RandPos.setX(x);
1048 RandPos.setY(y); 812 RandPos.setY(y);
1049 RandPos.setZ(z); 813 RandPos.setZ(z);
1050 814
1051 // Translate 815 // Translate
1052 // << 816 particle_position = CentreCoords + RandPos;
1053 pos = CentreCoords + RandPos; <<
1054 817
1055 if(verbosityLevel >= 1) 818 if(verbosityLevel >= 1)
1056 { <<
1057 if(verbosityLevel == 2) <<
1058 { 819 {
1059 G4cout << "Rotated position " << RandPo << 820 if(verbosityLevel == 2)
>> 821 G4cout << "Rotated position " << RandPos << G4endl;
>> 822 G4cout << "Rotated and translated position " << particle_position << G4endl;
1060 } 823 }
1061 G4cout << "Rotated and translated positio <<
1062 } <<
1063 if(verbosityLevel == 2) 824 if(verbosityLevel == 2)
1064 { << 825 {
1065 G4cout << "Reference vectors for cosine-l << 826 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
1066 << " " << td.CSideRefVec2 << " " < << 827 }
1067 } <<
1068 } 828 }
1069 829
1070 void G4SPSPosDistribution::GeneratePointsInVo << 830 void G4SPSPosDistribution::GeneratePointsInVolume()
1071 { 831 {
1072 G4ThreeVector RandPos; 832 G4ThreeVector RandPos;
1073 G4double tempx, tempy, tempz; 833 G4double tempx, tempy, tempz;
1074 G4double x, y, z; 834 G4double x, y, z;
1075 G4double expression; <<
1076 x = y = z = 0.; 835 x = y = z = 0.;
1077 <<
1078 if(SourcePosType != "Volume" && verbosityLe 836 if(SourcePosType != "Volume" && verbosityLevel >= 1)
1079 { <<
1080 G4cout << "Error SourcePosType not Volume 837 G4cout << "Error SourcePosType not Volume" << G4endl;
1081 } << 838 //Private method to create points in a volume
1082 <<
1083 // Private method to create points in a vol <<
1084 // <<
1085 if(Shape == "Sphere") 839 if(Shape == "Sphere")
1086 { << 840 {
1087 x = Radius*2.; << 841 x = Radius*2.;
1088 y = Radius*2.; << 842 y = Radius*2.;
1089 z = Radius*2.; << 843 z = Radius*2.;
1090 while(((x*x)+(y*y)+(z*z)) > (Radius*Radiu << 844 while(((x*x)+(y*y)+(z*z)) > (Radius*Radius))
1091 { << 845 {
1092 x = PosRndm->GenRandX(); << 846 x = posRndm->GenRandX();
1093 y = PosRndm->GenRandY(); << 847 y = posRndm->GenRandY();
1094 z = PosRndm->GenRandZ(); << 848 z = posRndm->GenRandZ();
1095 << 849
1096 x = (x*2.*Radius) - Radius; << 850 x = (x*2.*Radius) - Radius;
1097 y = (y*2.*Radius) - Radius; << 851 y = (y*2.*Radius) - Radius;
1098 z = (z*2.*Radius) - Radius; << 852 z = (z*2.*Radius) - Radius;
>> 853 }
1099 } 854 }
1100 } <<
1101 else if(Shape == "Ellipsoid") 855 else if(Shape == "Ellipsoid")
1102 { << 856 {
1103 G4double temp; << 857 G4double temp;
1104 temp = 100.; << 858 temp = 100.;
1105 while(temp > 1.) << 859 while(temp > 1.)
1106 { << 860 {
1107 x = PosRndm->GenRandX(); << 861 x = posRndm->GenRandX();
1108 y = PosRndm->GenRandY(); << 862 y = posRndm->GenRandY();
1109 z = PosRndm->GenRandZ(); << 863 z = posRndm->GenRandZ();
1110 << 864
1111 x = (x*2.*halfx) - halfx; << 865 x = (x*2.*halfx) - halfx;
1112 y = (y*2.*halfy) - halfy; << 866 y = (y*2.*halfy) - halfy;
1113 z = (z*2.*halfz) - halfz; << 867 z = (z*2.*halfz) - halfz;
1114 << 868
1115 temp = ((x*x)/(halfx*halfx))+((y*y)/(ha << 869 temp = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy))
>> 870 + ((z*z)/(halfz*halfz));
>> 871 }
1116 } 872 }
1117 } <<
1118 else if(Shape == "Cylinder") 873 else if(Shape == "Cylinder")
1119 { << 874 {
1120 x = Radius*2.; << 875 x = Radius*2.;
1121 y = Radius*2.; << 876 y = Radius*2.;
1122 while(((x*x)+(y*y)) > (Radius*Radius)) << 877 while(((x*x)+(y*y)) > (Radius*Radius))
1123 { << 878 {
1124 x = PosRndm->GenRandX(); << 879 x = posRndm->GenRandX();
1125 y = PosRndm->GenRandY(); << 880 y = posRndm->GenRandY();
1126 z = PosRndm->GenRandZ(); << 881 z = posRndm->GenRandZ();
1127 << 882
1128 x = (x*2.*Radius) - Radius; << 883 x = (x*2.*Radius) - Radius;
1129 y = (y*2.*Radius) - Radius; << 884 y = (y*2.*Radius) - Radius;
1130 z = (z*2.*halfz) - halfz; << 885 z = (z*2.*halfz) - halfz;
>> 886 }
1131 } 887 }
1132 } << 888 else if(Shape == "Para")
1133 else if(Shape == "EllipticCylinder") << 889 {
1134 { << 890 x = posRndm->GenRandX();
1135 expression = 20.; << 891 y = posRndm->GenRandY();
1136 while(expression > 1.) << 892 z = posRndm->GenRandZ();
1137 { <<
1138 x = PosRndm->GenRandX(); <<
1139 y = PosRndm->GenRandY(); <<
1140 z = PosRndm->GenRandZ(); <<
1141 <<
1142 x = (x*2.*halfx) - halfx; 893 x = (x*2.*halfx) - halfx;
1143 y = (y*2.*halfy) - halfy; 894 y = (y*2.*halfy) - halfy;
1144 z = (z*2.*halfz) - halfz; 895 z = (z*2.*halfz) - halfz;
1145 << 896 x = x + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
1146 expression = ((x*x)/(halfx*halfx)) + (( << 897 y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
>> 898 z = z;
1147 } 899 }
1148 } <<
1149 else if(Shape == "Para") <<
1150 { <<
1151 x = PosRndm->GenRandX(); <<
1152 y = PosRndm->GenRandY(); <<
1153 z = PosRndm->GenRandZ(); <<
1154 x = (x*2.*halfx) - halfx; <<
1155 y = (y*2.*halfy) - halfy; <<
1156 z = (z*2.*halfz) - halfz; <<
1157 x = x + z*std::tan(ParTheta)*std::cos(Par <<
1158 y = y + z*std::tan(ParTheta)*std::sin(Par <<
1159 // z = z; <<
1160 } <<
1161 else 900 else
1162 { << 901 G4cout << "Error: Volume Shape Doesnt Exist" << G4endl;
1163 G4cout << "Error: Volume Shape does not e <<
1164 } <<
1165 902
1166 RandPos.setX(x); 903 RandPos.setX(x);
1167 RandPos.setY(y); 904 RandPos.setY(y);
1168 RandPos.setZ(z); 905 RandPos.setZ(z);
1169 906
1170 // Apply Rotation Matrix 907 // Apply Rotation Matrix
1171 // x * Rotx, y * Roty and z * Rotz 908 // x * Rotx, y * Roty and z * Rotz
1172 // <<
1173 tempx = (x * Rotx.x()) + (y * Roty.x()) + ( 909 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
1174 tempy = (x * Rotx.y()) + (y * Roty.y()) + ( 910 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
1175 tempz = (x * Rotx.z()) + (y * Roty.z()) + ( 911 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
1176 912
1177 RandPos.setX(tempx); 913 RandPos.setX(tempx);
1178 RandPos.setY(tempy); 914 RandPos.setY(tempy);
1179 RandPos.setZ(tempz); 915 RandPos.setZ(tempz);
1180 916
1181 // Translate 917 // Translate
1182 // << 918 particle_position = CentreCoords + RandPos;
1183 pos = CentreCoords + RandPos; <<
1184 919
1185 if(verbosityLevel == 2) 920 if(verbosityLevel == 2)
1186 { << 921 {
1187 G4cout << "Raw position " << x << "," << << 922 G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
1188 G4cout << "Rotated position " << RandPos << 923 G4cout << "Rotated position " << RandPos << G4endl;
1189 } << 924 }
1190 if(verbosityLevel >= 1) 925 if(verbosityLevel >= 1)
1191 { << 926 G4cout << "Rotated and translated position " << particle_position << G4endl;
1192 G4cout << "Rotated and translated positio <<
1193 } <<
1194 927
1195 // Cosine-law (not a good idea to use this 928 // Cosine-law (not a good idea to use this here)
1196 // <<
1197 G4ThreeVector zdash(tempx,tempy,tempz); 929 G4ThreeVector zdash(tempx,tempy,tempz);
1198 zdash = zdash.unit(); 930 zdash = zdash.unit();
1199 G4ThreeVector xdash = Rotz.cross(zdash); 931 G4ThreeVector xdash = Rotz.cross(zdash);
1200 G4ThreeVector ydash = xdash.cross(zdash); 932 G4ThreeVector ydash = xdash.cross(zdash);
1201 thread_data_t& td = ThreadData.Get(); << 933 SideRefVec1 = xdash.unit();
1202 td.CSideRefVec1 = xdash.unit(); << 934 SideRefVec2 = ydash.unit();
1203 td.CSideRefVec2 = ydash.unit(); << 935 SideRefVec3 = zdash.unit();
1204 td.CSideRefVec3 = zdash.unit(); << 936
1205 if(verbosityLevel == 2) 937 if(verbosityLevel == 2)
1206 { << 938 {
1207 G4cout << "Reference vectors for cosine-l << 939 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
1208 << " " << td.CSideRefVec2 << " " < << 940 }
1209 } <<
1210 } 941 }
1211 942
1212 G4bool G4SPSPosDistribution::IsSourceConfined << 943 G4bool G4SPSPosDistribution::IsSourceConfined()
1213 { 944 {
1214 // Method to check point is within the volu 945 // Method to check point is within the volume specified
1215 << 946 if(Confine == false)
1216 if(!Confine) <<
1217 { <<
1218 G4cout << "Error: Confine is false" << G4 947 G4cout << "Error: Confine is false" << G4endl;
1219 } <<
1220 G4ThreeVector null(0.,0.,0.); 948 G4ThreeVector null(0.,0.,0.);
1221 G4ThreeVector* ptr = &null; << 949 G4ThreeVector *ptr;
>> 950 ptr = &null;
1222 951
1223 // Check position is within VolName, if so << 952 // Check particle_position is within VolName, if so true,
1224 // << 953 // else false
1225 G4VPhysicalVolume* theVolume; << 954 G4VPhysicalVolume *theVolume;
1226 G4Navigator* gNavigator = G4TransportationM << 955 theVolume=gNavigator->LocateGlobalPointAndSetup(particle_position,ptr,true);
1227 ->GetNavigatorForTr <<
1228 theVolume = gNavigator->LocateGlobalPointAn <<
1229 if(theVolume == nullptr) { return false; } <<
1230 G4String theVolName = theVolume->GetName(); 956 G4String theVolName = theVolume->GetName();
1231 <<
1232 if(theVolName == VolName) 957 if(theVolName == VolName)
1233 { <<
1234 if(verbosityLevel >= 1) <<
1235 { 958 {
1236 G4cout << "Particle is in volume " << V << 959 if(verbosityLevel >= 1)
>> 960 G4cout << "Particle is in volume " << VolName << G4endl;
>> 961 return(true);
1237 } 962 }
1238 return true; << 963 else
1239 } << 964 return(false);
1240 <<
1241 return false; <<
1242 <<
1243 } 965 }
1244 966
1245 G4ThreeVector G4SPSPosDistribution::GenerateO 967 G4ThreeVector G4SPSPosDistribution::GenerateOne()
1246 { 968 {
1247 G4ThreeVector localP; << 969 //
1248 G4bool srcconf = false; 970 G4bool srcconf = false;
1249 G4int LoopCount = 0; 971 G4int LoopCount = 0;
1250 while(!srcconf) << 972 while(srcconf == false)
1251 { << 973 {
1252 if(SourcePosType == "Point") << 974 if(SourcePosType == "Point")
1253 GeneratePointSource(localP); << 975 GeneratePointSource();
1254 else if(SourcePosType == "Beam") << 976 else if(SourcePosType == "Beam")
1255 GeneratePointsInBeam(localP); << 977 GeneratePointsInBeam();
1256 else if(SourcePosType == "Plane") << 978 else if(SourcePosType == "Plane")
1257 GeneratePointsInPlane(localP); << 979 GeneratePointsInPlane();
1258 else if(SourcePosType == "Surface") << 980 else if(SourcePosType == "Surface")
1259 GeneratePointsOnSurface(localP); << 981 GeneratePointsOnSurface();
1260 else if(SourcePosType == "Volume") << 982 else if(SourcePosType == "Volume")
1261 GeneratePointsInVolume(localP); << 983 GeneratePointsInVolume();
1262 else << 984 else
1263 { << 985 {
1264 G4ExceptionDescription msg; << 986 G4cout << "Error: SourcePosType undefined" << G4endl;
1265 msg << "Error: SourcePosType undefined\ << 987 G4cout << "Generating point source" << G4endl;
1266 msg << "Generating point source\n"; << 988 GeneratePointSource();
1267 G4Exception("G4SPSPosDistribution::Gene << 989 }
1268 "G4GPS001", JustWarning, ms << 990 if(Confine == true)
1269 GeneratePointSource(localP); << 991 {
1270 } << 992 srcconf = IsSourceConfined();
1271 if(Confine) << 993 // if source in confined srcconf = true terminating the loop
1272 { << 994 // if source isnt confined srcconf = false and loop continues
1273 srcconf = IsSourceConfined(localP); << 995 }
1274 // if source in confined srcconf = true << 996 else if(Confine == false)
1275 // if source isnt confined srcconf = fa << 997 srcconf = true; // terminate loop
1276 } << 998 LoopCount++;
1277 else if(!Confine) << 999 if(LoopCount == 100000)
1278 { << 1000 {
1279 srcconf = true; // terminate loop << 1001 G4cout << "*************************************" << G4endl;
1280 } << 1002 G4cout << "LoopCount = 100000" << G4endl;
1281 ++LoopCount; << 1003 G4cout << "Either the source distribution >> confinement" << G4endl;
1282 if(LoopCount == 100000) << 1004 G4cout << "or any confining volume may not overlap with" << G4endl;
1283 { << 1005 G4cout << "the source distribution or any confining volumes" << G4endl;
1284 G4ExceptionDescription msg; << 1006 G4cout << "may not exist"<< G4endl;
1285 msg << "LoopCount = 100000\n"; << 1007 G4cout << "If you have set confine then this will be ignored" <<G4endl;
1286 msg << "Either the source distribution << 1008 G4cout << "for this event." << G4endl;
1287 msg << "or any confining volume may not << 1009 G4cout << "*************************************" << G4endl;
1288 msg << "the source distribution or any << 1010 srcconf = true; //Avoids an infinite loop
1289 msg << "may not exist\n"<< G4endl; << 1011 }
1290 msg << "If you have set confine then th <<
1291 msg << "for this event.\n" << G4endl; <<
1292 G4Exception("G4SPSPosDistribution::Gene <<
1293 "G4GPS001", JustWarning, ms <<
1294 srcconf = true; // Avoids an infinite l <<
1295 } 1012 }
1296 } << 1013 return particle_position;
1297 ThreadData.Get().CParticlePos = localP; <<
1298 return localP; <<
1299 } 1014 }
>> 1015
>> 1016
>> 1017
1300 1018
1301 1019