| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 G4RandomRadi
39 // G4RandomPointOnEll
40 // 07.11.08 - P.Gumplinger, based on implement
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 (re
57 //
58 inline G4ThreeVector G4LambertianRand(const G4
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
78
79 return vect;
80 }
81
82 // -------------------------------------------
83 // Chooses a random vector within a plane give
84 //
85 inline G4ThreeVector G4PlaneVectorRand(const G
86 {
87 G4ThreeVector vec1 = normal.orthogonal();
88 G4ThreeVector vec2 = vec1.cross(normal);
89
90 G4double phi = CLHEP::twopi * G4UniformRa
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
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 * rm
109 }
110
111 // -------------------------------------------
112 // Returns a random point in ellipse (x/a)^2 +
113 // (rejection sampling)
114 //
115 inline G4TwoVector G4RandomPointInEllipse(G4do
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 +
131 // (rejection sampling)
132 //
133 inline G4TwoVector G4RandomPointOnEllipse(G4do
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 * G4UniformRan
143 x = std::cos(phi);
144 y = std::sin(phi);
145 G4double mu = std::sqrt((B * x) * (B * x)
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
154 // (rejection sampling)
155 //
156 inline G4ThreeVector G4RandomPointOnEllipsoid(
157
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
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 *
172 if(mu_max * G4UniformRand() <= mu)
173 break;
174 }
175 return G4ThreeVector(A * p.x(), B * p.y(), C
176 }
177
178 #endif /* G4RANDOMTOOLS_HH */
179