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