| Geant4 Cross Reference |
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 //
29 // ---------------------------------------------------------------------------
30 // GEANT 4 class header file
31 // ---------------------------------------------------------------------------
32 // Class description:
33 //
34 // Utility functions
35
36 // History:
37 //
38 // 24.08.17 - E.Tcherniaev, added G4RandomRadiusInRing, G4RandomPointInEllipse
39 // G4RandomPointOnEllipse, G4RandomPointOnEllipsoid
40 // 07.11.08 - P.Gumplinger, based on implementation in G4OpBoundaryProcess
41 //
42 // ---------------------------------------------------------------------------
43
44 #ifndef G4RANDOMTOOLS_HH
45 #define G4RANDOMTOOLS_HH
46
47 #include <CLHEP/Units/PhysicalConstants.h>
48
49 #include "G4RandomDirection.hh"
50 #include "G4ThreeVector.hh"
51 #include "G4TwoVector.hh"
52 #include "Randomize.hh"
53 #include "globals.hh"
54
55 // ---------------------------------------------------------------------------
56 // Returns a random lambertian unit vector (rejection sampling)
57 //
58 inline G4ThreeVector G4LambertianRand(const G4ThreeVector& normal)
59 {
60 G4ThreeVector vect;
61 G4double ndotv;
62 G4int count = 0;
63 const G4int max_trials = 1024;
64
65 do
66 {
67 ++count;
68 vect = G4RandomDirection();
69 ndotv = normal * vect;
70
71 if(ndotv < 0.0)
72 {
73 vect = -vect;
74 ndotv = -ndotv;
75 }
76
77 } while(!(G4UniformRand() < ndotv) && (count < max_trials));
78
79 return vect;
80 }
81
82 // ---------------------------------------------------------------------------
83 // Chooses a random vector within a plane given by the unit normal
84 //
85 inline G4ThreeVector G4PlaneVectorRand(const G4ThreeVector& normal)
86 {
87 G4ThreeVector vec1 = normal.orthogonal();
88 G4ThreeVector vec2 = vec1.cross(normal);
89
90 G4double phi = CLHEP::twopi * G4UniformRand();
91 G4double cosphi = std::cos(phi);
92 G4double sinphi = std::sin(phi);
93
94 return cosphi * vec1 + sinphi * vec2;
95 }
96
97 // ---------------------------------------------------------------------------
98 // Returns a random radius in annular ring
99 //
100 inline G4double G4RandomRadiusInRing(G4double rmin, G4double rmax)
101 {
102 if(rmin == rmax)
103 {
104 return rmin;
105 }
106 G4double k = G4UniformRand();
107 return (rmin <= 0) ? rmax * std::sqrt(k)
108 : std::sqrt(k * rmax * rmax + (1. - k) * rmin * rmin);
109 }
110
111 // ---------------------------------------------------------------------------
112 // Returns a random point in ellipse (x/a)^2 + (y/b)^2 = 1
113 // (rejection sampling)
114 //
115 inline G4TwoVector G4RandomPointInEllipse(G4double a, G4double b)
116 {
117 G4double aa = (a * a == 0) ? 0 : 1 / (a * a);
118 G4double bb = (b * b == 0) ? 0 : 1 / (b * b);
119 for(G4int i = 0; i < 1000; ++i)
120 {
121 G4double x = a * (2 * G4UniformRand() - 1);
122 G4double y = b * (2 * G4UniformRand() - 1);
123 if(x * x * aa + y * y * bb <= 1)
124 return G4TwoVector(x, y);
125 }
126 return G4TwoVector(0, 0);
127 }
128
129 // ---------------------------------------------------------------------------
130 // Returns a random point on ellipse (x/a)^2 + (y/b)^2 = 1
131 // (rejection sampling)
132 //
133 inline G4TwoVector G4RandomPointOnEllipse(G4double a, G4double b)
134 {
135 G4double A = std::abs(a);
136 G4double B = std::abs(b);
137 G4double mu_max = std::max(A, B);
138
139 G4double x, y;
140 for(G4int i = 0; i < 1000; ++i)
141 {
142 G4double phi = CLHEP::twopi * G4UniformRand();
143 x = std::cos(phi);
144 y = std::sin(phi);
145 G4double mu = std::sqrt((B * x) * (B * x) + (A * y) * (A * y));
146 if(mu_max * G4UniformRand() <= mu)
147 break;
148 }
149 return G4TwoVector(A * x, B * y);
150 }
151
152 // ---------------------------------------------------------------------------
153 // Returns a random point on ellipsoid (x/a)^2 + (y/b)^2 + (z/c)^2 = 1
154 // (rejection sampling)
155 //
156 inline G4ThreeVector G4RandomPointOnEllipsoid(G4double a, G4double b,
157 G4double c)
158 {
159 G4double A = std::abs(a);
160 G4double B = std::abs(b);
161 G4double C = std::abs(c);
162 G4double mu_max = std::max(std::max(A * B, A * C), B * C);
163
164 G4ThreeVector p;
165 for(G4int i = 0; i < 1000; ++i)
166 {
167 p = G4RandomDirection();
168 G4double xbc = p.x() * B * C;
169 G4double yac = p.y() * A * C;
170 G4double zab = p.z() * A * B;
171 G4double mu = std::sqrt(xbc * xbc + yac * yac + zab * zab);
172 if(mu_max * G4UniformRand() <= mu)
173 break;
174 }
175 return G4ThreeVector(A * p.x(), B * p.y(), C * p.z());
176 }
177
178 #endif /* G4RANDOMTOOLS_HH */
179