| 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/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(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(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(sqrt((x*x) + (y*y)) > Radius || 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->GenRandTheta();
539 theta = std::acos(1. - 2.*theta); // theta << 440 phi = posRndm->GenRandPhi();
540 phi = phi * 2. * pi; << 441
541 << 442 theta = acos(1. - 2.*theta); // theta isotropic
542 x = Radius * std::sin(theta) * std::cos(ph << 443 phi = phi * 2. * pi;
543 y = Radius * std::sin(theta) * std::sin(ph << 444 tantheta = tan(theta);
544 z = Radius * std::cos(theta); <<
545 <<
546 RandPos.setX(x); <<
547 RandPos.setY(y); <<
548 RandPos.setZ(z); <<
549 <<
550 // Cosine-law (not a good idea to use this <<
551 // <<
552 G4ThreeVector zdash(x,y,z); <<
553 zdash = zdash.unit(); <<
554 G4ThreeVector xdash = Rotz.cross(zdash); <<
555 G4ThreeVector ydash = xdash.cross(zdash); <<
556 td.CSideRefVec1 = xdash.unit(); <<
557 td.CSideRefVec2 = ydash.unit(); <<
558 td.CSideRefVec3 = zdash.unit(); <<
559 } <<
560 else if(Shape == "Ellipsoid") <<
561 { <<
562 G4double minphi, maxphi, middlephi; <<
563 G4double answer, constant; <<
564 <<
565 constant = pi/(halfx*halfx) + pi/(halfy*ha <<
566 <<
567 // Simplified approach <<
568 // <<
569 theta = PosRndm->GenRandPosTheta(); <<
570 phi = PosRndm->GenRandPosPhi(); <<
571 445
572 theta = std::acos(1. - 2.*theta); << 446 x = Radius * sin(theta) * cos(phi);
573 minphi = 0.; << 447 y = Radius * sin(theta) * sin(phi);
574 maxphi = twopi; << 448 z = Radius * cos(theta);
575 while(maxphi-minphi > 0.) <<
576 { <<
577 middlephi = (maxphi+minphi)/2.; <<
578 answer = (1./(halfx*halfx))*(middlephi/2 <<
579 + (1./(halfy*halfy))*(middlephi/2 <<
580 + middlephi/(halfz*halfz); <<
581 answer = answer/constant; <<
582 if(answer > phi) maxphi = middlephi; <<
583 if(answer < phi) minphi = middlephi; <<
584 if(std::fabs(answer-phi) <= 0.00001) <<
585 { <<
586 minphi = maxphi +1; <<
587 phi = middlephi; <<
588 } <<
589 } <<
590 <<
591 // x,y and z form a unit vector. Put this <<
592 // <<
593 x = std::sin(theta)*std::cos(phi); <<
594 y = std::sin(theta)*std::sin(phi); <<
595 z = std::cos(theta); <<
596 <<
597 G4double lhs; <<
598 <<
599 // Solve for x <<
600 // <<
601 G4double numYinX = y/x; <<
602 G4double numZinX = z/x; <<
603 G4double tempxvar; <<
604 tempxvar = 1./(halfx*halfx)+(numYinX*numYi <<
605 + (numZinX*numZinX)/(halfz*halfz) <<
606 tempxvar = 1./tempxvar; <<
607 G4double coordx = std::sqrt(tempxvar); <<
608 <<
609 // Solve for y <<
610 // <<
611 G4double numXinY = x/y; <<
612 G4double numZinY = z/y; <<
613 G4double tempyvar; <<
614 tempyvar = (numXinY*numXinY)/(halfx*halfx) <<
615 + (numZinY*numZinY)/(halfz*halfz) <<
616 tempyvar = 1./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 449
633 if(std::fabs(lhs-1.) > 0.001 && verbosityL << 450 RandPos.setX(x);
634 { << 451 RandPos.setY(y);
635 G4cout << "Error: theta, phi not really << 452 RandPos.setZ(z);
636 } <<
637 453
638 // coordx, coordy and coordz are all posit << 454 // Cosine-law (not a good idea to use this here)
639 // << 455 G4ThreeVector zdash(x,y,z);
640 G4double TestRandVar = G4UniformRand(); << 456 zdash = zdash.unit();
641 if(TestRandVar > 0.5) << 457 G4ThreeVector xdash = Rotz.cross(zdash);
642 { << 458 G4ThreeVector ydash = xdash.cross(zdash);
643 coordx = -coordx; << 459 SideRefVec1 = xdash.unit();
644 } << 460 SideRefVec2 = ydash.unit();
645 TestRandVar = G4UniformRand(); << 461 SideRefVec3 = zdash.unit();
646 if(TestRandVar > 0.5) <<
647 { <<
648 coordy = -coordy; <<
649 } 462 }
650 TestRandVar = G4UniformRand(); << 463 else if(Shape == "Ellipsoid")
651 if(TestRandVar > 0.5) <<
652 { 464 {
653 coordz = -coordz; << 465 G4double theta, phi, minphi, maxphi, middlephi;
654 } << 466 G4double answer, constant;
655 467
656 RandPos.setX(coordx); << 468 constant = pi/(halfx*halfx) + pi/(halfy*halfy) +
657 RandPos.setY(coordy); << 469 twopi/(halfz*halfz);
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 470
685 prob1 = AreaTop / AreaTotal; << 471 // simplified approach
686 prob2 = AreaBot / AreaTotal; << 472 theta = posRndm->GenRandTheta();
687 prob3 = 1.00 - prob1 - prob2; << 473 phi = posRndm->GenRandPhi();
688 if(std::fabs(prob3 - (AreaLat/AreaTotal)) << 474
689 { << 475 theta = acos(1. - 2.*theta);
690 if(verbosityLevel >= 1) << 476 minphi = 0.;
691 { << 477 maxphi = twopi;
692 G4cout << AreaLat/AreaTotal << " " << << 478 while(maxphi-minphi > 0.)
693 } << 479 {
694 G4cout << "Error in prob3" << G4endl; << 480 middlephi = (maxphi+minphi)/2.;
695 } << 481 answer = (1./(halfx*halfx))*(middlephi/2. + sin(2*middlephi)/4.)
696 << 482 + (1./(halfy*halfy))*(middlephi/2. - sin(2*middlephi)/4.)
697 // Decide surface to calculate point on. << 483 + middlephi/(halfz*halfz);
>> 484 answer = answer/constant;
>> 485 if(answer > phi) maxphi = middlephi;
>> 486 if(answer < phi) minphi = middlephi;
>> 487 if(fabs(answer-phi) <= 0.00001)
>> 488 {
>> 489 minphi = maxphi +1;
>> 490 phi = middlephi;
>> 491 }
>> 492 }
>> 493
>> 494 x = sin(theta)*cos(phi);
>> 495 y = sin(theta)*sin(phi);
>> 496 z = cos(theta);
>> 497 // x,y and z form a unit vector. Put this onto the ellipse.
>> 498 G4double lhs;
>> 499 // solve for x
>> 500 G4double numYinX = y/x;
>> 501 G4double numZinX = z/x;
>> 502 G4double tempxvar;
>> 503 tempxvar= 1./(halfx*halfx)+(numYinX*numYinX)/(halfy*halfy)
>> 504 + (numZinX*numZinX)/(halfz*halfz);
698 505
699 testrand = G4UniformRand(); << 506 tempxvar = 1./tempxvar;
700 if(testrand <= prob1) // Point on Top sur << 507 G4double coordx = sqrt(tempxvar);
701 { << 508
702 z = halfz; << 509 //solve for y
703 x = Radius + 100.; << 510 G4double numXinY = x/y;
704 y = Radius + 100.; << 511 G4double numZinY = z/y;
705 while(((x*x)+(y*y)) > (Radius*Radius)) << 512 G4double tempyvar;
706 { << 513 tempyvar=(numXinY*numXinY)/(halfx*halfx)+1./(halfy*halfy)
707 x = PosRndm->GenRandX(); << 514 +(numZinY*numZinY)/(halfz*halfz);
708 y = PosRndm->GenRandY(); << 515 tempyvar = 1./tempyvar;
709 << 516 G4double coordy = sqrt(tempyvar);
710 x = x * 2. * Radius; << 517
711 y = y * 2. * Radius; << 518 //solve for z
712 x = x - Radius; << 519 G4double numXinZ = x/z;
713 y = y - Radius; << 520 G4double numYinZ = y/z;
714 } << 521 G4double tempzvar;
715 // Cosine law << 522 tempzvar=(numXinZ*numXinZ)/(halfx*halfx)
716 // << 523 +(numYinZ*numYinZ)/(halfy*halfy)+1./(halfz*halfz);
717 td.CSideRefVec1 = Rotx; << 524 tempzvar = 1./tempzvar;
718 td.CSideRefVec2 = Roty; << 525 G4double coordz = sqrt(tempzvar);
719 td.CSideRefVec3 = Rotz; << 526
720 } << 527 lhs = sqrt((coordx*coordx)/(halfx*halfx) +
721 else if((testrand > prob1) && (testrand <= << 528 (coordy*coordy)/(halfy*halfy) +
722 { // Point on Bot << 529 (coordz*coordz)/(halfz*halfz));
723 z = -halfz; << 530
724 x = Radius + 100.; << 531 if(fabs(lhs-1.) > 0.001 && verbosityLevel >= 1)
725 y = Radius + 100.; << 532 G4cout << "Error: theta, phi not really on ellipsoid" << G4endl;
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 <<
749 x = Radius * std::cos(rand); <<
750 y = Radius * std::sin(rand); <<
751 533
752 z = PosRndm->GenRandZ(); << 534 // coordx, coordy and coordz are all positive
>> 535 G4double TestRandVar = G4UniformRand();
>> 536 if(TestRandVar > 0.5)
>> 537 {
>> 538 coordx = -coordx;
>> 539 }
>> 540 TestRandVar = G4UniformRand();
>> 541 if(TestRandVar > 0.5)
>> 542 {
>> 543 coordy = -coordy;
>> 544 }
>> 545 TestRandVar = G4UniformRand();
>> 546 if(TestRandVar > 0.5)
>> 547 {
>> 548 coordz = -coordz;
>> 549 }
>> 550
>> 551 RandPos.setX(coordx);
>> 552 RandPos.setY(coordy);
>> 553 RandPos.setZ(coordz);
753 554
754 z = z * 2. *halfz; << 555 // Cosine-law (not a good idea to use this here)
755 z = z - halfz; << 556 G4ThreeVector zdash(coordx,coordy,coordz);
756 <<
757 // Cosine law <<
758 // <<
759 G4ThreeVector zdash(x,y,0.); <<
760 zdash = zdash.unit(); 557 zdash = zdash.unit();
761 G4ThreeVector xdash = Rotz.cross(zdash); 558 G4ThreeVector xdash = Rotz.cross(zdash);
762 G4ThreeVector ydash = xdash.cross(zdash) 559 G4ThreeVector ydash = xdash.cross(zdash);
763 td.CSideRefVec1 = xdash.unit(); << 560 SideRefVec1 = xdash.unit();
764 td.CSideRefVec2 = ydash.unit(); << 561 SideRefVec2 = ydash.unit();
765 td.CSideRefVec3 = zdash.unit(); << 562 SideRefVec3 = zdash.unit();
766 } <<
767 else <<
768 { <<
769 G4cout << "Error: testrand " << testrand <<
770 } 563 }
771 << 564 else if(Shape == "Cylinder")
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 } <<
845 else if(testrand > (prob1+prob2)) // Poin <<
846 { 565 {
847 G4double rand; << 566 G4double AreaTop, AreaBot, AreaLat;
848 << 567 G4double AreaTotal, prob1, prob2, prob3;
849 rand = PosRndm->GenRandPosPhi(); << 568 G4double testrand;
850 rand = rand * 2. * pi; << 569
851 << 570 // User giver Radius and z-half length
852 x = halfx * std::cos(rand); << 571 // Calculate surface areas, maybe move this to
853 y = halfy * std::sin(rand); << 572 // a different routine.
854 << 573
855 z = PosRndm->GenRandZ(); << 574 AreaTop = pi * Radius * Radius;
>> 575 AreaBot = AreaTop;
>> 576 AreaLat = 2. * pi * Radius * 2. * halfz;
>> 577 AreaTotal = AreaTop + AreaBot + AreaLat;
>> 578
>> 579 prob1 = AreaTop / AreaTotal;
>> 580 prob2 = AreaBot / AreaTotal;
>> 581 prob3 = 1.00 - prob1 - prob2;
>> 582 if(fabs(prob3 - (AreaLat/AreaTotal)) >= 0.001)
>> 583 {
>> 584 if(verbosityLevel >= 1)
>> 585 G4cout << AreaLat/AreaTotal << " " << prob3<<G4endl;
>> 586 G4cout << "Error in prob3" << G4endl;
>> 587 }
>> 588
>> 589 // Decide surface to calculate point on.
>> 590
>> 591 testrand = G4UniformRand();
>> 592 if(testrand <= prob1)
>> 593 {
>> 594 //Point on Top surface
>> 595 z = halfz;
>> 596 x = Radius + 100.;
>> 597 y = Radius + 100.;
>> 598 while(((x*x)+(y*y)) > (Radius*Radius))
>> 599 {
>> 600 x = posRndm->GenRandX();
>> 601 y = posRndm->GenRandY();
>> 602
>> 603 x = x * 2. * Radius;
>> 604 y = y * 2. * Radius;
>> 605 x = x - Radius;
>> 606 y = y - Radius;
>> 607 }
>> 608 // Cosine law
>> 609 SideRefVec1 = Rotx;
>> 610 SideRefVec2 = Roty;
>> 611 SideRefVec3 = Rotz;
>> 612 }
>> 613 else if((testrand > prob1) && (testrand <= (prob1 + prob2)))
>> 614 {
>> 615 //Point on Bottom surface
>> 616 z = -halfz;
>> 617 x = Radius + 100.;
>> 618 y = Radius + 100.;
>> 619 while(((x*x)+(y*y)) > (Radius*Radius))
>> 620 {
>> 621 x = posRndm->GenRandX();
>> 622 y = posRndm->GenRandY();
>> 623
>> 624 x = x * 2. * Radius;
>> 625 y = y * 2. * Radius;
>> 626 x = x - Radius;
>> 627 y = y - Radius;
>> 628 }
>> 629 // Cosine law
>> 630 SideRefVec1 = Rotx;
>> 631 SideRefVec2 = -Roty;
>> 632 SideRefVec3 = -Rotz;
>> 633 }
>> 634 else if(testrand > (prob1+prob2))
>> 635 {
>> 636 G4double rand;
>> 637 //Point on Lateral Surface
>> 638
>> 639 rand = G4UniformRand();
>> 640 rand = rand * 2. * pi;
>> 641
>> 642 x = Radius * cos(rand);
>> 643 y = Radius * sin(rand);
>> 644
>> 645 z = posRndm->GenRandZ();
>> 646
>> 647 z = z * 2. * halfz;
>> 648 z = z - halfz;
>> 649
>> 650 // Cosine law
>> 651 G4ThreeVector zdash(x,y,0.);
>> 652 zdash = zdash.unit();
>> 653 G4ThreeVector xdash = Rotz.cross(zdash);
>> 654 G4ThreeVector ydash = xdash.cross(zdash);
>> 655 SideRefVec1 = xdash.unit();
>> 656 SideRefVec2 = ydash.unit();
>> 657 SideRefVec3 = zdash.unit();
>> 658 }
>> 659 else
>> 660 G4cout << "Error: testrand " << testrand << G4endl;
>> 661
>> 662 RandPos.setX(x);
>> 663 RandPos.setY(y);
>> 664 RandPos.setZ(z);
856 665
857 z = (z * 2. * halfz) - halfz; <<
858 } <<
859 else <<
860 { <<
861 G4cout << "Error: testrand " << testrand <<
862 } 666 }
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") 667 else if(Shape == "Para")
879 { << 668 {
880 // Right Parallelepiped. << 669 G4double testrand;
881 // User gives x,y,z half lengths and ParAl << 670 //Right Parallelepiped.
882 // ParTheta and ParPhi << 671 // User gives x,y,z half lengths and ParAlpha
883 // +x = <1, -x >1 & <2, +y >2 & <3, -y >3 << 672 // ParTheta and ParPhi
884 // +z >4 & < 5, -z >5 &<6 << 673 // +x = <1, -x >1 & <2, +y >2 & <3, -y >3 &<4
885 << 674 // +z >4 & < 5, -z >5 &<6.
886 G4double testrand = G4UniformRand(); << 675 testrand = G4UniformRand();
887 G4double AreaX = halfy * halfz * 4.; << 676 G4double AreaX = halfy * halfz * 4.;
888 G4double AreaY = halfx * halfz * 4.; << 677 G4double AreaY = halfx * halfz * 4.;
889 G4double AreaZ = halfx * halfy * 4.; << 678 G4double AreaZ = halfx * halfy * 4.;
890 G4double AreaTotal = 2*(AreaX + AreaY + Ar << 679 G4double AreaTotal = 2*(AreaX + AreaY + AreaZ);
891 G4double Probs[6]; << 680 G4double Probs[6];
892 Probs[0] = AreaX/AreaTotal; << 681 Probs[0] = AreaX/AreaTotal;
893 Probs[1] = Probs[0] + AreaX/AreaTotal; << 682 Probs[1] = Probs[0] + AreaX/AreaTotal;
894 Probs[2] = Probs[1] + AreaY/AreaTotal; << 683 Probs[2] = Probs[1] + AreaY/AreaTotal;
895 Probs[3] = Probs[2] + AreaY/AreaTotal; << 684 Probs[3] = Probs[2] + AreaY/AreaTotal;
896 Probs[4] = Probs[3] + AreaZ/AreaTotal; << 685 Probs[4] = Probs[3] + AreaZ/AreaTotal;
897 Probs[5] = Probs[4] + AreaZ/AreaTotal; << 686 Probs[5] = Probs[4] + AreaZ/AreaTotal;
898 687
899 x = PosRndm->GenRandX(); << 688 x = posRndm->GenRandX();
900 y = PosRndm->GenRandY(); << 689 y = posRndm->GenRandY();
901 z = PosRndm->GenRandZ(); << 690 z = posRndm->GenRandZ();
902 691
903 x = x * halfx * 2.; << 692 x = x * halfx * 2.;
904 x = x - halfx; << 693 x = x - halfx;
905 y = y * halfy * 2.; << 694 y = y * halfy * 2.;
906 y = y - halfy; << 695 y = y - halfy;
907 z = z * halfz * 2.; << 696 z = z * halfz * 2.;
908 z = z - halfz; << 697 z = z - halfz;
909 << 698 // Pick a side first
910 // Pick a side first << 699 if(testrand < Probs[0])
911 // << 700 {
912 if(testrand < Probs[0]) << 701 // side is +x
913 { << 702 x = halfx + z*tan(ParTheta)*cos(ParPhi) + y*tan(ParAlpha);
914 // side is +x << 703 y = y + z*tan(ParTheta)*sin(ParPhi);
915 << 704 z = z;
916 x = halfx + z*std::tan(ParTheta)*std::co << 705 // Cosine-law
917 y = y + z*std::tan(ParTheta)*std::sin(Pa << 706 G4ThreeVector xdash(halfz*tan(ParTheta)*cos(ParPhi),
918 // z = z; << 707 halfz*tan(ParTheta)*sin(ParPhi),
919 << 708 halfz/cos(ParPhi));
920 // Cosine-law << 709 G4ThreeVector ydash(halfy*tan(ParAlpha), -halfy, 0.0);
921 // << 710 xdash = xdash.unit();
922 G4ThreeVector xdash(halfz*std::tan(ParTh << 711 ydash = ydash.unit();
923 halfz*std::tan(ParTh << 712 G4ThreeVector zdash = xdash.cross(ydash);
924 halfz/std::cos(ParPh << 713 SideRefVec1 = xdash.unit();
925 G4ThreeVector ydash(halfy*std::tan(ParAl << 714 SideRefVec2 = ydash.unit();
926 xdash = xdash.unit(); << 715 SideRefVec3 = zdash.unit();
927 ydash = ydash.unit(); << 716 }
928 G4ThreeVector zdash = xdash.cross(ydash) << 717 else if(testrand >= Probs[0] && testrand < Probs[1])
929 td.CSideRefVec1 = xdash.unit(); << 718 {
930 td.CSideRefVec2 = ydash.unit(); << 719 // side is -x
931 td.CSideRefVec3 = zdash.unit(); << 720 x = -halfx + z*tan(ParTheta)*cos(ParPhi) + y*tan(ParAlpha);
932 } << 721 y = y + z*tan(ParTheta)*sin(ParPhi);
933 else if(testrand >= Probs[0] && testrand < << 722 z = z;
934 { << 723 // Cosine-law
935 // side is -x << 724 G4ThreeVector xdash(halfz*tan(ParTheta)*cos(ParPhi),
936 << 725 halfz*tan(ParTheta)*sin(ParPhi),
937 x = -halfx + z*std::tan(ParTheta)*std::c << 726 halfz/cos(ParPhi));
938 y = y + z*std::tan(ParTheta)*std::sin(Pa << 727 G4ThreeVector ydash(halfy*tan(ParAlpha), halfy, 0.0);
939 // z = z; << 728 xdash = xdash.unit();
940 << 729 ydash = ydash.unit();
941 // Cosine-law << 730 G4ThreeVector zdash = xdash.cross(ydash);
942 // << 731 SideRefVec1 = xdash.unit();
943 G4ThreeVector xdash(halfz*std::tan(ParTh << 732 SideRefVec2 = ydash.unit();
944 halfz*std::tan(ParTh << 733 SideRefVec3 = zdash.unit();
945 halfz/std::cos(ParPh << 734 }
946 G4ThreeVector ydash(halfy*std::tan(ParAl << 735 else if(testrand >= Probs[1] && testrand < Probs[2])
947 xdash = xdash.unit(); << 736 {
948 ydash = ydash.unit(); << 737 // side is +y
949 G4ThreeVector zdash = xdash.cross(ydash) << 738 x = x + z*tan(ParTheta)*cos(ParPhi) + halfy*tan(ParAlpha);
950 td.CSideRefVec1 = xdash.unit(); << 739 y = halfy + z*tan(ParTheta)*sin(ParPhi);
951 td.CSideRefVec2 = ydash.unit(); << 740 z = z;
952 td.CSideRefVec3 = zdash.unit(); << 741 // Cosine-law
953 } << 742 G4ThreeVector ydash(halfz*tan(ParTheta)*cos(ParPhi),
954 else if(testrand >= Probs[1] && testrand < << 743 halfz*tan(ParTheta)*sin(ParPhi),
955 { << 744 halfz/cos(ParPhi));
956 // side is +y << 745 ydash = ydash.unit();
957 << 746 G4ThreeVector xdash = Roty.cross(ydash);
958 x = x + z*std::tan(ParTheta)*std::cos(Pa << 747 G4ThreeVector zdash = xdash.cross(ydash);
959 y = halfy + z*std::tan(ParTheta)*std::si << 748 SideRefVec1 = xdash.unit();
960 // z = z; << 749 SideRefVec2 = -ydash.unit();
961 << 750 SideRefVec3 = -zdash.unit();
962 // Cosine-law << 751 }
963 // << 752 else if(testrand >= Probs[2] && testrand < Probs[3])
964 G4ThreeVector ydash(halfz*std::tan(ParTh << 753 {
965 halfz*std::tan(ParTh << 754 // side is -y
966 halfz/std::cos(ParPh << 755 x = x + z*tan(ParTheta)*cos(ParPhi) - halfy*tan(ParAlpha);
967 ydash = ydash.unit(); << 756 y = -halfy + z*tan(ParTheta)*sin(ParPhi);
968 G4ThreeVector xdash = Roty.cross(ydash); << 757 z = z;
969 G4ThreeVector zdash = xdash.cross(ydash) << 758 // Cosine-law
970 td.CSideRefVec1 = xdash.unit(); << 759 G4ThreeVector ydash(halfz*tan(ParTheta)*cos(ParPhi),
971 td.CSideRefVec2 = -ydash.unit(); << 760 halfz*tan(ParTheta)*sin(ParPhi),
972 td.CSideRefVec3 = -zdash.unit(); << 761 halfz/cos(ParPhi));
973 } << 762 ydash = ydash.unit();
974 else if(testrand >= Probs[2] && testrand < << 763 G4ThreeVector xdash = Roty.cross(ydash);
975 { << 764 G4ThreeVector zdash = xdash.cross(ydash);
976 // side is -y << 765 SideRefVec1 = xdash.unit();
977 << 766 SideRefVec2 = ydash.unit();
978 x = x + z*std::tan(ParTheta)*std::cos(Pa << 767 SideRefVec3 = zdash.unit();
979 y = -halfy + z*std::tan(ParTheta)*std::s << 768 }
980 // z = z; << 769 else if(testrand >= Probs[3] && testrand < Probs[4])
981 << 770 {
982 // Cosine-law << 771 // side is +z
983 // << 772 z = halfz;
984 G4ThreeVector ydash(halfz*std::tan(ParTh << 773 y = y + halfz*sin(ParPhi)*tan(ParTheta);
985 halfz*std::tan(ParTh << 774 x = x + halfz*cos(ParPhi)*tan(ParTheta) + y*tan(ParAlpha);
986 halfz/std::cos(ParPh << 775 // Cosine-law
987 ydash = ydash.unit(); << 776 SideRefVec1 = Rotx;
988 G4ThreeVector xdash = Roty.cross(ydash); << 777 SideRefVec2 = Roty;
989 G4ThreeVector zdash = xdash.cross(ydash) << 778 SideRefVec3 = Rotz;
990 td.CSideRefVec1 = xdash.unit(); << 779 }
991 td.CSideRefVec2 = ydash.unit(); << 780 else if(testrand >= Probs[4] && testrand < Probs[5])
992 td.CSideRefVec3 = zdash.unit(); << 781 {
993 } << 782 // side is -z
994 else if(testrand >= Probs[3] && testrand < << 783 z = -halfz;
995 { << 784 y = y - halfz*sin(ParPhi)*tan(ParTheta);
996 // side is +z << 785 x = x - halfz*cos(ParPhi)*tan(ParTheta) + y*tan(ParAlpha);
997 << 786 // Cosine-law
998 z = halfz; << 787 SideRefVec1 = Rotx;
999 y = y + halfz*std::sin(ParPhi)*std::tan( << 788 SideRefVec2 = -Roty;
1000 x = x + halfz*std::cos(ParPhi)*std::tan << 789 SideRefVec3 = -Rotz;
1001 << 790 }
1002 // Cosine-law << 791 else
1003 // << 792 {
1004 td.CSideRefVec1 = Rotx; << 793 G4cout << "Error: testrand out of range" << G4endl;
1005 td.CSideRefVec2 = Roty; << 794 if(verbosityLevel >= 1)
1006 td.CSideRefVec3 = Rotz; << 795 G4cout << "testrand=" << testrand << " Probs[5]=" << Probs[5] <<G4endl;
1007 } << 796 }
1008 else if(testrand >= Probs[4] && testrand << 797
1009 { << 798 RandPos.setX(x);
1010 // side is -z << 799 RandPos.setY(y);
1011 << 800 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 } 801 }
1030 802
1031 RandPos.setX(x); <<
1032 RandPos.setY(y); <<
1033 RandPos.setZ(z); <<
1034 } <<
1035 <<
1036 // Apply Rotation Matrix 803 // Apply Rotation Matrix
1037 // x * Rotx, y * Roty and z * Rotz 804 // x * Rotx, y * Roty and z * Rotz
1038 // <<
1039 if(verbosityLevel == 2) 805 if(verbosityLevel == 2)
1040 { <<
1041 G4cout << "Raw position " << RandPos << G 806 G4cout << "Raw position " << RandPos << G4endl;
1042 } << 807
1043 x=(RandPos.x()*Rotx.x())+(RandPos.y()*Roty. 808 x=(RandPos.x()*Rotx.x())+(RandPos.y()*Roty.x())+(RandPos.z()*Rotz.x());
1044 y=(RandPos.x()*Rotx.y())+(RandPos.y()*Roty. 809 y=(RandPos.x()*Rotx.y())+(RandPos.y()*Roty.y())+(RandPos.z()*Rotz.y());
1045 z=(RandPos.x()*Rotx.z())+(RandPos.y()*Roty. 810 z=(RandPos.x()*Rotx.z())+(RandPos.y()*Roty.z())+(RandPos.z()*Rotz.z());
1046 811
1047 RandPos.setX(x); 812 RandPos.setX(x);
1048 RandPos.setY(y); 813 RandPos.setY(y);
1049 RandPos.setZ(z); 814 RandPos.setZ(z);
1050 815
1051 // Translate 816 // Translate
1052 // << 817 particle_position = CentreCoords + RandPos;
1053 pos = CentreCoords + RandPos; <<
1054 818
1055 if(verbosityLevel >= 1) 819 if(verbosityLevel >= 1)
1056 { <<
1057 if(verbosityLevel == 2) <<
1058 { 820 {
1059 G4cout << "Rotated position " << RandPo << 821 if(verbosityLevel == 2)
>> 822 G4cout << "Rotated position " << RandPos << G4endl;
>> 823 G4cout << "Rotated and translated position " << particle_position << G4endl;
1060 } 824 }
1061 G4cout << "Rotated and translated positio <<
1062 } <<
1063 if(verbosityLevel == 2) 825 if(verbosityLevel == 2)
1064 { << 826 {
1065 G4cout << "Reference vectors for cosine-l << 827 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
1066 << " " << td.CSideRefVec2 << " " < << 828 }
1067 } <<
1068 } 829 }
1069 830
1070 void G4SPSPosDistribution::GeneratePointsInVo << 831 void G4SPSPosDistribution::GeneratePointsInVolume()
1071 { 832 {
1072 G4ThreeVector RandPos; 833 G4ThreeVector RandPos;
1073 G4double tempx, tempy, tempz; 834 G4double tempx, tempy, tempz;
1074 G4double x, y, z; 835 G4double x, y, z;
1075 G4double expression; <<
1076 x = y = z = 0.; 836 x = y = z = 0.;
1077 <<
1078 if(SourcePosType != "Volume" && verbosityLe 837 if(SourcePosType != "Volume" && verbosityLevel >= 1)
1079 { <<
1080 G4cout << "Error SourcePosType not Volume 838 G4cout << "Error SourcePosType not Volume" << G4endl;
1081 } << 839 //Private method to create points in a volume
1082 <<
1083 // Private method to create points in a vol <<
1084 // <<
1085 if(Shape == "Sphere") 840 if(Shape == "Sphere")
1086 { << 841 {
1087 x = Radius*2.; << 842 x = Radius*2.;
1088 y = Radius*2.; << 843 y = Radius*2.;
1089 z = Radius*2.; << 844 z = Radius*2.;
1090 while(((x*x)+(y*y)+(z*z)) > (Radius*Radiu << 845 while(((x*x)+(y*y)+(z*z)) > (Radius*Radius))
1091 { << 846 {
1092 x = PosRndm->GenRandX(); << 847 x = posRndm->GenRandX();
1093 y = PosRndm->GenRandY(); << 848 y = posRndm->GenRandY();
1094 z = PosRndm->GenRandZ(); << 849 z = posRndm->GenRandZ();
1095 << 850
1096 x = (x*2.*Radius) - Radius; << 851 x = (x*2.*Radius) - Radius;
1097 y = (y*2.*Radius) - Radius; << 852 y = (y*2.*Radius) - Radius;
1098 z = (z*2.*Radius) - Radius; << 853 z = (z*2.*Radius) - Radius;
>> 854 }
1099 } 855 }
1100 } <<
1101 else if(Shape == "Ellipsoid") 856 else if(Shape == "Ellipsoid")
1102 { << 857 {
1103 G4double temp; << 858 G4double temp;
1104 temp = 100.; << 859 temp = 100.;
1105 while(temp > 1.) << 860 while(temp > 1.)
1106 { << 861 {
1107 x = PosRndm->GenRandX(); << 862 x = posRndm->GenRandX();
1108 y = PosRndm->GenRandY(); << 863 y = posRndm->GenRandY();
1109 z = PosRndm->GenRandZ(); << 864 z = posRndm->GenRandZ();
1110 << 865
1111 x = (x*2.*halfx) - halfx; << 866 x = (x*2.*halfx) - halfx;
1112 y = (y*2.*halfy) - halfy; << 867 y = (y*2.*halfy) - halfy;
1113 z = (z*2.*halfz) - halfz; << 868 z = (z*2.*halfz) - halfz;
1114 << 869
1115 temp = ((x*x)/(halfx*halfx))+((y*y)/(ha << 870 temp = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy))
>> 871 + ((z*z)/(halfz*halfz));
>> 872 }
1116 } 873 }
1117 } <<
1118 else if(Shape == "Cylinder") 874 else if(Shape == "Cylinder")
1119 { << 875 {
1120 x = Radius*2.; << 876 x = Radius*2.;
1121 y = Radius*2.; << 877 y = Radius*2.;
1122 while(((x*x)+(y*y)) > (Radius*Radius)) << 878 while(((x*x)+(y*y)) > (Radius*Radius))
1123 { << 879 {
1124 x = PosRndm->GenRandX(); << 880 x = posRndm->GenRandX();
1125 y = PosRndm->GenRandY(); << 881 y = posRndm->GenRandY();
1126 z = PosRndm->GenRandZ(); << 882 z = posRndm->GenRandZ();
1127 << 883
1128 x = (x*2.*Radius) - Radius; << 884 x = (x*2.*Radius) - Radius;
1129 y = (y*2.*Radius) - Radius; << 885 y = (y*2.*Radius) - Radius;
1130 z = (z*2.*halfz) - halfz; << 886 z = (z*2.*halfz) - halfz;
>> 887 }
1131 } 888 }
1132 } << 889 else if(Shape == "Para")
1133 else if(Shape == "EllipticCylinder") << 890 {
1134 { << 891 x = posRndm->GenRandX();
1135 expression = 20.; << 892 y = posRndm->GenRandY();
1136 while(expression > 1.) << 893 z = posRndm->GenRandZ();
1137 { <<
1138 x = PosRndm->GenRandX(); <<
1139 y = PosRndm->GenRandY(); <<
1140 z = PosRndm->GenRandZ(); <<
1141 <<
1142 x = (x*2.*halfx) - halfx; 894 x = (x*2.*halfx) - halfx;
1143 y = (y*2.*halfy) - halfy; 895 y = (y*2.*halfy) - halfy;
1144 z = (z*2.*halfz) - halfz; 896 z = (z*2.*halfz) - halfz;
1145 << 897 x = x + z*tan(ParTheta)*cos(ParPhi) + y*tan(ParAlpha);
1146 expression = ((x*x)/(halfx*halfx)) + (( << 898 y = y + z*tan(ParTheta)*sin(ParPhi);
>> 899 z = z;
1147 } 900 }
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 901 else
1162 { << 902 G4cout << "Error: Volume Shape Doesnt Exist" << G4endl;
1163 G4cout << "Error: Volume Shape does not e <<
1164 } <<
1165 903
1166 RandPos.setX(x); 904 RandPos.setX(x);
1167 RandPos.setY(y); 905 RandPos.setY(y);
1168 RandPos.setZ(z); 906 RandPos.setZ(z);
1169 907
1170 // Apply Rotation Matrix 908 // Apply Rotation Matrix
1171 // x * Rotx, y * Roty and z * Rotz 909 // x * Rotx, y * Roty and z * Rotz
1172 // <<
1173 tempx = (x * Rotx.x()) + (y * Roty.x()) + ( 910 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
1174 tempy = (x * Rotx.y()) + (y * Roty.y()) + ( 911 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
1175 tempz = (x * Rotx.z()) + (y * Roty.z()) + ( 912 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
1176 913
1177 RandPos.setX(tempx); 914 RandPos.setX(tempx);
1178 RandPos.setY(tempy); 915 RandPos.setY(tempy);
1179 RandPos.setZ(tempz); 916 RandPos.setZ(tempz);
1180 917
1181 // Translate 918 // Translate
1182 // << 919 particle_position = CentreCoords + RandPos;
1183 pos = CentreCoords + RandPos; <<
1184 920
1185 if(verbosityLevel == 2) 921 if(verbosityLevel == 2)
1186 { << 922 {
1187 G4cout << "Raw position " << x << "," << << 923 G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
1188 G4cout << "Rotated position " << RandPos << 924 G4cout << "Rotated position " << RandPos << G4endl;
1189 } << 925 }
1190 if(verbosityLevel >= 1) 926 if(verbosityLevel >= 1)
1191 { << 927 G4cout << "Rotated and translated position " << particle_position << G4endl;
1192 G4cout << "Rotated and translated positio <<
1193 } <<
1194 928
1195 // Cosine-law (not a good idea to use this 929 // Cosine-law (not a good idea to use this here)
1196 // <<
1197 G4ThreeVector zdash(tempx,tempy,tempz); 930 G4ThreeVector zdash(tempx,tempy,tempz);
1198 zdash = zdash.unit(); 931 zdash = zdash.unit();
1199 G4ThreeVector xdash = Rotz.cross(zdash); 932 G4ThreeVector xdash = Rotz.cross(zdash);
1200 G4ThreeVector ydash = xdash.cross(zdash); 933 G4ThreeVector ydash = xdash.cross(zdash);
1201 thread_data_t& td = ThreadData.Get(); << 934 SideRefVec1 = xdash.unit();
1202 td.CSideRefVec1 = xdash.unit(); << 935 SideRefVec2 = ydash.unit();
1203 td.CSideRefVec2 = ydash.unit(); << 936 SideRefVec3 = zdash.unit();
1204 td.CSideRefVec3 = zdash.unit(); << 937
1205 if(verbosityLevel == 2) 938 if(verbosityLevel == 2)
1206 { << 939 {
1207 G4cout << "Reference vectors for cosine-l << 940 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
1208 << " " << td.CSideRefVec2 << " " < << 941 }
1209 } <<
1210 } 942 }
1211 943
1212 G4bool G4SPSPosDistribution::IsSourceConfined << 944 G4bool G4SPSPosDistribution::IsSourceConfined()
1213 { 945 {
1214 // Method to check point is within the volu 946 // Method to check point is within the volume specified
1215 << 947 if(Confine == false)
1216 if(!Confine) <<
1217 { <<
1218 G4cout << "Error: Confine is false" << G4 948 G4cout << "Error: Confine is false" << G4endl;
1219 } <<
1220 G4ThreeVector null(0.,0.,0.); 949 G4ThreeVector null(0.,0.,0.);
1221 G4ThreeVector* ptr = &null; << 950 G4ThreeVector *ptr;
>> 951 ptr = &null;
1222 952
1223 // Check position is within VolName, if so << 953 // Check particle_position is within VolName, if so true,
1224 // << 954 // else false
1225 G4VPhysicalVolume* theVolume; << 955 G4VPhysicalVolume *theVolume;
1226 G4Navigator* gNavigator = G4TransportationM << 956 theVolume=gNavigator->LocateGlobalPointAndSetup(particle_position,ptr,true);
1227 ->GetNavigatorForTr <<
1228 theVolume = gNavigator->LocateGlobalPointAn <<
1229 if(theVolume == nullptr) { return false; } <<
1230 G4String theVolName = theVolume->GetName(); 957 G4String theVolName = theVolume->GetName();
1231 <<
1232 if(theVolName == VolName) 958 if(theVolName == VolName)
1233 { <<
1234 if(verbosityLevel >= 1) <<
1235 { 959 {
1236 G4cout << "Particle is in volume " << V << 960 if(verbosityLevel >= 1)
>> 961 G4cout << "Particle is in volume " << VolName << G4endl;
>> 962 return(true);
1237 } 963 }
1238 return true; << 964 else
1239 } << 965 return(false);
1240 <<
1241 return false; <<
1242 <<
1243 } 966 }
1244 967
1245 G4ThreeVector G4SPSPosDistribution::GenerateO 968 G4ThreeVector G4SPSPosDistribution::GenerateOne()
1246 { 969 {
1247 G4ThreeVector localP; << 970 //
1248 G4bool srcconf = false; 971 G4bool srcconf = false;
1249 G4int LoopCount = 0; 972 G4int LoopCount = 0;
1250 while(!srcconf) << 973 while(srcconf == false)
1251 { << 974 {
1252 if(SourcePosType == "Point") << 975 if(SourcePosType == "Point")
1253 GeneratePointSource(localP); << 976 GeneratePointSource();
1254 else if(SourcePosType == "Beam") << 977 else if(SourcePosType == "Beam")
1255 GeneratePointsInBeam(localP); << 978 GeneratePointsInBeam();
1256 else if(SourcePosType == "Plane") << 979 else if(SourcePosType == "Plane")
1257 GeneratePointsInPlane(localP); << 980 GeneratePointsInPlane();
1258 else if(SourcePosType == "Surface") << 981 else if(SourcePosType == "Surface")
1259 GeneratePointsOnSurface(localP); << 982 GeneratePointsOnSurface();
1260 else if(SourcePosType == "Volume") << 983 else if(SourcePosType == "Volume")
1261 GeneratePointsInVolume(localP); << 984 GeneratePointsInVolume();
1262 else << 985 else
1263 { << 986 {
1264 G4ExceptionDescription msg; << 987 G4cout << "Error: SourcePosType undefined" << G4endl;
1265 msg << "Error: SourcePosType undefined\ << 988 G4cout << "Generating point source" << G4endl;
1266 msg << "Generating point source\n"; << 989 GeneratePointSource();
1267 G4Exception("G4SPSPosDistribution::Gene << 990 }
1268 "G4GPS001", JustWarning, ms << 991 if(Confine == true)
1269 GeneratePointSource(localP); << 992 {
1270 } << 993 srcconf = IsSourceConfined();
1271 if(Confine) << 994 // if source in confined srcconf = true terminating the loop
1272 { << 995 // if source isnt confined srcconf = false and loop continues
1273 srcconf = IsSourceConfined(localP); << 996 }
1274 // if source in confined srcconf = true << 997 else if(Confine == false)
1275 // if source isnt confined srcconf = fa << 998 srcconf = true; // terminate loop
1276 } << 999 LoopCount++;
1277 else if(!Confine) << 1000 if(LoopCount == 100000)
1278 { << 1001 {
1279 srcconf = true; // terminate loop << 1002 G4cout << "*************************************" << G4endl;
1280 } << 1003 G4cout << "LoopCount = 100000" << G4endl;
1281 ++LoopCount; << 1004 G4cout << "Either the source distribution >> confinement" << G4endl;
1282 if(LoopCount == 100000) << 1005 G4cout << "or any confining volume may not overlap with" << G4endl;
1283 { << 1006 G4cout << "the source distribution or any confining volumes" << G4endl;
1284 G4ExceptionDescription msg; << 1007 G4cout << "may not exist"<< G4endl;
1285 msg << "LoopCount = 100000\n"; << 1008 G4cout << "If you have set confine then this will be ignored" <<G4endl;
1286 msg << "Either the source distribution << 1009 G4cout << "for this event." << G4endl;
1287 msg << "or any confining volume may not << 1010 G4cout << "*************************************" << G4endl;
1288 msg << "the source distribution or any << 1011 srcconf = true; //Avoids an infinite loop
1289 msg << "may not exist\n"<< G4endl; << 1012 }
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 } 1013 }
1296 } << 1014 return particle_position;
1297 ThreadData.Get().CParticlePos = localP; <<
1298 return localP; <<
1299 } 1015 }
>> 1016
>> 1017
>> 1018
1300 1019
1301 1020