ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RandomTools.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RandomTools.hh
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 
48 
49 #include "globals.hh"
50 #include "Randomize.hh"
51 #include "G4TwoVector.hh"
52 #include "G4ThreeVector.hh"
53 #include "G4RandomDirection.hh"
54 
55 // ---------------------------------------------------------------------------
56 // Returns a random lambertian unit vector (rejection sampling)
57 //
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 //
86 {
87  G4ThreeVector vec1 = normal.orthogonal();
88  G4ThreeVector vec2 = vec1.cross(normal);
89 
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 //
101 {
102  if (rmin == rmax)
103  {
104  return rmin;
105  }
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 //
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) return G4TwoVector(x,y);
124  }
125  return G4TwoVector(0,0);
126 }
127 
128 // ---------------------------------------------------------------------------
129 // Returns a random point on ellipse (x/a)^2 + (y/b)^2 = 1
130 // (rejection sampling)
131 //
133 {
134  G4double A = std::abs(a);
135  G4double B = std::abs(b);
136  G4double mu_max = std::max(A,B);
137 
138  G4double x,y;
139  for (G4int i=0; i<1000; ++i)
140  {
142  x = std::cos(phi);
143  y = std::sin(phi);
144  G4double mu = std::sqrt((B*x)*(B*x) + (A*y)*(A*y));
145  if (mu_max*G4UniformRand() <= mu) break;
146  }
147  return G4TwoVector(A*x,B*y);
148 }
149 
150 // ---------------------------------------------------------------------------
151 // Returns a random point on ellipsoid (x/a)^2 + (y/b)^2 + (z/c)^2 = 1
152 // (rejection sampling)
153 //
154 inline
156 {
157  G4double A = std::abs(a);
158  G4double B = std::abs(b);
159  G4double C = std::abs(c);
160  G4double mu_max = std::max(std::max(A*B,A*C),B*C);
161 
163  for (G4int i=0; i<1000; ++i)
164  {
165  p = G4RandomDirection();
166  G4double xbc = p.x()*B*C;
167  G4double yac = p.y()*A*C;
168  G4double zab = p.z()*A*B;
169  G4double mu = std::sqrt(xbc*xbc + yac*yac + zab*zab);
170  if (mu_max*G4UniformRand() <= mu) break;
171  }
172  return G4ThreeVector(A*p.x(),B*p.y(),C*p.z());
173 }
174 
175 #endif /* G4RANDOMTOOLS_HH */