include/Geant4/G4RandomTools.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 //
0029 // ---------------------------------------------------------------------------
0030 //      GEANT 4 class header file
0031 // ---------------------------------------------------------------------------
0032 // Class description:
0033 //
0034 // Utility functions
0035
0036 // History:
0037 //
0038 // 24.08.17 - E.Tcherniaev, added G4RandomRadiusInRing, G4RandomPointInEllipse
0039 //                          G4RandomPointOnEllipse, G4RandomPointOnEllipsoid
0040 // 07.11.08 - P.Gumplinger, based on implementation in G4OpBoundaryProcess
0041 //
0042 // ---------------------------------------------------------------------------
0043
0044 #ifndef G4RANDOMTOOLS_HH
0045 #define G4RANDOMTOOLS_HH
0046
0047 #include <CLHEP/Units/PhysicalConstants.h>
0048
0049 #include "G4RandomDirection.hh"
0050 #include "G4ThreeVector.hh"
0051 #include "G4TwoVector.hh"
0052 #include "Randomize.hh"
0053 #include "globals.hh"
0054
0055 // ---------------------------------------------------------------------------
0056 // Returns a random lambertian unit vector (rejection sampling)
0057 //
0058 inline G4ThreeVector G4LambertianRand(const G4ThreeVector& normal)
0059 {
0060   G4ThreeVector vect;
0061   G4double ndotv;
0062   G4int count            = 0;
0063   const G4int max_trials = 1024;
0064
0065   do
0066   {
0067     ++count;
0068     vect  = G4RandomDirection();
0069     ndotv = normal * vect;
0070
0071     if(ndotv < 0.0)
0072     {
0073       vect  = -vect;
0074       ndotv = -ndotv;
0075     }
0076
0077   } while(!(G4UniformRand() < ndotv) && (count < max_trials));
0078
0079   return vect;
0080 }
0081
0082 // ---------------------------------------------------------------------------
0083 // Chooses a random vector within a plane given by the unit normal
0084 //
0085 inline G4ThreeVector G4PlaneVectorRand(const G4ThreeVector& normal)
0086 {
0087   G4ThreeVector vec1 = normal.orthogonal();
0088   G4ThreeVector vec2 = vec1.cross(normal);
0089
0090   G4double phi    = CLHEP::twopi * G4UniformRand();
0091   G4double cosphi = std::cos(phi);
0092   G4double sinphi = std::sin(phi);
0093
0094   return cosphi * vec1 + sinphi * vec2;
0095 }
0096
0097 // ---------------------------------------------------------------------------
0098 // Returns a random radius in annular ring
0099 //
0100 inline G4double G4RandomRadiusInRing(G4double rmin, G4double rmax)
0101 {
0102   if(rmin == rmax)
0103   {
0104     return rmin;
0105   }
0106   G4double k = G4UniformRand();
0107   return (rmin <= 0) ? rmax * std::sqrt(k)
0108                      : std::sqrt(k * rmax * rmax + (1. - k) * rmin * rmin);
0109 }
0110
0111 // ---------------------------------------------------------------------------
0112 // Returns a random point in ellipse (x/a)^2 + (y/b)^2 = 1
0113 // (rejection sampling)
0114 //
0115 inline G4TwoVector G4RandomPointInEllipse(G4double a, G4double b)
0116 {
0117   G4double aa = (a * a == 0) ? 0 : 1 / (a * a);
0118   G4double bb = (b * b == 0) ? 0 : 1 / (b * b);
0119   for(G4int i = 0; i < 1000; ++i)
0120   {
0121     G4double x = a * (2 * G4UniformRand() - 1);
0122     G4double y = b * (2 * G4UniformRand() - 1);
0123     if(x * x * aa + y * y * bb <= 1)
0124       return G4TwoVector(x, y);
0125   }
0126   return G4TwoVector(0, 0);
0127 }
0128
0129 // ---------------------------------------------------------------------------
0130 // Returns a random point on ellipse (x/a)^2 + (y/b)^2 = 1
0131 // (rejection sampling)
0132 //
0133 inline G4TwoVector G4RandomPointOnEllipse(G4double a, G4double b)
0134 {
0135   G4double A      = std::abs(a);
0136   G4double B      = std::abs(b);
0137   G4double mu_max = std::max(A, B);
0138
0139   G4double x, y;
0140   for(G4int i = 0; i < 1000; ++i)
0141   {
0142     G4double phi = CLHEP::twopi * G4UniformRand();
0143     x            = std::cos(phi);
0144     y            = std::sin(phi);
0145     G4double mu  = std::sqrt((B * x) * (B * x) + (A * y) * (A * y));
0146     if(mu_max * G4UniformRand() <= mu)
0147       break;
0148   }
0149   return G4TwoVector(A * x, B * y);
0150 }
0151
0152 // ---------------------------------------------------------------------------
0153 // Returns a random point on ellipsoid (x/a)^2 + (y/b)^2 + (z/c)^2 = 1
0154 // (rejection sampling)
0155 //
0156 inline G4ThreeVector G4RandomPointOnEllipsoid(G4double a, G4double b,
0157                                               G4double c)
0158 {
0159   G4double A      = std::abs(a);
0160   G4double B      = std::abs(b);
0161   G4double C      = std::abs(c);
0162   G4double mu_max = std::max(std::max(A * B, A * C), B * C);
0163
0164   G4ThreeVector p;
0165   for(G4int i = 0; i < 1000; ++i)
0166   {
0167     p            = G4RandomDirection();
0168     G4double xbc = p.x() * B * C;
0169     G4double yac = p.y() * A * C;
0170     G4double zab = p.z() * A * B;
0171     G4double mu  = std::sqrt(xbc * xbc + yac * yac + zab * zab);
0172     if(mu_max * G4UniformRand() <= mu)
0173       break;
0174   }
0175   return G4ThreeVector(A * p.x(), B * p.y(), C * p.z());
0176 }
0177
0178 #endif /* G4RANDOMTOOLS_HH */