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