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