| 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 // ------------------------------------------- 27 // --------------------------------------------------------------------
28 // GEANT 4 class header file 28 // GEANT 4 class header file
29 // 29 //
30 // 30 //
31 // FastAerosol 31 // FastAerosol
32 // 32 //
33 // Class description: 33 // Class description:
34 // 34 //
35 // A FastAerosol is a collection of points i 35 // A FastAerosol is a collection of points in a voxelized
36 // arbitrarily-shaped volume with methods im 36 // arbitrarily-shaped volume with methods implementing population of
37 // grids/voxels and for efficiently finding 37 // grids/voxels and for efficiently finding the nearest point.
38 // 38 //
39 // Author: A.Knaian (ara@nklabs.com), N.MacFad 39 // Author: A.Knaian (ara@nklabs.com), N.MacFadden (natemacfadden@gmail.com)
40 // ------------------------------------------- 40 // --------------------------------------------------------------------
41 41
42 #ifndef FastAerosol_h 42 #ifndef FastAerosol_h
43 #define FastAerosol_h 43 #define FastAerosol_h
44 44
45 #include "globals.hh" 45 #include "globals.hh"
46 #include "Randomize.hh" 46 #include "Randomize.hh"
47 #include "G4ThreeVector.hh" 47 #include "G4ThreeVector.hh"
48 #include "G4VSolid.hh" 48 #include "G4VSolid.hh"
49 49
50 // rotations and number density distribution 50 // rotations and number density distribution
51 #include <functional> 51 #include <functional>
52 #include "G4RotationMatrix.hh" 52 #include "G4RotationMatrix.hh"
53 #include <atomic> 53 #include <atomic>
54 54
55 //using namespace std; 55 //using namespace std;
56 56
57 class FastAerosol << 57 class FastAerosol {
58 { << 58 public:
59 public: << 59 // Constructor; creates a random cloud of droplets
60 // Constructor; creates a random cloud of dro << 60 FastAerosol(const G4String& pName, G4VSolid* pCloud,
61 FastAerosol(const G4String& pName, G4VSolid* << 61 G4double pR, G4double pMinD, G4double pAvgNumDens, G4double pdR,
62 G4double pR, G4double pMinD, << 62 std::function<G4double (G4ThreeVector)> pNumDensDistribution);
63 G4double pAvgNumDens, G4double pdR, <<
64 std::function<G4double (G4ThreeVector)> <<
65 <<
66 FastAerosol(const G4String& pName, G4VSolid* p <<
67 G4double pR, G4double pMinD, <<
68 G4double pNumDens, G4double pdR); <<
69 63
70 FastAerosol(const G4String& pName, G4VSolid* p << 64 FastAerosol(const G4String& pName, G4VSolid* pCloud,
71 G4double pR, G4double pMinD, G4double pN << 65 G4double pR, G4double pMinD, G4double pNumDens, G4double pdR);
>> 66
>> 67 FastAerosol(const G4String& pName, G4VSolid* pCloud,
>> 68 G4double pR, G4double pMinD, G4double pNumDens);
72 69
73 ~FastAerosol()=default; << 70 // Destructor; frees memory
>> 71 ~FastAerosol();
74 72
75 // Populate all grids. Otherwise, they are pop << 73 // Populate all grids. Otherwise, they are populated on-the-fly
76 void PopulateAllGrids(); << 74 void PopulateAllGrids();
77 75
78 // Save locations of droplets to a file for vi << 76 // Save locations of droplets to a file for visualization/analysis purposes
79 void SaveToFile(const char *filename); << 77 void SaveToFile(const char *filename);
80 78
81 // Get absolutely nearest droplet - must be pu << 79 // Get absolutely nearest droplet - must be public as FastAerosolSolid uses it
82 bool GetNearestDroplet(const G4ThreeVector &p, << 80 bool GetNearestDroplet(const G4ThreeVector &p, G4ThreeVector ¢er, G4double &closestDistance, G4double stepLim, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation);
83 81
84 // Get nearest droplet along a vector - must b << 82 // Get nearest droplet along a vector - must be public as FastAerosolSolid uses it
85 bool GetNearestDroplet(const G4ThreeVector &p, << 83 bool GetNearestDroplet(const G4ThreeVector &p, const G4ThreeVector &v, G4ThreeVector ¢er, G4double &closestDistance, G4double stepLim, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation);
86 84
87 // ====== << 85 // ======
88 // Inline << 86 // Inline
89 // ====== << 87 // ======
90 // Input quantities << 88 // Input quantities
91 inline G4String GetName() const; //fasterosol << 89 inline G4String GetName() const; // aerosol name
92 inline G4VSolid* GetBulk() const; // bulk shap << 90 inline G4VSolid* GetBulk() const; // bulk shape
93 inline G4double GetRadius() const; // droplet << 91 inline G4double GetRadius() const; // droplet radius
94 inline G4double GetAvgNumDens() const; // dro << 92 inline G4double GetAvgNumDens() const; // droplet number density
95 //inline G4double GetPitch() const; // grid << 93 //inline G4double GetPitch() const; // grid pitch
96 <<
97 inline G4int GetNumDroplets() const; <<
98 // in case the absolute number is more relevan <<
99 <<
100 // Bulk quantities <<
101 inline G4double GetXHalfLength() const; <<
102 inline G4double GetYHalfLength() const; <<
103 inline G4double GetZHalfLength() const; <<
104 inline void GetBoundingLimits(G4ThreeVector &p <<
105 inline G4double GetCubicVolume() const; <<
106 <<
107 inline G4double DistanceToCloud(const G4ThreeV <<
108 inline G4double DistanceToCloud(const G4ThreeV <<
109 <<
110 // Misc getters and setters <<
111 inline long GetSeed(); <<
112 inline void SetSeed(long seed); <<
113 <<
114 inline G4int GetNumPlacementTries(); <<
115 inline void SetNumPlacementTries(G4int numTrie <<
116 <<
117 inline G4int GetPreSphereR(); <<
118 inline void SetPreSphereR(G4int fPreSphereRIn) <<
119 94
120 inline std::function<G4double (G4ThreeVector)> << 95 inline G4int GetNumDroplets() const; // in case the absolute number is more relevant than density
121 96
122 inline G4double GetDropletsPerVoxel(); << 97 // Bulk quantities
123 inline void SetDropletsPerVoxel(G4double newDr << 98 inline G4double GetXHalfLength() const;
>> 99 inline G4double GetYHalfLength() const;
>> 100 inline G4double GetZHalfLength() const;
>> 101 inline void GetBoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const;
>> 102 inline G4double GetCubicVolume() const;
124 103
125 // Printing diagnostic tool << 104 inline G4double DistanceToCloud(const G4ThreeVector &p);
126 inline void PrintPopulationReport(); << 105 inline G4double DistanceToCloud(const G4ThreeVector &p, const G4ThreeVector &v);
127 <<
128 private: <<
129 G4double kCarTolerance; <<
130 106
131 // Parameters, set in constructor << 107 // Misc getters and setters
132 G4String fName; << 108 inline long GetSeed();
133 G4VSolid* fCloud; // Solid volume of the clou << 109 inline void SetSeed(long seed);
134 G4double fDx, fDy, fDz; // Half widths <<
135 G4double fR; // Bounding radius of each dropl <<
136 G4double fR2; // Bounding radius squared of <<
137 G4double fdR; // Uncertainty in DistanceToIn <<
138 110
139 G4double fMinD; // Minimum distance allowed be << 111 inline G4int GetNumPlacementTries();
>> 112 inline void SetNumPlacementTries(G4int numTries);
140 113
141 std::function<G4double (G4ThreeVector)> fDistr << 114 inline G4int GetPreSphereR();
142 G4double fAvgNumDens; // Average droplet numbe << 115 inline void SetPreSphereR(G4int fPreSphereRIn);
143 116
144 long int fNumDroplets = 0; << 117 inline std::function<G4double (G4ThreeVector)> GetDistribution(); // the droplet number density distribution
145 // Number of droplets that have been created <<
146 118
147 G4double fGridPitch; << 119 inline G4double GetDropletsPerVoxel();
148 // Pitch of collision detection grid. Must be << 120 inline void SetDropletsPerVoxel(G4double newDropletsPerVoxel);
149 121
150 // Ramdom engine << 122 // Printing diagnostic tool
151 CLHEP::HepJamesRandom fCloudEngine; << 123 inline void PrintPopulationReport();
152 long fSeed = 0; // Global random seed << 124
>> 125 private:
>> 126 G4double kCarTolerance;
153 127
154 G4double fDropletsPerVoxel = 4.0; << 128 // Parameters, set in constructor
155 // Expected number of droplets per voxel << 129 G4String fName;
156 130
157 // How far the voxel center must be inside the << 131 G4VSolid* fCloud; // Solid volume of the cloud
158 //order for there to be no risk of placing a d << 132 G4double fDx, fDy, fDz; // Half widths
159 G4double fEdgeDistance; <<
160 133
161 // Grid variables << 134 G4double fR; // Bounding radius of each droplet
162 std::vector<std::vector<G4ThreeVector>> fGrid; << 135 G4double fR2; // Bounding radius squared of each droplet
163 // Grid of lists of inidices to grid points, <<
164 //used for fast collsion checking <<
165 136
166 std::vector<G4double> fGridMean; << 137 G4double fdR; // Uncertainty in DistanceToIn droplet when just using knowledge of droplet center
167 // Array listing mean count for each voxel <<
168 138
169 std::atomic<bool> *fGridValid; << 139 G4double fMinD; // Minimum distance allowed between faces of droplets when constructing random array of droplets
170 // Array listing validity of each grid. uses a <<
171 140
172 G4int fNx, fNy, fNz; // Number of x, y, and z << 141 std::function<G4double (G4ThreeVector)> fDistribution;
>> 142 G4double fAvgNumDens; // Average droplet number density
>> 143 long int fNumDroplets = 0; // Number of droplets that have been created
173 144
174 G4int fNxy; // Cached fNx*fNy << 145 G4double fGridPitch; // Pitch of collision detection grid. Must be greater than diameter of droplets for correctness of collision detection.
175 long int fNumGridCells; // Cached fNx*fNy*fNz <<
176 146
177 G4double fCollisionLimit2; << 147 // Ramdom engine
178 // Threshold distance squared when checking fo << 148 CLHEP::HepJamesRandom fCloudEngine;
179 << 149 long fSeed = 0; // Global random seed
180 G4int fNumNewPointTries = 100; <<
181 // How many times we try to place droplets <<
182 <<
183 G4double fMaxDropPercent = 1.0; <<
184 // The maximal percentage of skipped droplets <<
185 150
186 G4int fMaxDropCount; << 151 G4double fDropletsPerVoxel = 4.0; // Expected number of droplets per voxel
187 // The maximal number of skipped droplets befo <<
188 <<
189 G4int fNumDropped = 0; <<
190 // Number of skipped droplets due to collision <<
191 152
192 G4int fNumCollisions = 0; << 153 // How far the voxel center must be inside the bulk order for there to be no risk of placing a droplet outside
193 // How many collisions occured when attempting << 154 G4double fEdgeDistance;
194 155
195 // Droplet search variables << 156 // Grid variables
196 G4int fVectorSearchRadius; << 157 std::vector<std::vector<G4ThreeVector>> fGrid;// Grid of lists of inidices to grid points, used for fast collsion checking
197 // maximum vector search radius << 158 std::vector<G4double> fGridMean; // Array listing mean count for each voxel
>> 159 std::atomic<bool> *fGridValid; // Array listing validity of each grid. uses atomic variables
198 160
199 // Droplet placement functions << 161 G4int fNx, fNy, fNz; // Number of x, y, and z elements in fGrid
200 // =========================== << 162 G4int fNxy; // Cached fNx*fNy
201 void InitializeGrid(); << 163 long int fNumGridCells; // Cached fNx*fNy*fNz
202 164
203 G4bool FindNewPoint(G4bool edgeVoxel, G4double <<
204 165
205 G4double VoxelOverlap(G4ThreeVector voxelCente << 166 G4double fCollisionLimit2; // Threshold distance squared when checking for collsion
>> 167 G4int fNumNewPointTries = 100; // How many times we try to place droplets
>> 168 G4double fMaxDropPercent = 1.0; // The maximal percentage of skipped droplets before crashing0
>> 169 G4int fMaxDropCount; // The maximal number of skipped droplets before crashing
>> 170 G4int fNumDropped = 0; // Number of skipped droplets due to collisions/out of bulk placement
206 171
207 bool CheckCollision(G4double x, G4double y, G4 << 172 G4int fNumCollisions = 0; // How many collisions occured when attempting to place
208 bool CheckCollisionInsideGrid(G4double x, G4do <<
209 bool CheckCollisionWithDroplet(G4double x, G4d <<
210 173
211 // Droplet distance functions << 174 // Droplet search variables
212 // ========================== << 175 G4int fVectorSearchRadius; // maximum vector search radius
213 void SearchSphere(G4int searchRad, G4double &m <<
214 176
215 void GetNearestDropletInsideRegion(G4double &m << 177 // Droplet placement functions
>> 178 // ===========================
>> 179 void InitializeGrid();
216 180
217 void GetNearestDropletInsideGrid(G4double &min << 181 G4bool FindNewPoint(G4bool edgeVoxel, G4double dX, G4double dY, G4double dZ, G4double minX, G4double minY, G4double minZ, G4ThreeVector &foundVec);
218 182
219 void GetNearestDropletInsideGrid(G4double &min << 183 G4double VoxelOverlap(G4ThreeVector voxelCenter, G4int nStat, G4double epsilon);
220 184
221 // Voxelized sphere methods << 185 bool CheckCollision(G4double x, G4double y, G4double z);
222 // ======================== << 186 bool CheckCollisionInsideGrid(G4double x, G4double y, G4double z, unsigned int xi, unsigned int yi, unsigned int zi);
223 // a collection of points as in {{x1,y1},{x2,y << 187 bool CheckCollisionWithDroplet(G4double x, G4double y, G4double z, G4ThreeVector p);
224 typedef std::vector<std::vector<int>> fCircleT <<
225 188
226 // a collection of points describing a spheric << 189 // Droplet distance functions
227 // with points (x,y,z)=(i-R,j-R,sphere[i][j][k << 190 // ==========================
228 // that is, first index gives x-position, seco << 191 void SearchSphere(G4int searchRad, G4double &minDistance, std::vector<G4ThreeVector> &candidates, std::vector<G4double> &distances2, G4int xGrid, G4int yGrid, G4int zGrid, const G4ThreeVector &p);
229 // gives y-position, and the value gives z-pos << 192 void GetNearestDropletInsideRegion(G4double &minDistance, G4ThreeVector ¢er, int xGrid, int yGrid, int zGrid, int xWidth, int yWidth, int zWidth, const G4ThreeVector &p, const G4ThreeVector &v, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation);
230 // << 193 void GetNearestDropletInsideGrid(G4double &minDistance, std::vector<G4ThreeVector> &candidates, std::vector<G4double> &distances2, unsigned int xGrid, unsigned int yGrid, unsigned int zGrid, const G4ThreeVector &p);
231 // this is done so that searching may be optim << 194 void GetNearestDropletInsideGrid(G4double &minDistance, G4ThreeVector ¢er, unsigned int xGrid, unsigned int yGrid, unsigned int zGrid, const G4ThreeVector &p, const G4ThreeVector &v, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation);
232 // if searching some x=i-R that is outside the <<
233 // aerosol's bounding box, immediately increme <<
234 // (similar for y). <<
235 typedef std::vector<std::vector<std::vector<in <<
236 195
237 G4int fMaxCircleR; << 196 // Voxelized sphere methods
238 G4int fMaxSphereR; << 197 // ========================
239 G4int fPreSphereR = 20; << 198 // a collection of points as in {{x1,y1},{x2,y2},...}
240 std::vector<fCircleType> fCircleCollection; << 199 typedef std::vector<std::vector<int>> fCircleType;
241 std::vector<fSphereType> fSphereCollection; <<
242 fSphereType MakeSphere(G4int R); <<
243 fCircleType MakeCircle(G4int R); <<
244 fCircleType MakeHalfCircle(G4int R); <<
245 200
246 void PopulateGrid(unsigned int xi, unsigned in << 201 // a collection of points describing a spherical shell
>> 202 // with points (x,y,z)=(i-R,j-R,sphere[i][j][k])
>> 203 // that is, first index gives x-position, second index
>> 204 // gives y-position, and the value gives z-position
>> 205 //
>> 206 // this is done so that searching may be optimized:
>> 207 // if searching some x=i-R that is outside the
>> 208 // aerosol's bounding box, immediately increment i
>> 209 // (similar for y).
>> 210 typedef std::vector<std::vector<std::vector<int>>> fSphereType;
247 211
248 // ====== << 212 G4int fMaxCircleR;
249 // Inline << 213 G4int fMaxSphereR;
250 // ====== << 214 G4int fPreSphereR = 20;
251 inline bool GetGrid(const G4ThreeVector &p, G4 << 215 std::vector<fCircleType> fCircleCollection;
>> 216 std::vector<fSphereType> fSphereCollection;
>> 217
>> 218 fSphereType MakeSphere(G4int R);
>> 219 fCircleType MakeCircle(G4int R);
>> 220 fCircleType MakeHalfCircle(G4int R);
252 221
253 inline bool AnyIndexOutOfBounds(G4int xGrid, G << 222 void PopulateGrid(unsigned int xi, unsigned int yi, unsigned int zi, unsigned int& gi);
254 223
255 inline unsigned int GetGridIndex(unsigned int << 224 // ======
>> 225 // Inline
>> 226 // ======
>> 227 inline bool GetGrid(const G4ThreeVector &p, G4int &xGrid, G4int &yGrid, G4int &zGrid);
>> 228 inline bool AnyIndexOutOfBounds(G4int xGrid, G4int yGrid, G4int zGrid);
256 229
257 inline G4ThreeVector GetIndexCoord(G4int index << 230 inline unsigned int GetGridIndex(unsigned int xi, unsigned int yi, unsigned int zi);
>> 231 inline G4ThreeVector GetIndexCoord(G4int index);
258 232
259 inline std::pair<G4int, G4int> GetMinMaxSide(G << 233 inline std::pair<G4int, G4int> GetMinMaxSide(G4int index, G4int numGrids);
260 }; 234 };
261 235
262 #include "FastAerosol.icc" 236 #include "FastAerosol.icc"
263 237
264 #endif 238 #endif
265 239