Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // >> 27 // $Id: G4RandomTools.hh 105368 2017-07-24 09:44:20Z gcosmo $ 27 // 28 // 28 // << 29 // 29 // ------------------------------------------- 30 // --------------------------------------------------------------------------- 30 // GEANT 4 class header file << 31 // GEANT 4 class header file 31 // ------------------------------------------- 32 // --------------------------------------------------------------------------- 32 // Class description: 33 // Class description: 33 // 34 // 34 // Utility functions << 35 // Utility functions 35 36 36 // History: 37 // History: 37 // 38 // 38 // 24.08.17 - E.Tcherniaev, added G4RandomRadi 39 // 24.08.17 - E.Tcherniaev, added G4RandomRadiusInRing, G4RandomPointInEllipse 39 // G4RandomPointOnEll 40 // G4RandomPointOnEllipse, G4RandomPointOnEllipsoid 40 // 07.11.08 - P.Gumplinger, based on implement 41 // 07.11.08 - P.Gumplinger, based on implementation in G4OpBoundaryProcess 41 // 42 // 42 // ------------------------------------------- 43 // --------------------------------------------------------------------------- 43 44 44 #ifndef G4RANDOMTOOLS_HH 45 #ifndef G4RANDOMTOOLS_HH 45 #define G4RANDOMTOOLS_HH 46 #define G4RANDOMTOOLS_HH 46 47 47 #include <CLHEP/Units/PhysicalConstants.h> 48 #include <CLHEP/Units/PhysicalConstants.h> 48 49 49 #include "G4RandomDirection.hh" << 50 #include "G4ThreeVector.hh" << 51 #include "G4TwoVector.hh" << 52 #include "Randomize.hh" << 53 #include "globals.hh" 50 #include "globals.hh" >> 51 #include "Randomize.hh" >> 52 #include "G4TwoVector.hh" >> 53 #include "G4ThreeVector.hh" >> 54 #include "G4RandomDirection.hh" 54 55 55 // ------------------------------------------- 56 // --------------------------------------------------------------------------- 56 // Returns a random lambertian unit vector (re 57 // Returns a random lambertian unit vector (rejection sampling) 57 // 58 // 58 inline G4ThreeVector G4LambertianRand(const G4 59 inline G4ThreeVector G4LambertianRand(const G4ThreeVector& normal) 59 { 60 { 60 G4ThreeVector vect; 61 G4ThreeVector vect; 61 G4double ndotv; 62 G4double ndotv; 62 G4int count = 0; << 63 G4int count=0; 63 const G4int max_trials = 1024; 64 const G4int max_trials = 1024; 64 65 65 do 66 do 66 { 67 { 67 ++count; 68 ++count; 68 vect = G4RandomDirection(); << 69 vect = G4RandomDirection(); 69 ndotv = normal * vect; 70 ndotv = normal * vect; 70 71 71 if(ndotv < 0.0) << 72 if (ndotv < 0.0) 72 { 73 { 73 vect = -vect; << 74 vect = -vect; 74 ndotv = -ndotv; 75 ndotv = -ndotv; 75 } 76 } 76 77 77 } while(!(G4UniformRand() < ndotv) && (count << 78 } while (!(G4UniformRand() < ndotv) && (count < max_trials)); 78 79 79 return vect; 80 return vect; 80 } 81 } 81 82 82 // ------------------------------------------- 83 // --------------------------------------------------------------------------- 83 // Chooses a random vector within a plane give 84 // Chooses a random vector within a plane given by the unit normal 84 // 85 // 85 inline G4ThreeVector G4PlaneVectorRand(const G 86 inline G4ThreeVector G4PlaneVectorRand(const G4ThreeVector& normal) 86 { 87 { 87 G4ThreeVector vec1 = normal.orthogonal(); 88 G4ThreeVector vec1 = normal.orthogonal(); 88 G4ThreeVector vec2 = vec1.cross(normal); 89 G4ThreeVector vec2 = vec1.cross(normal); 89 90 90 G4double phi = CLHEP::twopi * G4UniformRa << 91 G4double phi = CLHEP::twopi*G4UniformRand(); 91 G4double cosphi = std::cos(phi); 92 G4double cosphi = std::cos(phi); 92 G4double sinphi = std::sin(phi); 93 G4double sinphi = std::sin(phi); 93 94 94 return cosphi * vec1 + sinphi * vec2; 95 return cosphi * vec1 + sinphi * vec2; 95 } 96 } 96 97 97 // ------------------------------------------- 98 // --------------------------------------------------------------------------- 98 // Returns a random radius in annular ring 99 // Returns a random radius in annular ring 99 // 100 // 100 inline G4double G4RandomRadiusInRing(G4double 101 inline G4double G4RandomRadiusInRing(G4double rmin, G4double rmax) 101 { 102 { 102 if(rmin == rmax) << 103 if (rmin == rmax) 103 { 104 { 104 return rmin; 105 return rmin; 105 } 106 } 106 G4double k = G4UniformRand(); 107 G4double k = G4UniformRand(); 107 return (rmin <= 0) ? rmax * std::sqrt(k) << 108 return (rmin <= 0) ? rmax*std::sqrt(k) 108 : std::sqrt(k * rmax * rm << 109 : std::sqrt(k*rmax*rmax + (1.-k)*rmin*rmin); 109 } 110 } 110 111 111 // ------------------------------------------- 112 // --------------------------------------------------------------------------- 112 // Returns a random point in ellipse (x/a)^2 + 113 // Returns a random point in ellipse (x/a)^2 + (y/b)^2 = 1 113 // (rejection sampling) 114 // (rejection sampling) 114 // 115 // 115 inline G4TwoVector G4RandomPointInEllipse(G4do 116 inline G4TwoVector G4RandomPointInEllipse(G4double a, G4double b) 116 { 117 { 117 G4double aa = (a * a == 0) ? 0 : 1 / (a * a) << 118 G4double aa = (a*a == 0) ? 0 : 1/(a*a); 118 G4double bb = (b * b == 0) ? 0 : 1 / (b * b) << 119 G4double bb = (b*b == 0) ? 0 : 1/(b*b); 119 for(G4int i = 0; i < 1000; ++i) << 120 for (G4int i=0; i<1000; ++i) 120 { 121 { 121 G4double x = a * (2 * G4UniformRand() - 1) << 122 G4double x = a*(2*G4UniformRand() - 1); 122 G4double y = b * (2 * G4UniformRand() - 1) << 123 G4double y = b*(2*G4UniformRand() - 1); 123 if(x * x * aa + y * y * bb <= 1) << 124 if (x*x*aa + y*y*bb <= 1) return G4TwoVector(x,y); 124 return G4TwoVector(x, y); << 125 } 125 } 126 return G4TwoVector(0, 0); << 126 return G4TwoVector(0,0); 127 } 127 } 128 128 129 // ------------------------------------------- 129 // --------------------------------------------------------------------------- 130 // Returns a random point on ellipse (x/a)^2 + 130 // Returns a random point on ellipse (x/a)^2 + (y/b)^2 = 1 131 // (rejection sampling) 131 // (rejection sampling) 132 // 132 // 133 inline G4TwoVector G4RandomPointOnEllipse(G4do 133 inline G4TwoVector G4RandomPointOnEllipse(G4double a, G4double b) 134 { 134 { 135 G4double A = std::abs(a); << 135 G4double A = std::abs(a); 136 G4double B = std::abs(b); << 136 G4double B = std::abs(b); 137 G4double mu_max = std::max(A, B); << 137 G4double mu_max = std::max(A,B); 138 138 139 G4double x, y; << 139 G4double x,y; 140 for(G4int i = 0; i < 1000; ++i) << 140 for (G4int i=0; i<1000; ++i) 141 { 141 { 142 G4double phi = CLHEP::twopi * G4UniformRan << 142 G4double phi = CLHEP::twopi*G4UniformRand(); 143 x = std::cos(phi); << 143 x = std::cos(phi); 144 y = std::sin(phi); << 144 y = std::sin(phi); 145 G4double mu = std::sqrt((B * x) * (B * x) << 145 G4double mu = std::sqrt((B*x)*(B*x) + (A*y)*(A*y)); 146 if(mu_max * G4UniformRand() <= mu) << 146 if (mu_max*G4UniformRand() <= mu) break; 147 break; << 148 } 147 } 149 return G4TwoVector(A * x, B * y); << 148 return G4TwoVector(A*x,B*y); 150 } 149 } 151 150 152 // ------------------------------------------- 151 // --------------------------------------------------------------------------- 153 // Returns a random point on ellipsoid (x/a)^2 152 // Returns a random point on ellipsoid (x/a)^2 + (y/b)^2 + (z/c)^2 = 1 154 // (rejection sampling) 153 // (rejection sampling) 155 // 154 // 156 inline G4ThreeVector G4RandomPointOnEllipsoid( << 155 inline 157 << 156 G4ThreeVector G4RandomPointOnEllipsoid(G4double a, G4double b, G4double c) 158 { 157 { 159 G4double A = std::abs(a); << 158 G4double A = std::abs(a); 160 G4double B = std::abs(b); << 159 G4double B = std::abs(b); 161 G4double C = std::abs(c); << 160 G4double C = std::abs(c); 162 G4double mu_max = std::max(std::max(A * B, A << 161 G4double mu_max = std::max(std::max(A*B,A*C),B*C); 163 162 164 G4ThreeVector p; 163 G4ThreeVector p; 165 for(G4int i = 0; i < 1000; ++i) << 164 for (G4int i=0; i<1000; ++i) 166 { 165 { 167 p = G4RandomDirection(); << 166 p = G4RandomDirection(); 168 G4double xbc = p.x() * B * C; << 167 G4double xbc = p.x()*B*C; 169 G4double yac = p.y() * A * C; << 168 G4double yac = p.y()*A*C; 170 G4double zab = p.z() * A * B; << 169 G4double zab = p.z()*A*B; 171 G4double mu = std::sqrt(xbc * xbc + yac * << 170 G4double mu = std::sqrt(xbc*xbc + yac*yac + zab*zab); 172 if(mu_max * G4UniformRand() <= mu) << 171 if (mu_max*G4UniformRand() <= mu) break; 173 break; << 174 } 172 } 175 return G4ThreeVector(A * p.x(), B * p.y(), C << 173 return G4ThreeVector(A*p.x(),B*p.y(),C*p.z()); 176 } 174 } 177 175 178 #endif /* G4RANDOMTOOLS_HH */ << 176 #endif /* G4RANDOMTOOLS_HH */ 179 177