source: trunk/source/event/src/G4SPSPosDistribution.cc@ 1337

Last change on this file since 1337 was 816, checked in by garnier, 17 years ago

import all except CVS
File size: 27.6 KB
RevLine 
[816]1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
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///////////////////////////////////////////////////////////////////////////////
46//
47#include "Randomize.hh"
48#include "G4TransportationManager.hh"
49#include "G4VPhysicalVolume.hh"
50#include "G4PhysicalVolumeStore.hh"
51
52#include "G4SPSPosDistribution.hh"
53
54G4SPSPosDistribution::G4SPSPosDistribution()
55{
56
57 // Initialise all variables
58 // Position distribution Variables
59
60 SourcePosType = "Point";
61 Shape = "NULL";
62 halfx = 0.;
63 halfy = 0.;
64 halfz = 0.;
65 Radius = 0.;
66 Radius0 = 0.;
67 SR = 0.;
68 SX = 0.;
69 SY = 0.;
70 ParAlpha = 0.;
71 ParTheta = 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
78 VolName = "NULL";
79 SideRefVec1 = CLHEP::HepXHat; // x-axis
80 SideRefVec2 = CLHEP::HepYHat; // y-axis
81 SideRefVec3 = CLHEP::HepZHat; // z-axis
82 verbosityLevel = 0 ;
83 gNavigator = G4TransportationManager::GetTransportationManager()
84 ->GetNavigatorForTracking();
85}
86
87G4SPSPosDistribution::~G4SPSPosDistribution()
88{
89}
90
91void G4SPSPosDistribution::SetPosDisType(G4String PosType)
92{
93 SourcePosType = PosType;
94}
95
96void G4SPSPosDistribution::SetPosDisShape(G4String shapeType)
97{
98 Shape = shapeType;
99}
100
101void G4SPSPosDistribution::SetCentreCoords(G4ThreeVector coordsOfCentre)
102{
103 CentreCoords = coordsOfCentre;
104}
105
106void G4SPSPosDistribution::SetPosRot1(G4ThreeVector posrot1)
107{
108 // This should be x'
109 Rotx = posrot1;
110 if(verbosityLevel == 2)
111 {
112 G4cout << "Vector x' " << Rotx << G4endl;
113 }
114 GenerateRotationMatrices();
115}
116
117void G4SPSPosDistribution::SetPosRot2(G4ThreeVector posrot2)
118{
119 // This is a vector in the plane x'y' but need not
120 // be y'
121 Roty = posrot2;
122 if(verbosityLevel == 2)
123 {
124 G4cout << "The vector in the x'-y' plane " << Roty << G4endl;
125 }
126 GenerateRotationMatrices();
127}
128
129void G4SPSPosDistribution::SetHalfX(G4double xhalf)
130{
131 halfx = xhalf;
132}
133
134void G4SPSPosDistribution::SetHalfY(G4double yhalf)
135{
136 halfy = yhalf;
137}
138
139void G4SPSPosDistribution::SetHalfZ(G4double zhalf)
140{
141 halfz = zhalf;
142}
143
144void G4SPSPosDistribution::SetRadius(G4double rad)
145{
146 Radius = rad;
147}
148
149void G4SPSPosDistribution::SetRadius0(G4double rad)
150{
151 Radius0 = rad;
152}
153
154void G4SPSPosDistribution::SetBeamSigmaInR(G4double r)
155{
156 SR = r;
157 SX = SY = r/std::sqrt(2.);
158}
159
160void G4SPSPosDistribution::SetBeamSigmaInX(G4double r)
161{
162 SX = r;
163}
164
165void G4SPSPosDistribution::SetBeamSigmaInY(G4double r)
166{
167 SY = r;
168}
169
170void G4SPSPosDistribution::SetParAlpha(G4double paralp)
171{
172 ParAlpha = paralp;
173}
174
175void G4SPSPosDistribution::SetParTheta(G4double parthe)
176{
177 ParTheta = parthe;
178}
179
180void G4SPSPosDistribution::SetParPhi(G4double parphi)
181{
182 ParPhi = parphi;
183}
184
185void G4SPSPosDistribution::GenerateRotationMatrices()
186{
187 // This takes in 2 vectors, x' and one in the plane x'-y',
188 // and from these takes a cross product to calculate z'.
189 // Then a cross product is taken between x' and z' to give
190 // y'.
191 Rotx = Rotx.unit(); // x'
192 Roty = Roty.unit(); // vector in x'y' plane
193 Rotz = Rotx.cross(Roty); // z'
194 Roty = Rotz.cross(Rotx); // y'
195 if(verbosityLevel == 2)
196 {
197 G4cout << "The new axes, x', y', z' " << Rotx << " " << Roty << " " << Rotz << G4endl;
198 }
199}
200
201void G4SPSPosDistribution::ConfineSourceToVolume(G4String Vname)
202{
203 VolName = Vname;
204 if(verbosityLevel == 2)
205 G4cout << VolName << G4endl;
206 G4VPhysicalVolume *tempPV = NULL;
207 G4PhysicalVolumeStore *PVStore = 0;
208 G4String theRequiredVolumeName = VolName;
209 PVStore = G4PhysicalVolumeStore::GetInstance();
210 G4int i = 0;
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++;}
221 }
222 // found = true then the volume exists else it doesnt.
223 if(found == true)
224 {
225 if(verbosityLevel >= 1)
226 G4cout << "Volume " << VolName << " exists" << G4endl;
227 Confine = true;
228 }
229 else
230 {
231 G4cout << " **** Error: Volume does not exist **** " << G4endl;
232 G4cout << " Ignoring confine condition" << G4endl;
233 Confine = false;
234 VolName = "NULL";
235 }
236
237}
238
239void G4SPSPosDistribution::GeneratePointSource()
240{
241 // Generates Points given the point source.
242 if(SourcePosType == "Point")
243 particle_position = CentreCoords;
244 else
245 if(verbosityLevel >= 1)
246 G4cout << "Error SourcePosType is not set to Point" << G4endl;
247}
248
249void G4SPSPosDistribution::GeneratePointsInBeam()
250{
251 G4double x, y, z;
252
253 G4ThreeVector RandPos;
254 G4double tempx, tempy, tempz;
255 z = 0.;
256
257 // Private Method to create points in a plane
258 if(Shape == "Circle")
259 {
260 x = Radius + 100.;
261 y = Radius + 100.;
262 while(std::sqrt((x*x) + (y*y)) > Radius)
263 {
264 x = posRndm->GenRandX();
265 y = posRndm->GenRandY();
266
267 x = (x*2.*Radius) - Radius;
268 y = (y*2.*Radius) - Radius;
269 }
270 x += G4RandGauss::shoot(0.0,SX) ;
271 y += G4RandGauss::shoot(0.0,SY) ;
272 }
273 else
274 {
275 // all other cases default to Rectangle case
276 x = posRndm->GenRandX();
277 y = posRndm->GenRandY();
278 x = (x*2.*halfx) - halfx;
279 y = (y*2.*halfy) - halfy;
280 x += G4RandGauss::shoot(0.0,SX);
281 y += G4RandGauss::shoot(0.0,SY);
282 }
283 // Apply Rotation Matrix
284 // x * Rotx, y * Roty and z * Rotz
285 if(verbosityLevel >= 2)
286 {
287 G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
288 }
289 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
290 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
291 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
292
293 RandPos.setX(tempx);
294 RandPos.setY(tempy);
295 RandPos.setZ(tempz);
296
297 // Translate
298 particle_position = CentreCoords + RandPos;
299 if(verbosityLevel >= 1)
300 {
301 if(verbosityLevel >= 2)
302 {
303 G4cout << "Rotated Position " << RandPos << G4endl;
304 }
305 G4cout << "Rotated and Translated position " << particle_position << G4endl;
306 }
307}
308
309void G4SPSPosDistribution::GeneratePointsInPlane()
310{
311 G4double x, y, z;
312 G4double expression;
313 G4ThreeVector RandPos;
314 G4double tempx, tempy, tempz;
315 x = y = z = 0.;
316
317 if(SourcePosType != "Plane" && verbosityLevel >= 1)
318 G4cout << "Error: SourcePosType is not Plane" << G4endl;
319
320 // Private Method to create points in a plane
321 if(Shape == "Circle")
322 {
323 x = Radius + 100.;
324 y = Radius + 100.;
325 while(std::sqrt((x*x) + (y*y)) > Radius)
326 {
327 x = posRndm->GenRandX();
328 y = posRndm->GenRandY();
329
330 x = (x*2.*Radius) - Radius;
331 y = (y*2.*Radius) - Radius;
332 }
333 }
334 else if(Shape == "Annulus")
335 {
336 x = Radius + 100.;
337 y = Radius + 100.;
338 while(std::sqrt((x*x) + (y*y)) > Radius || std::sqrt((x*x) + (y*y)) < Radius0 )
339 {
340 x = posRndm->GenRandX();
341 y = posRndm->GenRandY();
342
343 x = (x*2.*Radius) - Radius;
344 y = (y*2.*Radius) - Radius;
345 }
346 }
347 else if(Shape == "Ellipse")
348 {
349 expression = 20.;
350 while(expression > 1.)
351 {
352 x = posRndm->GenRandX();
353 y = posRndm->GenRandY();
354
355 x = (x*2.*halfx) - halfx;
356 y = (y*2.*halfy) - halfy;
357
358 expression = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy));
359 }
360 }
361 else if(Shape == "Square")
362 {
363 x = posRndm->GenRandX();
364 y = posRndm->GenRandY();
365 x = (x*2.*halfx) - halfx;
366 y = (y*2.*halfy) - halfy;
367 }
368 else if(Shape == "Rectangle")
369 {
370 x = posRndm->GenRandX();
371 y = posRndm->GenRandY();
372 x = (x*2.*halfx) - halfx;
373 y = (y*2.*halfy) - halfy;
374 }
375 else
376 G4cout << "Shape not one of the plane types" << G4endl;
377
378 // Apply Rotation Matrix
379 // x * Rotx, y * Roty and z * Rotz
380 if(verbosityLevel == 2)
381 {
382 G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
383 }
384 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
385 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
386 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
387
388 RandPos.setX(tempx);
389 RandPos.setY(tempy);
390 RandPos.setZ(tempz);
391
392 // Translate
393 particle_position = CentreCoords + RandPos;
394 if(verbosityLevel >= 1)
395 {
396 if(verbosityLevel == 2)
397 {
398 G4cout << "Rotated Position " << RandPos << G4endl;
399 }
400 G4cout << "Rotated and Translated position " << particle_position << G4endl;
401 }
402
403 // For Cosine-Law make SideRefVecs = to Rotation matrix vectors
404 SideRefVec1 = Rotx;
405 SideRefVec2 = Roty;
406 SideRefVec3 = Rotz;
407 // If rotation matrix z' point to origin then invert the matrix
408 // So that SideRefVecs point away.
409 if((CentreCoords.x() > 0. && Rotz.x() < 0.)
410 || (CentreCoords.x() < 0. && Rotz.x() > 0.)
411 || (CentreCoords.y() > 0. && Rotz.y() < 0.)
412 || (CentreCoords.y() < 0. && Rotz.y() > 0.)
413 || (CentreCoords.z() > 0. && Rotz.z() < 0.)
414 || (CentreCoords.z() < 0. && Rotz.z() > 0.))
415 {
416 // Invert y and z.
417 SideRefVec2 = -SideRefVec2;
418 SideRefVec3 = -SideRefVec3;
419 }
420 if(verbosityLevel == 2)
421 {
422 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
423 }
424}
425
426void G4SPSPosDistribution::GeneratePointsOnSurface()
427{
428 //Private method to create points on a surface
429 G4double theta, phi;
430 G4double x, y, z;
431 x = y = z = 0.;
432 G4ThreeVector RandPos;
433 // G4double tempx, tempy, tempz;
434
435 if(SourcePosType != "Surface" && verbosityLevel >= 1)
436 G4cout << "Error SourcePosType not Surface" << G4endl;
437
438 if(Shape == "Sphere")
439 {
440 G4double tantheta;
441 theta = posRndm->GenRandPosTheta();
442 phi = posRndm->GenRandPosPhi();
443 theta = std::acos(1. - 2.*theta); // theta isotropic
444 phi = phi * 2. * pi;
445 tantheta = std::tan(theta);
446
447 x = Radius * std::sin(theta) * std::cos(phi);
448 y = Radius * std::sin(theta) * std::sin(phi);
449 z = Radius * std::cos(theta);
450
451 RandPos.setX(x);
452 RandPos.setY(y);
453 RandPos.setZ(z);
454
455 // Cosine-law (not a good idea to use this here)
456 G4ThreeVector zdash(x,y,z);
457 zdash = zdash.unit();
458 G4ThreeVector xdash = Rotz.cross(zdash);
459 G4ThreeVector ydash = xdash.cross(zdash);
460 SideRefVec1 = xdash.unit();
461 SideRefVec2 = ydash.unit();
462 SideRefVec3 = zdash.unit();
463 }
464 else if(Shape == "Ellipsoid")
465 {
466 G4double theta, phi, minphi, maxphi, middlephi;
467 G4double answer, constant;
468
469 constant = pi/(halfx*halfx) + pi/(halfy*halfy) +
470 twopi/(halfz*halfz);
471
472 // simplified approach
473 theta = posRndm->GenRandPosTheta();
474 phi = posRndm->GenRandPosPhi();
475
476 theta = std::acos(1. - 2.*theta);
477 minphi = 0.;
478 maxphi = twopi;
479 while(maxphi-minphi > 0.)
480 {
481 middlephi = (maxphi+minphi)/2.;
482 answer = (1./(halfx*halfx))*(middlephi/2. + std::sin(2*middlephi)/4.)
483 + (1./(halfy*halfy))*(middlephi/2. - std::sin(2*middlephi)/4.)
484 + middlephi/(halfz*halfz);
485 answer = answer/constant;
486 if(answer > phi) maxphi = middlephi;
487 if(answer < phi) minphi = middlephi;
488 if(std::fabs(answer-phi) <= 0.00001)
489 {
490 minphi = maxphi +1;
491 phi = middlephi;
492 }
493 }
494
495 x = std::sin(theta)*std::cos(phi);
496 y = std::sin(theta)*std::sin(phi);
497 z = std::cos(theta);
498 // x,y and z form a unit vector. Put this onto the ellipse.
499 G4double lhs;
500 // solve for x
501 G4double numYinX = y/x;
502 G4double numZinX = z/x;
503 G4double tempxvar;
504 tempxvar= 1./(halfx*halfx)+(numYinX*numYinX)/(halfy*halfy)
505 + (numZinX*numZinX)/(halfz*halfz);
506
507 tempxvar = 1./tempxvar;
508 G4double coordx = std::sqrt(tempxvar);
509
510 //solve for y
511 G4double numXinY = x/y;
512 G4double numZinY = z/y;
513 G4double tempyvar;
514 tempyvar=(numXinY*numXinY)/(halfx*halfx)+1./(halfy*halfy)
515 +(numZinY*numZinY)/(halfz*halfz);
516 tempyvar = 1./tempyvar;
517 G4double coordy = std::sqrt(tempyvar);
518
519 //solve for z
520 G4double numXinZ = x/z;
521 G4double numYinZ = y/z;
522 G4double tempzvar;
523 tempzvar=(numXinZ*numXinZ)/(halfx*halfx)
524 +(numYinZ*numYinZ)/(halfy*halfy)+1./(halfz*halfz);
525 tempzvar = 1./tempzvar;
526 G4double coordz = std::sqrt(tempzvar);
527
528 lhs = std::sqrt((coordx*coordx)/(halfx*halfx) +
529 (coordy*coordy)/(halfy*halfy) +
530 (coordz*coordz)/(halfz*halfz));
531
532 if(std::fabs(lhs-1.) > 0.001 && verbosityLevel >= 1)
533 G4cout << "Error: theta, phi not really on ellipsoid" << G4endl;
534
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);
555
556 // Cosine-law (not a good idea to use this here)
557 G4ThreeVector zdash(coordx,coordy,coordz);
558 zdash = zdash.unit();
559 G4ThreeVector xdash = Rotz.cross(zdash);
560 G4ThreeVector ydash = xdash.cross(zdash);
561 SideRefVec1 = xdash.unit();
562 SideRefVec2 = ydash.unit();
563 SideRefVec3 = zdash.unit();
564 }
565 else if(Shape == "Cylinder")
566 {
567 G4double AreaTop, AreaBot, AreaLat;
568 G4double AreaTotal, prob1, prob2, prob3;
569 G4double testrand;
570
571 // User giver Radius and z-half length
572 // Calculate surface areas, maybe move this to
573 // a different routine.
574
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);
666
667 }
668 else if(Shape == "Para")
669 {
670 G4double testrand;
671 //Right Parallelepiped.
672 // User gives x,y,z half lengths and ParAlpha
673 // ParTheta and ParPhi
674 // +x = <1, -x >1 & <2, +y >2 & <3, -y >3 &<4
675 // +z >4 & < 5, -z >5 &<6.
676 testrand = G4UniformRand();
677 G4double AreaX = halfy * halfz * 4.;
678 G4double AreaY = halfx * halfz * 4.;
679 G4double AreaZ = halfx * halfy * 4.;
680 G4double AreaTotal = 2*(AreaX + AreaY + AreaZ);
681 G4double Probs[6];
682 Probs[0] = AreaX/AreaTotal;
683 Probs[1] = Probs[0] + AreaX/AreaTotal;
684 Probs[2] = Probs[1] + AreaY/AreaTotal;
685 Probs[3] = Probs[2] + AreaY/AreaTotal;
686 Probs[4] = Probs[3] + AreaZ/AreaTotal;
687 Probs[5] = Probs[4] + AreaZ/AreaTotal;
688
689 x = posRndm->GenRandX();
690 y = posRndm->GenRandY();
691 z = posRndm->GenRandZ();
692
693 x = x * halfx * 2.;
694 x = x - halfx;
695 y = y * halfy * 2.;
696 y = y - halfy;
697 z = z * halfz * 2.;
698 z = z - halfz;
699 // Pick a side first
700 if(testrand < Probs[0])
701 {
702 // side is +x
703 x = halfx + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
704 y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
705 z = z;
706 // Cosine-law
707 G4ThreeVector xdash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
708 halfz*std::tan(ParTheta)*std::sin(ParPhi),
709 halfz/std::cos(ParPhi));
710 G4ThreeVector ydash(halfy*std::tan(ParAlpha), -halfy, 0.0);
711 xdash = xdash.unit();
712 ydash = ydash.unit();
713 G4ThreeVector zdash = xdash.cross(ydash);
714 SideRefVec1 = xdash.unit();
715 SideRefVec2 = ydash.unit();
716 SideRefVec3 = zdash.unit();
717 }
718 else if(testrand >= Probs[0] && testrand < Probs[1])
719 {
720 // side is -x
721 x = -halfx + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
722 y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
723 z = z;
724 // Cosine-law
725 G4ThreeVector xdash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
726 halfz*std::tan(ParTheta)*std::sin(ParPhi),
727 halfz/std::cos(ParPhi));
728 G4ThreeVector ydash(halfy*std::tan(ParAlpha), halfy, 0.0);
729 xdash = xdash.unit();
730 ydash = ydash.unit();
731 G4ThreeVector zdash = xdash.cross(ydash);
732 SideRefVec1 = xdash.unit();
733 SideRefVec2 = ydash.unit();
734 SideRefVec3 = zdash.unit();
735 }
736 else if(testrand >= Probs[1] && testrand < Probs[2])
737 {
738 // side is +y
739 x = x + z*std::tan(ParTheta)*std::cos(ParPhi) + halfy*std::tan(ParAlpha);
740 y = halfy + z*std::tan(ParTheta)*std::sin(ParPhi);
741 z = z;
742 // Cosine-law
743 G4ThreeVector ydash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
744 halfz*std::tan(ParTheta)*std::sin(ParPhi),
745 halfz/std::cos(ParPhi));
746 ydash = ydash.unit();
747 G4ThreeVector xdash = Roty.cross(ydash);
748 G4ThreeVector zdash = xdash.cross(ydash);
749 SideRefVec1 = xdash.unit();
750 SideRefVec2 = -ydash.unit();
751 SideRefVec3 = -zdash.unit();
752 }
753 else if(testrand >= Probs[2] && testrand < Probs[3])
754 {
755 // side is -y
756 x = x + z*std::tan(ParTheta)*std::cos(ParPhi) - halfy*std::tan(ParAlpha);
757 y = -halfy + z*std::tan(ParTheta)*std::sin(ParPhi);
758 z = z;
759 // Cosine-law
760 G4ThreeVector ydash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
761 halfz*std::tan(ParTheta)*std::sin(ParPhi),
762 halfz/std::cos(ParPhi));
763 ydash = ydash.unit();
764 G4ThreeVector xdash = Roty.cross(ydash);
765 G4ThreeVector zdash = xdash.cross(ydash);
766 SideRefVec1 = xdash.unit();
767 SideRefVec2 = ydash.unit();
768 SideRefVec3 = zdash.unit();
769 }
770 else if(testrand >= Probs[3] && testrand < Probs[4])
771 {
772 // side is +z
773 z = halfz;
774 y = y + halfz*std::sin(ParPhi)*std::tan(ParTheta);
775 x = x + halfz*std::cos(ParPhi)*std::tan(ParTheta) + y*std::tan(ParAlpha);
776 // Cosine-law
777 SideRefVec1 = Rotx;
778 SideRefVec2 = Roty;
779 SideRefVec3 = Rotz;
780 }
781 else if(testrand >= Probs[4] && testrand < Probs[5])
782 {
783 // side is -z
784 z = -halfz;
785 y = y - halfz*std::sin(ParPhi)*std::tan(ParTheta);
786 x = x - halfz*std::cos(ParPhi)*std::tan(ParTheta) + y*std::tan(ParAlpha);
787 // Cosine-law
788 SideRefVec1 = Rotx;
789 SideRefVec2 = -Roty;
790 SideRefVec3 = -Rotz;
791 }
792 else
793 {
794 G4cout << "Error: testrand out of range" << G4endl;
795 if(verbosityLevel >= 1)
796 G4cout << "testrand=" << testrand << " Probs[5]=" << Probs[5] <<G4endl;
797 }
798
799 RandPos.setX(x);
800 RandPos.setY(y);
801 RandPos.setZ(z);
802 }
803
804 // Apply Rotation Matrix
805 // x * Rotx, y * Roty and z * Rotz
806 if(verbosityLevel == 2)
807 G4cout << "Raw position " << RandPos << G4endl;
808
809 x=(RandPos.x()*Rotx.x())+(RandPos.y()*Roty.x())+(RandPos.z()*Rotz.x());
810 y=(RandPos.x()*Rotx.y())+(RandPos.y()*Roty.y())+(RandPos.z()*Rotz.y());
811 z=(RandPos.x()*Rotx.z())+(RandPos.y()*Roty.z())+(RandPos.z()*Rotz.z());
812
813 RandPos.setX(x);
814 RandPos.setY(y);
815 RandPos.setZ(z);
816
817 // Translate
818 particle_position = CentreCoords + RandPos;
819
820 if(verbosityLevel >= 1)
821 {
822 if(verbosityLevel == 2)
823 G4cout << "Rotated position " << RandPos << G4endl;
824 G4cout << "Rotated and translated position " << particle_position << G4endl;
825 }
826 if(verbosityLevel == 2)
827 {
828 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
829 }
830}
831
832void G4SPSPosDistribution::GeneratePointsInVolume()
833{
834 G4ThreeVector RandPos;
835 G4double tempx, tempy, tempz;
836 G4double x, y, z;
837 x = y = z = 0.;
838 if(SourcePosType != "Volume" && verbosityLevel >= 1)
839 G4cout << "Error SourcePosType not Volume" << G4endl;
840 //Private method to create points in a volume
841 if(Shape == "Sphere")
842 {
843 x = Radius*2.;
844 y = Radius*2.;
845 z = Radius*2.;
846 while(((x*x)+(y*y)+(z*z)) > (Radius*Radius))
847 {
848 x = posRndm->GenRandX();
849 y = posRndm->GenRandY();
850 z = posRndm->GenRandZ();
851
852 x = (x*2.*Radius) - Radius;
853 y = (y*2.*Radius) - Radius;
854 z = (z*2.*Radius) - Radius;
855 }
856 }
857 else if(Shape == "Ellipsoid")
858 {
859 G4double temp;
860 temp = 100.;
861 while(temp > 1.)
862 {
863 x = posRndm->GenRandX();
864 y = posRndm->GenRandY();
865 z = posRndm->GenRandZ();
866
867 x = (x*2.*halfx) - halfx;
868 y = (y*2.*halfy) - halfy;
869 z = (z*2.*halfz) - halfz;
870
871 temp = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy))
872 + ((z*z)/(halfz*halfz));
873 }
874 }
875 else if(Shape == "Cylinder")
876 {
877 x = Radius*2.;
878 y = Radius*2.;
879 while(((x*x)+(y*y)) > (Radius*Radius))
880 {
881 x = posRndm->GenRandX();
882 y = posRndm->GenRandY();
883 z = posRndm->GenRandZ();
884
885 x = (x*2.*Radius) - Radius;
886 y = (y*2.*Radius) - Radius;
887 z = (z*2.*halfz) - halfz;
888 }
889 }
890 else if(Shape == "Para")
891 {
892 x = posRndm->GenRandX();
893 y = posRndm->GenRandY();
894 z = posRndm->GenRandZ();
895 x = (x*2.*halfx) - halfx;
896 y = (y*2.*halfy) - halfy;
897 z = (z*2.*halfz) - halfz;
898 x = x + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
899 y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
900 z = z;
901 }
902 else
903 G4cout << "Error: Volume Shape Doesnt Exist" << G4endl;
904
905 RandPos.setX(x);
906 RandPos.setY(y);
907 RandPos.setZ(z);
908
909 // Apply Rotation Matrix
910 // x * Rotx, y * Roty and z * Rotz
911 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
912 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
913 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
914
915 RandPos.setX(tempx);
916 RandPos.setY(tempy);
917 RandPos.setZ(tempz);
918
919 // Translate
920 particle_position = CentreCoords + RandPos;
921
922 if(verbosityLevel == 2)
923 {
924 G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
925 G4cout << "Rotated position " << RandPos << G4endl;
926 }
927 if(verbosityLevel >= 1)
928 G4cout << "Rotated and translated position " << particle_position << G4endl;
929
930 // Cosine-law (not a good idea to use this here)
931 G4ThreeVector zdash(tempx,tempy,tempz);
932 zdash = zdash.unit();
933 G4ThreeVector xdash = Rotz.cross(zdash);
934 G4ThreeVector ydash = xdash.cross(zdash);
935 SideRefVec1 = xdash.unit();
936 SideRefVec2 = ydash.unit();
937 SideRefVec3 = zdash.unit();
938
939 if(verbosityLevel == 2)
940 {
941 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
942 }
943}
944
945G4bool G4SPSPosDistribution::IsSourceConfined()
946{
947 // Method to check point is within the volume specified
948 if(Confine == false)
949 G4cout << "Error: Confine is false" << G4endl;
950 G4ThreeVector null(0.,0.,0.);
951 G4ThreeVector *ptr;
952 ptr = &null;
953
954 // Check particle_position is within VolName, if so true,
955 // else false
956 G4VPhysicalVolume *theVolume;
957 theVolume=gNavigator->LocateGlobalPointAndSetup(particle_position,ptr,true);
958 G4String theVolName = theVolume->GetName();
959 if(theVolName == VolName)
960 {
961 if(verbosityLevel >= 1)
962 G4cout << "Particle is in volume " << VolName << G4endl;
963 return(true);
964 }
965 else
966 return(false);
967}
968
969G4ThreeVector G4SPSPosDistribution::GenerateOne()
970{
971 //
972 G4bool srcconf = false;
973 G4int LoopCount = 0;
974 while(srcconf == false)
975 {
976 if(SourcePosType == "Point")
977 GeneratePointSource();
978 else if(SourcePosType == "Beam")
979 GeneratePointsInBeam();
980 else if(SourcePosType == "Plane")
981 GeneratePointsInPlane();
982 else if(SourcePosType == "Surface")
983 GeneratePointsOnSurface();
984 else if(SourcePosType == "Volume")
985 GeneratePointsInVolume();
986 else
987 {
988 G4cout << "Error: SourcePosType undefined" << G4endl;
989 G4cout << "Generating point source" << G4endl;
990 GeneratePointSource();
991 }
992 if(Confine == true)
993 {
994 srcconf = IsSourceConfined();
995 // if source in confined srcconf = true terminating the loop
996 // if source isnt confined srcconf = false and loop continues
997 }
998 else if(Confine == false)
999 srcconf = true; // terminate loop
1000 LoopCount++;
1001 if(LoopCount == 100000)
1002 {
1003 G4cout << "*************************************" << G4endl;
1004 G4cout << "LoopCount = 100000" << G4endl;
1005 G4cout << "Either the source distribution >> confinement" << G4endl;
1006 G4cout << "or any confining volume may not overlap with" << G4endl;
1007 G4cout << "the source distribution or any confining volumes" << G4endl;
1008 G4cout << "may not exist"<< G4endl;
1009 G4cout << "If you have set confine then this will be ignored" <<G4endl;
1010 G4cout << "for this event." << G4endl;
1011 G4cout << "*************************************" << G4endl;
1012 srcconf = true; //Avoids an infinite loop
1013 }
1014 }
1015 return particle_position;
1016}
1017
1018
1019
1020
Note: See TracBrowser for help on using the repository browser.