| 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 // Implementation for FastAerosol class 28 // Implementation for FastAerosol class
29 // Author: A.Knaian (ara@nklabs.com), N.MacFad 29 // Author: A.Knaian (ara@nklabs.com), N.MacFadden (natemacfadden@gmail.com)
30 // ------------------------------------------- 30 // --------------------------------------------------------------------
31 31
32 #include "FastAerosol.hh" 32 #include "FastAerosol.hh"
33 33
34 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh"
35 #include "G4GeometryTolerance.hh" 35 #include "G4GeometryTolerance.hh"
36 36
37 #include <numeric> // for summing vectors wit 37 #include <numeric> // for summing vectors with accumulate
38 38
39 // multithreaded safety 39 // multithreaded safety
40 //#include <atomic> 40 //#include <atomic>
41 #include "G4AutoLock.hh" 41 #include "G4AutoLock.hh"
42 namespace 42 namespace
43 { 43 {
44 G4Mutex gridMutex = G4MUTEX_INITIALIZER; << 44 G4Mutex gridMutex = G4MUTEX_INITIALIZER;
45 G4Mutex sphereMutex = G4MUTEX_INITIALIZER; << 45 G4Mutex sphereMutex = G4MUTEX_INITIALIZER;
46 } 46 }
47 47
48 ////////////////////////////////////////////// 48 ///////////////////////////////////////////////////////////////////////////////
49 // 49 //
50 // Constructor - check inputs 50 // Constructor - check inputs
51 // - initialize grid 51 // - initialize grid
52 // - cache values 52 // - cache values
53 // 53 //
54 FastAerosol::FastAerosol(const G4String& pName << 54 FastAerosol::FastAerosol(const G4String& pName,
55 G4double pR, G4double pMinD, G4do << 55 G4VSolid* pCloud,
56 G4double pdR, << 56 G4double pR,
57 std::function<G4double (G4ThreeV << 57 G4double pMinD,
>> 58 G4double pAvgNumDens,
>> 59 G4double pdR,
>> 60 std::function<G4double (G4ThreeVector)> pDistribution)
58 : fName(pName), fCloud(pCloud), fMinD(pMinD) 61 : fName(pName), fCloud(pCloud), fMinD(pMinD), fDistribution(pDistribution)
59 { 62 {
60 kCarTolerance = G4GeometryTolerance::GetInsta << 63 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
61 64
62 // get and set bounding box dimensions << 65
63 G4ThreeVector minBounding, maxBounding; << 66 // get and set bounding box dimensions
64 fCloud->BoundingLimits(minBounding, maxBoundi << 67 G4ThreeVector minBounding, maxBounding;
65 G4ThreeVector halfSizeVec = 0.5*(maxBounding << 68 fCloud->BoundingLimits(minBounding, maxBounding);
66 << 69 G4ThreeVector halfSizeVec = 0.5*(maxBounding - minBounding);
67 G4double pX = halfSizeVec[0]; << 70
68 G4double pY = halfSizeVec[1]; << 71 G4double pX = halfSizeVec[0];
69 G4double pZ = halfSizeVec[2]; << 72 G4double pY = halfSizeVec[1];
70 << 73 G4double pZ = halfSizeVec[2];
71 if (pX < 2*kCarTolerance || << 74
72 pY < 2*kCarTolerance || << 75 if (pX < 2*kCarTolerance ||
73 pZ < 2*kCarTolerance) // limit to thicknes << 76 pY < 2*kCarTolerance ||
>> 77 pZ < 2*kCarTolerance) // limit to thickness of surfaces
74 { 78 {
75 std::ostringstream message; << 79 std::ostringstream message;
76 message << "Dimensions too small for cloud: << 80 message << "Dimensions too small for cloud: " << GetName() << "!" << G4endl
77 << GetName() << "!" << G4endl << 81 << " fDx, fDy, fDz = " << pX << ", " << pY << ", " << pZ;
78 << " fDx, fDy, fDz = " << pX << ", " <<
79 G4Exception("FastAerosol::FastAerosol()", 82 G4Exception("FastAerosol::FastAerosol()", "GeomSolids0002",
80 FatalException, message); 83 FatalException, message);
81 } 84 }
82 fDx = pX; 85 fDx = pX;
83 fDy = pY; 86 fDy = pY;
84 fDz = pZ; 87 fDz = pZ;
85 88
86 89
87 // check and set droplet radius 90 // check and set droplet radius
88 if (pR < 0.0) 91 if (pR < 0.0)
89 { 92 {
90 std::ostringstream message; 93 std::ostringstream message;
91 message << "Invalid radius for cloud: " << 94 message << "Invalid radius for cloud: " << GetName() << "!" << G4endl
92 << " Radius must be positive." 95 << " Radius must be positive."
93 << " Inputs: pR = " << pR; 96 << " Inputs: pR = " << pR;
94 G4Exception("FastAerosol::FastAerosol()", 97 G4Exception("FastAerosol::FastAerosol()", "GeomSolids0002",
95 FatalErrorInArgument, message); 98 FatalErrorInArgument, message);
96 } 99 }
97 fR = pR; 100 fR = pR;
98 fR2 = fR*fR; 101 fR2 = fR*fR;
99 102
100 103
101 // check and set droplet radius safety 104 // check and set droplet radius safety
102 if (pdR < 0.0 || pdR > fR) 105 if (pdR < 0.0 || pdR > fR)
103 { 106 {
104 std::ostringstream message; 107 std::ostringstream message;
105 message << "Invalid sphericality measure f 108 message << "Invalid sphericality measure for cloud: " << GetName() << "!" << G4endl
106 << " Radius uncertainty must be 109 << " Radius uncertainty must be between 0 and fR."
107 << " Inputs: pdR = " << pdR 110 << " Inputs: pdR = " << pdR
108 << " Inputs: pR = " << pR; 111 << " Inputs: pR = " << pR;
109 G4Exception("FastAerosol::FastAerosol()", 112 G4Exception("FastAerosol::FastAerosol()", "GeomSolids0002",
110 FatalErrorInArgument, message); 113 FatalErrorInArgument, message);
111 } 114 }
112 fdR = pdR; 115 fdR = pdR;
113 116
114 117
115 // check and set number density 118 // check and set number density
116 if (pAvgNumDens <= 0) 119 if (pAvgNumDens <= 0)
117 { 120 {
118 std::ostringstream message; 121 std::ostringstream message;
119 message << "Invalid average number density 122 message << "Invalid average number density for cloud: " << GetName() << "!" << G4endl
120 << " pAvgNumDens = " << pAvgNumDen 123 << " pAvgNumDens = " << pAvgNumDens;
121 G4Exception("FastAerosol::FastAerosol()", 124 G4Exception("FastAerosol::FastAerosol()", "GeomSolids0002",
122 FatalException, message); 125 FatalException, message);
123 } 126 }
124 fAvgNumDens = pAvgNumDens; 127 fAvgNumDens = pAvgNumDens;
125 128
126 129
127 // Set collision limit for collsion between 130 // Set collision limit for collsion between equal sized balls with fMinD between them
128 // no droplets will be placed closer than th 131 // no droplets will be placed closer than this distance forom each other
129 fCollisionLimit2 = (2*fR + fMinD)*(2*fR + fM 132 fCollisionLimit2 = (2*fR + fMinD)*(2*fR + fMinD);
130 133
131 // set maximum number of droplet that we are 134 // set maximum number of droplet that we are allowed to skip before halting with an error
132 fMaxDropCount = (G4int)floor(fAvgNumDens*(fC 135 fMaxDropCount = (G4int)floor(fAvgNumDens*(fCloud->GetCubicVolume())*(0.01*fMaxDropPercent));
133 136
134 // initialize grid variables 137 // initialize grid variables
135 InitializeGrid(); 138 InitializeGrid();
136 139
137 140
138 // begin cache of voxelized circles and sphe 141 // begin cache of voxelized circles and spheres
139 // see header for more details on these data 142 // see header for more details on these data structures
140 G4AutoLock lockSphere(&sphereMutex); // loc 143 G4AutoLock lockSphere(&sphereMutex); // lock for multithreaded safety. Likely not needed here, but doesn't hurt
141 144
142 fCircleCollection.push_back({{0, 0}}); // t 145 fCircleCollection.push_back({{0, 0}}); // the R=0 circle only has one point {{x1,y1}} = {{0,0}}
143 fSphereCollection.push_back({{{0}}}); // the 146 fSphereCollection.push_back({{{0}}}); // the R=0 sphere only has one point {{{z1}}} = {{{0}}}
144 147
145 fMaxCircleR = 0; 148 fMaxCircleR = 0;
146 fMaxSphereR = 0; 149 fMaxSphereR = 0;
147 150
148 MakeSphere(fPreSphereR); 151 MakeSphere(fPreSphereR);
149 152
150 lockSphere.unlock(); // unlock 153 lockSphere.unlock(); // unlock
151 154
152 // vector search radius. In terms of voxel w 155 // vector search radius. In terms of voxel width, how far do you search for droplets in vector search
153 // you need to search a larger area if fR is 156 // you need to search a larger area if fR is larger than one grid (currently disabled)
154 fVectorSearchRadius = (G4int)ceill(fR/fGridP 157 fVectorSearchRadius = (G4int)ceill(fR/fGridPitch);
155 } 158 }
156 159
157 ////////////////////////////////////////////// 160 ///////////////////////////////////////////////////////////////////////////////
158 // 161 //
159 // Alternative Constructor 162 // Alternative Constructor
160 // 163 //
161 // Same as standard constructor except with a 164 // Same as standard constructor except with a uniform droplet distribution
162 // 165 //
163 FastAerosol::FastAerosol(const G4String& pName << 166 FastAerosol::FastAerosol(const G4String& pName,
>> 167 G4VSolid* pCloud,
>> 168 G4double pR,
>> 169 G4double pMinD,
>> 170 G4double pNumDens,
>> 171 G4double pdR):
164 FastAerosol(pName, pCloud, pR, pMinD, pNumDe 172 FastAerosol(pName, pCloud, pR, pMinD, pNumDens, pdR,
165 [](G4ThreeVector) {return 1.0;}) 173 [](G4ThreeVector) {return 1.0;})
166 {} 174 {}
167 175
168 ////////////////////////////////////////////// 176 ///////////////////////////////////////////////////////////////////////////////
169 // 177 //
170 // Alternative Constructor 178 // Alternative Constructor
171 // 179 //
172 // Same as standard constructor except with a 180 // Same as standard constructor except with a uniform droplet distribution and
173 // assuming no uncertainty in sphericality 181 // assuming no uncertainty in sphericality
174 // 182 //
175 FastAerosol::FastAerosol(const G4String& pName << 183 FastAerosol::FastAerosol(const G4String& pName,
>> 184 G4VSolid* pCloud,
>> 185 G4double pR,
>> 186 G4double pMinD,
>> 187 G4double pNumDens):
>> 188 FastAerosol(pName, pCloud, pR, pMinD, pNumDens, 0.0,
>> 189 [](G4ThreeVector) {return 1.0;})
176 {} 190 {}
177 191
178 ////////////////////////////////////////////// 192 ///////////////////////////////////////////////////////////////////////////////
179 // 193 //
>> 194 // Destructor
>> 195 //
>> 196 FastAerosol::~FastAerosol() {
>> 197 }
>> 198
>> 199 ///////////////////////////////////////////////////////////////////////////////
>> 200 //
180 // Initialize grid 201 // Initialize grid
181 // 202 //
182 // Sets grids to initial values and calculates 203 // Sets grids to initial values and calculates expected number of droplets
183 // for each voxel 204 // for each voxel
184 // 205 //
185 void FastAerosol::InitializeGrid() << 206 void FastAerosol::InitializeGrid() {
186 { << 207 // set pitch so, on average, fDropletsPerVoxel droplets are in each voxel
187 // set pitch so, on average, fDropletsPerVoxe << 208 fGridPitch = std::pow(fDropletsPerVoxel/fAvgNumDens,1.0/3.0);
188 fGridPitch = std::pow(fDropletsPerVoxel/fAvgN << 209
189 << 210 // if a voxel has center farther than this distance from the bulk outside,
190 // if a voxel has center farther than this di << 211 // we know it is fully contained in the bulk
191 // we know it is fully contained in the bulk << 212 fEdgeDistance = fGridPitch*std::sqrt(3.0)/2.0 + fR;
192 fEdgeDistance = fGridPitch*std::sqrt(3.0)/2.0 << 213
193 << 214 // set number of grid cells
194 // set number of grid cells << 215 fNx = (G4int)ceill(2*fDx / fGridPitch);
195 fNx = (G4int)ceill(2*fDx / fGridPitch); << 216 fNy = (G4int)ceill(2*fDy / fGridPitch);
196 fNy = (G4int)ceill(2*fDy / fGridPitch); << 217 fNz = (G4int)ceill(2*fDz / fGridPitch);
197 fNz = (G4int)ceill(2*fDz / fGridPitch); << 218 fNumGridCells = (long) fNx*fNy*fNz;
198 fNumGridCells = (long) fNx*fNy*fNz; << 219 fNxy = fNx*fNy;
199 fNxy = fNx*fNy; << 220
200 << 221 // try... catch... because vectors can only be so long
201 // try... catch... because vectors can only b << 222 // thus you can't set grid too fine
202 // thus you can't set grid too fine << 223 try
203 try << 224 {
204 { << 225 std::vector<G4ThreeVector> emptyVoxel{};
205 std::vector<G4ThreeVector> emptyVoxel{}; << 226 fGrid.resize(fNumGridCells, emptyVoxel);
206 fGrid.resize(fNumGridCells, emptyVoxel); << 227 fGridMean.resize(fNumGridCells, 0);
207 fGridMean.resize(fNumGridCells, 0); << 228 fGridValid = new std::atomic<bool>[fNumGridCells];
208 fGridValid = new std::atomic<bool>[fNumGrid << 229 }
209 } << 230 catch ( const std::bad_alloc&)
210 catch ( const std::bad_alloc&) << 231 {
211 { << 232 std::ostringstream message;
212 std::ostringstream message; << 233 message << "Out of memory! Grid pitch too small for cloud: " << GetName() << "!" << G4endl
213 message << "Out of memory! Grid pitch too s << 234 << " Asked for fNumGridCells = " << fNumGridCells << G4endl
214 << GetName() << "!" << G4endl << 235 << " each with element of size " << sizeof(std::vector<G4ThreeVector>) << G4endl
215 << " Asked for fNumGridCells = " << 236 << " each with element of size " << sizeof(bool) << G4endl
216 << fNumGridCells << G4endl << 237 << " but the max size is " << fGrid.max_size() << "!";
217 << " each with element of size " << 238 G4Exception("FastAerosol::FastAerosol()", "GeomSolids0002",
218 << sizeof(std::vector<G4ThreeVector>) < << 239 FatalErrorInArgument, message);
219 << " each with element of size " < << 240 }
220 << G4endl <<
221 << " but the max size is " << fGri <<
222 G4Exception("FastAerosol::FastAerosol()", " <<
223 FatalErrorInArgument, message); } <<
224 241
225 // initialize grid cache << 242 // initialize grid cache
226 G4ThreeVector voxelCenter; << 243 G4ThreeVector voxelCenter;
227 G4bool valid; << 244 G4bool valid;
228 << 245
229 for (int i=0; i<fNumGridCells; i++) << 246 for (int i=0; i<fNumGridCells; i++)
230 { <<
231 voxelCenter = fGridPitch*(GetIndexCoord(i)+G <<
232 voxelCenter -= G4ThreeVector(fDx,fDy,fDz); <<
233 <<
234 //shift all voxels so that aerosol center is <<
235 // whether or not the grid is 'finished' <<
236 // fGridValid[i] is true if we don't plan on <<
237 // more droplets in the voxel <<
238 // false if otherwise <<
239 valid = (fCloud->Inside(voxelCenter) != kInsi <<
240 <<
241 if (valid) // will not populate the voxel <<
242 { <<
243 // have to store 'valid' in this way so tha <<
244 fGridValid[i].store(true, std::memory_order <<
245 } <<
246 else // might need to populate the voxel <<
247 { <<
248 fGridValid[i].store(false, std::memory_orde <<
249 <<
250 // find the overlap of the voxel with the b <<
251 G4double volScaling=1.0; <<
252 G4bool edgeVoxel = ( kInside != fCloud->In <<
253 if (edgeVoxel) <<
254 { 247 {
255 volScaling = VoxelOverlap(voxelCenter, 100, << 248 voxelCenter = fGridPitch*(GetIndexCoord(i)+G4ThreeVector(0.5,0.5,0.5)); // center of grid at index i with assuming minimum voxel is at [0,fGridPitch]x[0,fGridPitch]x[0,fGridPitch]
>> 249 voxelCenter -= G4ThreeVector(fDx,fDy,fDz); // shift all voxels so that aerosol center is properly at the origin
>> 250
>> 251 // whether or not the grid is 'finished'
>> 252 // fGridValid[i] is true if we don't plan on populating more droplets in the voxel
>> 253 // false if otherwise
>> 254 valid = (fCloud->Inside(voxelCenter) != kInside && fCloud->DistanceToIn(voxelCenter) >= fEdgeDistance);
>> 255
>> 256 if (valid) // will not populate the voxel
>> 257 {
>> 258 // have to store 'valid' in this way so that the value is atomic and respects multithreadedness
>> 259 fGridValid[i].store(true, std::memory_order_release);
>> 260 }
>> 261 else // might need to populate the voxel
>> 262 {
>> 263 fGridValid[i].store(false, std::memory_order_release);
>> 264
>> 265 // find the overlap of the voxel with the bulk
>> 266 G4double volScaling=1.0;
>> 267 G4bool edgeVoxel = ( kInside != fCloud->Inside(voxelCenter) || fCloud->DistanceToOut(voxelCenter) < fEdgeDistance );
>> 268 if (edgeVoxel)
>> 269 {
>> 270 volScaling = VoxelOverlap(voxelCenter, 100, 0.0);
256 } 271 }
257 // calculates number of droplets based off << 272
258 fGridMean[i] = std::max(0.0, volScaling*fDis << 273 // calculates number of droplets based off of voxel center - not ideal
>> 274 fGridMean[i] = std::max(0.0, volScaling*fDistribution(voxelCenter));
259 } 275 }
260 } 276 }
261 277
262 // must scale fGridMean[i] so that the desired << 278 // must scale fGridMean[i] so that the desired number densities are actualy achieved
263 // this is because no restrictions are applied << 279 // this is because no restrictions are applied to fDistribution
264 G4double tempScaling = (fCloud->GetCubicVolume << 280 G4double tempScaling = (fCloud->GetCubicVolume())*fAvgNumDens/accumulate(fGridMean.begin(), fGridMean.end(), 0.0);
265 for (int i=0; i<fNumGridCells; i++) {fGridMe 281 for (int i=0; i<fNumGridCells; i++) {fGridMean[i] *= tempScaling;}
266 } 282 }
267 283
268 ////////////////////////////////////////////// 284 ///////////////////////////////////////////////////////////////////////////////
269 // 285 //
270 // Calculate the overlap of the voxel with the 286 // Calculate the overlap of the voxel with the aerosol bulk
271 // 287 //
272 // The method is largely copied from G4VSolid: 288 // The method is largely copied from G4VSolid::EstimateCubicVolume
273 // 289 //
274 G4double FastAerosol::VoxelOverlap(G4ThreeVect << 290 G4double FastAerosol::VoxelOverlap(G4ThreeVector voxelCenter, G4int nStat, G4double epsilon) {
275 { << 291 G4double cenX = voxelCenter.x();
276 G4double cenX = voxelCenter.x(); << 292 G4double cenY = voxelCenter.y();
277 G4double cenY = voxelCenter.y(); << 293 G4double cenZ = voxelCenter.z();
278 G4double cenZ = voxelCenter.z(); <<
279 294
280 G4int iInside=0; << 295 G4int iInside=0;
281 G4double px,py,pz; << 296 G4double px,py,pz;
282 G4ThreeVector p; << 297 G4ThreeVector p;
283 G4bool in; << 298 G4bool in;
284 << 299
285 if(nStat < 100) nStat = 100; << 300 if(nStat < 100) nStat = 100;
286 if(epsilon > 0.01) epsilon = 0.01; << 301 if(epsilon > 0.01) epsilon = 0.01;
287 << 302
288 for(G4int i = 0; i < nStat; i++ ) << 303 for(G4int i = 0; i < nStat; i++ )
289 { << 304 {
290 px = cenX+(fGridPitch-epsilon)*(G4UniformRan << 305 px = cenX+(fGridPitch-epsilon)*(G4UniformRand()-0.5);
291 py = cenY+(fGridPitch-epsilon)*(G4UniformRan << 306 py = cenY+(fGridPitch-epsilon)*(G4UniformRand()-0.5);
292 pz = cenZ+(fGridPitch-epsilon)*(G4UniformRan << 307 pz = cenZ+(fGridPitch-epsilon)*(G4UniformRand()-0.5);
293 p = G4ThreeVector(px,py,pz); << 308 p = G4ThreeVector(px,py,pz);
294 in = (fCloud->Inside(p) == kInside) && (fClo << 309 in = (fCloud->Inside(p) == kInside) && (fCloud->DistanceToOut(p) >= fR);
295 if(in) iInside++; << 310 if(in) iInside++;
296 } << 311 }
297 << 312
298 return (double)iInside/nStat; << 313 return (double)iInside/nStat;
299 } 314 }
300 315
>> 316
301 ////////////////////////////////////////////// 317 ///////////////////////////////////////////////////////////////////////////////
302 // 318 //
303 // Populate all voxels 319 // Populate all voxels
304 // 320 //
305 // Allows for partially populated clouds. 321 // Allows for partially populated clouds.
306 // 322 //
307 void FastAerosol::PopulateAllGrids() << 323 void FastAerosol::PopulateAllGrids() {
308 { << 324 unsigned int gi;
309 unsigned int gi; <<
310 325
311 for (int xi=0; xi<fNx; xi++) << 326 for (int xi=0; xi<fNx; xi++)
312 { <<
313 for (int yi=0; yi<fNy; yi++) <<
314 { <<
315 for (int zi=0; zi<fNz; zi++) <<
316 { 327 {
317 // PopulateGrid takes unsigned arguments << 328 for (int yi=0; yi<fNy; yi++)
318 // fNx, fNy, and fNz are signed so that s << 329 {
319 // thus we iterate with xi, yi, and zi sig << 330 for (int zi=0; zi<fNz; zi++)
320 PopulateGrid((unsigned)xi,(unsigned)yi,(un << 331 {
321 } << 332 // PopulateGrid takes unsigned arguments
322 } << 333 // fNx, fNy, and fNz are signed so that some other algebra works
323 } << 334 // thus we iterate with xi, yi, and zi signed and then cast them
>> 335 PopulateGrid((unsigned)xi,(unsigned)yi,(unsigned)zi,gi);
>> 336 }
>> 337 }
>> 338 }
324 } 339 }
325 340
326 ////////////////////////////////////////////// 341 ///////////////////////////////////////////////////////////////////////////////
327 // 342 //
328 // Populate a single grid 343 // Populate a single grid
329 // 344 //
330 // If grid is already populated, does nothing. 345 // If grid is already populated, does nothing.
331 // 346 //
332 void FastAerosol::PopulateGrid(unsigned int xi << 347 void FastAerosol::PopulateGrid(unsigned int xi, unsigned int yi, unsigned int zi, unsigned int& gi) {
333 { << 348 gi = GetGridIndex(xi,yi,zi);
334 gi = GetGridIndex(xi,yi,zi); <<
335 349
336 // Check if this grid needs update << 350 // Check if this grid needs update
337 bool tmpValid = fGridValid[gi].load(std::memo << 351 bool tmpValid = fGridValid[gi].load(std::memory_order_acquire); // have to do this weirdly because we are using atomic variables so that multithreadedness is safe
338 // have to do this weirdly because we are usi << 352 std::atomic_thread_fence(std::memory_order_acquire);
339 std::atomic_thread_fence(std::memory_order_ac <<
340 353
341 if (!tmpValid) // if true then either outside << 354 if (!tmpValid) // if true then either outside the bulk or in a voxel that is already populated
342 { << 355 {
343 G4AutoLock lockGrid(&gridMutex); << 356 G4AutoLock lockGrid(&gridMutex);
344 // Check again now that we have the lock - in << 357
345 tmpValid = fGridValid[gi].load(std::memory_o << 358 // Check again now that we have the lock - in case grid became valid after first check and before lock
346 << 359 tmpValid = fGridValid[gi].load(std::memory_order_acquire);
347 if (!tmpValid) << 360
348 { << 361 if (!tmpValid)
349 // uniquely set the seed to randomly place dro << 362 {
350 // changing global seed gaurantees a totally n << 363 // uniquely set the seed to randomly place droplets
351 fCloudEngine.setSeed((long)gi + (long)fNum << 364 // changing global seed gaurantees a totally new batch of seeds
352 << 365 fCloudEngine.setSeed((long)gi + (long)fNumGridCells*(long)fSeed);
353 // find if the voxel is near the bulk edge. In << 366
354 // if not an edge voxel, we can place droplets << 367 // find if the voxel is near the bulk edge. In that case, we need to check whether a placed droplet is inside the bulk
355 G4ThreeVector voxelCenter = fGridPitch*G4Thr << 368 // if not an edge voxel, we can place droplets by only considering interference between with other droplets
356 G4bool edgeVoxel = ( kInside != fCloud->Insid << 369 G4ThreeVector voxelCenter = fGridPitch*G4ThreeVector(xi+0.5,yi+0.5,zi+0.5)-G4ThreeVector(fDx,fDy,fDz);
>> 370 G4bool edgeVoxel = ( kInside != fCloud->Inside(voxelCenter) || fCloud->DistanceToOut(voxelCenter) < fEdgeDistance );
357 371
358 // number of droplets to place << 372 // number of droplets to place
359 unsigned int numDropletsToPlace = CLHEP::Ran << 373 unsigned int numDropletsToPlace = CLHEP::RandPoisson::shoot(&fCloudEngine, fGridMean[gi]);
>> 374
>> 375 // actually add the points
>> 376 G4ThreeVector point;
>> 377
>> 378 while (numDropletsToPlace > 0)
>> 379 {
>> 380 // find a new point not overlapping with the others
>> 381 if (FindNewPoint(edgeVoxel, fGridPitch, fGridPitch, fGridPitch, (G4double)xi*fGridPitch-fDx, (G4double)yi*fGridPitch-fDy, (G4double)zi*fGridPitch-fDz, point) )
>> 382 {
>> 383 fGrid[gi].push_back(point);
>> 384 fNumDroplets++;
>> 385 }
360 386
361 // actually add the points << 387
362 G4ThreeVector point; << 388 numDropletsToPlace--;
>> 389 }
363 390
364 while (numDropletsToPlace > 0) << 391 // Memory fence to ensure sequential consistency,
365 { << 392 // because fGridValid is read without a lock
366 // find a new point not overlapping with th << 393 // Voxel data update must be complete before updating fGridValid
367 if (FindNewPoint(edgeVoxel, fGridPitch, fGr <<
368 { <<
369 fGrid[gi].push_back(point); <<
370 fNumDroplets++; <<
371 } <<
372 numDropletsToPlace--; <<
373 } <<
374 <<
375 // Memory fence to ensure sequential consisten <<
376 // because fGridValid is read without a lock <<
377 // Voxel data update must be complete before u <<
378 std::atomic_thread_fence(std::memory_ord 394 std::atomic_thread_fence(std::memory_order_release);
379 395
380 fGridValid[gi].store(true, std::memory_order_r << 396 fGridValid[gi].store(true, std::memory_order_release); // set the grid as populated
381 } << 397 }
382 398
383 lockGrid.unlock(); << 399 lockGrid.unlock();
384 } << 400 }
385 } 401 }
386 402
>> 403
387 ////////////////////////////////////////////// 404 ///////////////////////////////////////////////////////////////////////////////
388 // 405 //
389 // Find a new droplet position in a box of ful 406 // Find a new droplet position in a box of full (not half) widths dX, dY, and dZ
390 // with minimum corner at (minX, minY, minZ). 407 // with minimum corner at (minX, minY, minZ).
391 // 408 //
392 // ------------------------------------------- 409 // ----------------------------------------------------------------------------
393 // 410 //
394 // Ex: FindNewPoint(fGridPitch, fGridPitch, fG 411 // Ex: FindNewPoint(fGridPitch, fGridPitch, fGridPitch,
395 // (G4double)xi*fGridPitch-f 412 // (G4double)xi*fGridPitch-fDx, (G4double)yi*fGridPitch-fDy,
396 // (G4double)zi*fGridPitch-f 413 // (G4double)zi*fGridPitch-fDz);
397 // 414 //
398 // 1) CLHEP::RandFlat::shoot(&fCloudEngine)*fG 415 // 1) CLHEP::RandFlat::shoot(&fCloudEngine)*fGridPitch shoots a point in range
399 // [0, fGridPitch] 416 // [0, fGridPitch]
400 // 417 //
401 // 1.5) Note that the minimum x value of the ( 418 // 1.5) Note that the minimum x value of the (xi, yi, zi) grid is
402 // xi*fGridPitch-fDx 419 // xi*fGridPitch-fDx
403 // 420 //
404 // 2) xi*fGridPitch-fDx + (1) adds the minimum 421 // 2) xi*fGridPitch-fDx + (1) adds the minimum x value of the (xi, yi, zi)
405 // voxel 422 // voxel
406 // 423 //
407 // ------------------------------------------- 424 // ----------------------------------------------------------------------------
408 // 425 //
409 // Ex: FindNewPoint(2.0*fDx, 2.0*fDy, 2.0*fDz, 426 // Ex: FindNewPoint(2.0*fDx, 2.0*fDy, 2.0*fDz, -fDx, -fDy, -fDz);
410 // 427 //
411 // 1) CLHEP::RandFlat::shoot(&fCloudEngine)*2. 428 // 1) CLHEP::RandFlat::shoot(&fCloudEngine)*2.0*fDx shoots a point in
412 // range [0, 2.0*fDx] 429 // range [0, 2.0*fDx]
413 // 430 //
414 // 2) add -fDx to change range into [-fDx, fDx 431 // 2) add -fDx to change range into [-fDx, fDx]
415 // 432 //
416 G4bool FastAerosol::FindNewPoint(G4bool edgeVo << 433 G4bool FastAerosol::FindNewPoint(G4bool edgeVoxel, G4double dX, G4double dY, G4double dZ, G4double minX, G4double minY, G4double minZ, G4ThreeVector &foundPt) {
417 { << 434 G4int tries = 0; // counter of tries. Give up after fNumNewPointTries (you likely asked for too dense in this case)
418 G4int tries = 0; << 435
419 // counter of tries. Give up after fNumNewPoi << 436 G4double x, y, z;
420 << 437
421 G4double x, y, z; << 438 G4ThreeVector point;
422 << 439 G4bool placedOutside = false; // variable whether or not the point was placed outside. Used for edgeVoxel checking
423 G4ThreeVector point; << 440
424 G4bool placedOutside = false; << 441 // Generate a droplet and check if it overlaps with existing droplets
425 // variable whether or not the point was plac << 442 do {
426 << 443 tries++;
427 // Generate a droplet and check if it overlap << 444
428 do { << 445 if (tries > fNumNewPointTries) // skip if we tried more than fNumNewPointTries
429 tries++; << 446 {
430 if (tries > fNumNewPointTries) << 447 fNumDropped++;
431 // skip if we tried more than fNumNewPoi << 448 if (fNumDropped < fMaxDropCount) // return error if we skipped more than fMaxDropCount droplets
432 { << 449 {
433 fNumDropped++; << 450 return false;
434 if (fNumDropped < fMaxDropCount) << 451 }
435 // return error if we skipped more than fMa << 452
436 { << 453 std::ostringstream message;
437 return false; << 454 message << "Threw out too many droplets for cloud: " << GetName() << G4endl
438 } << 455 << " Tried to place individual droplest " << fNumNewPointTries << " times." << G4endl
439 << 456 << " This failed for " << fMaxDropCount << " droplets.";
440 std::ostringstream message; << 457 G4Exception("FastAerosol::FindNewPoint()", "GeomSolids0002",
441 message << "Threw out too many droplets for << 458 FatalErrorInArgument, message);
442 << GetName() << G4endl << 459 }
443 << " Tried to place individual drop << 460
444 << fNumNewPointTries << " times." << G4end << 461 x = minX + CLHEP::RandFlat::shoot(&fCloudEngine)*dX;
445 << " This failed for " << fMaxDropC << 462 y = minY + CLHEP::RandFlat::shoot(&fCloudEngine)*dY;
446 << " droplets."; << 463 z = minZ + CLHEP::RandFlat::shoot(&fCloudEngine)*dZ;
447 G4Exception("FastAerosol::FindNewPoint()", << 464
448 FatalErrorInArgument, message); << 465 if (edgeVoxel)
449 } << 466 {
450 << 467 point = G4ThreeVector(x,y,z);
451 x = minX + CLHEP::RandFlat::shoot(&fCloudE << 468 placedOutside = (fCloud->Inside(point) != kInside) || (fCloud->DistanceToOut(point) < fR);
452 y = minY + CLHEP::RandFlat::shoot(&fCloudE <<
453 z = minZ + CLHEP::RandFlat::shoot(&fCloudE <<
454 <<
455 if (edgeVoxel) <<
456 { <<
457 point = G4ThreeVector(x,y,z); <<
458 placedOutside = (fCloud->Inside(point) != <<
459 } 469 }
460 } 470 }
461 while (CheckCollision(x,y,z) || placedOutsid 471 while (CheckCollision(x,y,z) || placedOutside);
462 472
463 foundPt = G4ThreeVector(x,y,z); 473 foundPt = G4ThreeVector(x,y,z);
464 return true; 474 return true;
465 } 475 }
466 476
467 ////////////////////////////////////////////// 477 ///////////////////////////////////////////////////////////////////////////////
468 // 478 //
469 // Get absolutely nearest droplet 479 // Get absolutely nearest droplet
470 // 480 //
471 // Maximum search radius is maxSearch 481 // Maximum search radius is maxSearch
472 // 482 //
473 // Start in grid of point p, check if any sper 483 // Start in grid of point p, check if any speres are in there
474 // and then go out one more (spherically-shape 484 // and then go out one more (spherically-shaped) level to the surrounding grids
475 // and so on, until the whole volume has been 485 // and so on, until the whole volume has been checked or a droplet has been found
476 // in general as you go out, you always need t 486 // in general as you go out, you always need to check one more level out to make sure that
477 // the one you have (at the closer level) is t 487 // the one you have (at the closer level) is the actually the nearest one
478 // 488 //
479 bool FastAerosol::GetNearestDroplet(const G4Th 489 bool FastAerosol::GetNearestDroplet(const G4ThreeVector &p, G4ThreeVector ¢er, G4double &minRealDistance, G4double maxSearch, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation)
480 { 490 {
481 G4double cloudDistance = fCloud->DistanceToIn << 491 G4double cloudDistance = fCloud->DistanceToIn(p);
>> 492
>> 493 // starting radius/diam of voxel layer
>> 494 G4int searchRad = (G4int)floor(0.5+cloudDistance/fGridPitch); // no reason to search shells totally outside cloud
482 495
483 // starting radius/diam of voxel layer << 496 // find the voxel containing p
484 G4int searchRad = (G4int)floor(0.5+cloudDista << 497 int xGrid, yGrid, zGrid;
>> 498 GetGrid(p, xGrid, yGrid, zGrid); // it is OK if the starting point is outside the volume
485 499
486 // find the voxel containing p << 500 // initialize the containers for all candidate closest droplets and their respective distances
487 int xGrid, yGrid, zGrid; << 501 std::vector<G4ThreeVector> candidates;
488 GetGrid(p, xGrid, yGrid, zGrid); << 502 std::vector<G4double> distances;
489 // it is OK if the starting point is outside <<
490 <<
491 // initialize the containers for all candidat <<
492 std::vector<G4ThreeVector> candidates; <<
493 std::vector<G4double> distances; <<
494 // initialize distances to indicate that no d <<
495 G4double minDistance = kInfinity; <<
496 G4double maxCheckDistance = kInfinity; <<
497 <<
498 bool unsafe = true; <<
499 <<
500 while (unsafe) <<
501 { <<
502 // limit maximum search radius <<
503 if (searchRad*fGridPitch > maxSearch+fGridPi <<
504 503
505 SearchSphere(searchRad, minDistance, candid << 504 // initialize distances to indicate that no droplet has yet been found
506 maxCheckDistance = minDistance+fdR; << 505 G4double minDistance = kInfinity;
>> 506 G4double maxCheckDistance = kInfinity;
507 507
508 // theory says that, to safely have searched << 508 bool unsafe = true;
>> 509
>> 510
>> 511 while (unsafe)
>> 512 {
>> 513 // limit maximum search radius
>> 514 if (searchRad*fGridPitch > maxSearch+fGridPitch) { return false; }
>> 515
>> 516 SearchSphere(searchRad, minDistance, candidates, distances, xGrid, yGrid, zGrid, p);
>> 517 maxCheckDistance = minDistance+fdR;
>> 518
>> 519 // theory says that, to safely have searched for centers up to some maxCheckDistance, we must search voxelized spheres at least up to ceil(0.25+maxCheckDistance/fGridPitch)
509 // *** unsure if fully accurate. Calculati 520 // *** unsure if fully accurate. Calculations for 2D, not 3D... ***
510 unsafe = searchRad < std::ceil(0.25+maxCheckD << 521 unsafe = searchRad < std::ceil(0.25+maxCheckDistance/fGridPitch);
511 522
512 searchRad++; << 523 searchRad++;
513 } << 524 }
514 525
515 // delete any collected droplets that are too <<
516 unsigned int index = 0; <<
517 526
518 G4ThreeVector tempCenter; << 527 // delete any collected droplets that are too far out. Check all candidates to see what is closest
519 G4double tempDistance; << 528 unsigned int index = 0;
520 <<
521 minRealDistance = kInfinity; <<
522 <<
523 while (index < distances.size()) <<
524 { <<
525 if (distances[index]>maxCheckDistance) <<
526 { <<
527 candidates.erase(candidates.begin()+index) <<
528 distances.erase(distances.begin()+index); <<
529 } <<
530 else <<
531 { <<
532 tempCenter = candidates[index]; <<
533 tempDistance = droplet->DistanceToIn( rotat <<
534 <<
535 if (tempDistance < minRealDistance) <<
536 { <<
537 minRealDistance = tempDistance; <<
538 center = tempCenter; <<
539 } <<
540 index++; // move onto next droplet <<
541 } <<
542 } <<
543 529
544 return true; << 530 G4ThreeVector tempCenter;
>> 531 G4double tempDistance;
>> 532
>> 533 minRealDistance = kInfinity;
>> 534
>> 535
>> 536 while (index < distances.size())
>> 537 {
>> 538
>> 539 if (distances[index]>maxCheckDistance)
>> 540 {
>> 541 candidates.erase(candidates.begin()+index);
>> 542 distances.erase(distances.begin()+index);
>> 543 }
>> 544 else
>> 545 {
>> 546 tempCenter = candidates[index];
>> 547 tempDistance = droplet->DistanceToIn( rotation(tempCenter).inverse()*(p - tempCenter) );
>> 548
>> 549 if (tempDistance < minRealDistance)
>> 550 {
>> 551 minRealDistance = tempDistance;
>> 552 center = tempCenter;
>> 553 }
>> 554
>> 555 index++; // move onto next droplet
>> 556 }
>> 557 }
>> 558
>> 559 return true;
545 } 560 }
546 561
547 ////////////////////////////////////////////// 562 ///////////////////////////////////////////////////////////////////////////////
548 // 563 //
549 // Search sphere of radius R for droplets 564 // Search sphere of radius R for droplets
550 // 565 //
551 // Saves any droplets to "center" if they are 566 // Saves any droplets to "center" if they are closer than minDistance
552 // (updating this minimum found distance at ea 567 // (updating this minimum found distance at each interval).
553 // 568 //
554 // Returns if minDistance<infinity (i.e., if w 569 // Returns if minDistance<infinity (i.e., if we have yet found a droplet)
555 // 570 //
556 void FastAerosol::SearchSphere(G4int searchRad 571 void FastAerosol::SearchSphere(G4int searchRad, G4double &minDistance, std::vector<G4ThreeVector> &candidates, std::vector<G4double> &distances, G4int xGrid, G4int yGrid, G4int zGrid, const G4ThreeVector &p)
557 { 572 {
558 // search variables << 573 // search variables
559 G4int xSearch, ySearch, zSearch; << 574 G4int xSearch, ySearch, zSearch; // x, y, and z coordinates of the currently searching voxel in the shell
560 // x, y, and z coordinates of the currently s << 575
561 // in the shell voxel layer variables << 576 // voxel layer variables
562 fSphereType shell; // the shell that we are s << 577 fSphereType shell; // the shell that we are searching for droplets
563 <<
564 // we pre-calculate spheres up to radius fPreS <<
565 // any sphere smaller than that does not need <<
566 if (searchRad > fPreSphereR) <<
567 { <<
568 // need to consider making a sphere <<
569 G4AutoLock lockSphere(&sphereMutex); <<
570 if (searchRad > fMaxSphereR) <<
571 { <<
572 // actually need to make a sphere <<
573 shell = MakeSphere(searchRad); <<
574 } <<
575 else <<
576 { <<
577 // we have previously made a sphere large <<
578 shell = fSphereCollection[searchRad]; <<
579 } <<
580 lockSphere.unlock(); <<
581 } <<
582 else <<
583 { <<
584 shell = fSphereCollection[searchRad]; <<
585 } <<
586 <<
587 // search spherical voxel layer <<
588 for (int i=0; i<(2*searchRad+1); i++) <<
589 { <<
590 xSearch = xGrid+i-searchRad; <<
591 if (0<=xSearch && xSearch<fNx) <<
592 { <<
593 for (int j=0; j<(2*searchRad+1); j++) <<
594 { <<
595 ySearch = yGrid+j-searchRad; <<
596 578
597 if (0<=ySearch && ySearch<fNy) << 579 // we pre-calculate spheres up to radius fPreSphereR to speed up calculations
>> 580 // any sphere smaller than that does not need to use locks
>> 581 if (searchRad > fPreSphereR)
598 { 582 {
599 for (auto it = shell[i][j].begin(); it != << 583 // need to consider making a sphere
600 zSearch = *it+zGrid; << 584 G4AutoLock lockSphere(&sphereMutex);
601 585
602 if (0<=zSearch && zSearch<fNz) << 586 if (searchRad > fMaxSphereR) {
603 { GetNearestDropletInsideGrid( << 587 // actually need to make a sphere
604 } << 588 shell = MakeSphere(searchRad);
>> 589 }
>> 590 else
>> 591 {
>> 592 // we have previously made a sphere large enough
>> 593 shell = fSphereCollection[searchRad];
>> 594 }
>> 595
>> 596 lockSphere.unlock();
605 } 597 }
>> 598 else
>> 599 {
>> 600 shell = fSphereCollection[searchRad];
>> 601 }
>> 602
>> 603 // search spherical voxel layer
>> 604 for (int i=0; i<(2*searchRad+1); i++) {
>> 605 xSearch = xGrid+i-searchRad;
>> 606
>> 607 if (0<=xSearch && xSearch<fNx)
>> 608 {
>> 609 for (int j=0; j<(2*searchRad+1); j++) {
>> 610 ySearch = yGrid+j-searchRad;
>> 611
>> 612 if (0<=ySearch && ySearch<fNy)
>> 613 {
>> 614 for (auto it = shell[i][j].begin(); it != shell[i][j].end(); ++it) {
>> 615 zSearch = *it+zGrid;
>> 616
>> 617 if (0<=zSearch && zSearch<fNz)
>> 618 {
>> 619 GetNearestDropletInsideGrid(minDistance, candidates, distances, (unsigned)xSearch, (unsigned)ySearch, (unsigned)zSearch, p);
>> 620 }
>> 621 }
>> 622 }
>> 623 }
>> 624 }
606 } 625 }
607 } <<
608 } <<
609 } <<
610 } 626 }
611 627
612 ////////////////////////////////////////////// 628 ///////////////////////////////////////////////////////////////////////////////
613 // 629 //
614 // Make voxelized spheres up to radius R for f 630 // Make voxelized spheres up to radius R for fSphereCollection
615 // 631 //
616 // These spheres are used for searching for dr 632 // These spheres are used for searching for droplets
617 // 633 //
618 // To create them, we first create a helper ha 634 // To create them, we first create a helper half-circle "zr" for which each
619 // voxel/point in zr represents a layer of our 635 // voxel/point in zr represents a layer of our sphere. The 2nd coordinate of
620 // zr represents the radius of the layer and t 636 // zr represents the radius of the layer and the 1st coordinate denotes the
621 // layer's displacement from the origin along 637 // layer's displacement from the origin along the z-axis (layers are perp to
622 // z in our convention). 638 // z in our convention).
623 // 639 //
624 std::vector<std::vector<std::vector<int>>> Fas 640 std::vector<std::vector<std::vector<int>>> FastAerosol::MakeSphere(G4int R) {
625 // inductively make smaller spheres since fSph << 641 // inductively make smaller spheres since fSphereCollection organizes by index
626 if (fMaxSphereR<(R-1)) { MakeSphere(R-1);} << 642 if (fMaxSphereR<(R-1)) {
627 << 643 MakeSphere(R-1);
628 std::vector<int> z; << 644 }
629 // z-coordinates of voxels in (x,y) of sphere <<
630 <<
631 std::vector<std::vector<int>> y(2*R+1, z); <<
632 // plane YZ of sphere <<
633 <<
634 fSphereType x(2*R+1, y); // entire sphere <<
635 <<
636 x[R][R].push_back(R);// add (0,0,R) <<
637 x[R][R].push_back(-R);// add (0,0,-R) <<
638 fCircleType zr = MakeHalfCircle(R); <<
639 // break sphere into a collection of circles c <<
640 <<
641 for (auto it1 = zr.begin(); it1 != zr.end(); + <<
642 { <<
643 std::vector<int> pt1 = *it1; <<
644 fCircleType circ = MakeCircle(pt1[1]); // mak <<
645 645
646 for (auto it2 = circ.begin(); it2 != circ.en << 646 std::vector<int> z; // z-coordinates of voxels in (x,y) of sphere
647 { << 647 std::vector<std::vector<int>> y(2*R+1, z); // plane YZ of sphere
648 std::vector<int> pt2 = *it2; << 648 fSphereType x(2*R+1, y); // entire sphere
649 x[pt2[0]+R][pt2[1]+R].push_back(pt1[0]); << 649
650 } << 650 x[R][R].push_back(R); // add (0,0,R)
651 } << 651 x[R][R].push_back(-R); // add (0,0,-R)
>> 652
>> 653 fCircleType zr = MakeHalfCircle(R); // break sphere into a collection of circles centered at (0,0,z) for different -R<=z<=R
>> 654
>> 655 for (auto it1 = zr.begin(); it1 != zr.end(); ++it1) {
>> 656 std::vector<int> pt1 = *it1;
>> 657 fCircleType circ = MakeCircle(pt1[1]); // make the circle
>> 658
>> 659 for (auto it2 = circ.begin(); it2 != circ.end(); ++it2) {
>> 660 std::vector<int> pt2 = *it2;
>> 661 x[pt2[0]+R][pt2[1]+R].push_back(pt1[0]);
>> 662 }
>> 663 }
652 664
653 fSphereCollection.push_back(x); << 665 fSphereCollection.push_back(x);
654 fMaxSphereR = R; << 666 fMaxSphereR = R;
655 return x; << 667 return x;
656 } 668 }
657 669
658 ////////////////////////////////////////////// 670 ///////////////////////////////////////////////////////////////////////////////
659 // 671 //
660 // Make voxelized circles up to radius R for f 672 // Make voxelized circles up to radius R for fCircleCollection
661 // 673 //
662 // These circles are used to create voxelized 674 // These circles are used to create voxelized spheres (these spheres are then
663 // used for searching for droplets). This just 675 // used for searching for droplets). This just ports the calculation to making
664 // a half-circle, adding the points (R,0) and 676 // a half-circle, adding the points (R,0) and (0,R), and reflecting the half-
665 // circle along the x-axis. 677 // circle along the x-axis.
666 // 678 //
667 std::vector<std::vector<int>> FastAerosol::Mak << 679 std::vector<std::vector<int>> FastAerosol::MakeCircle(G4int R) {
668 { << 680 // inductively make smaller circles since fCircleCollection organizes by index
669 // inductively make smaller circles since fCir << 681 if (fMaxCircleR<R) {
670 if (fMaxCircleR<R) << 682 if (fMaxCircleR<(R-1)) {
671 { << 683 MakeCircle(R-1);
672 if (fMaxCircleR<(R-1)) {MakeCircle(R-1);} << 684 }
673 fCircleType voxels = {{R, 0}, {-R, 0}}; << 685
674 // add (R,0) and (-R,0) since half circle e << 686 fCircleType voxels = {{R, 0}, {-R, 0}}; // add (R,0) and (-R,0) since half circle excludes them
675 fCircleType halfVoxels = MakeHalfCircle(R); << 687 fCircleType halfVoxels = MakeHalfCircle(R);
676 << 688
677 // join the voxels, halfVoxels, and -halfVo << 689 // join the voxels, halfVoxels, and -halfVoxels
678 for (auto it = halfVoxels.begin(); it != ha << 690 for (auto it = halfVoxels.begin(); it != halfVoxels.end(); ++it) {
679 { << 691 std::vector<int> pt = *it;
680 std::vector<int> pt = *it; << 692 voxels.push_back(pt);
681 voxels.push_back(pt); << 693 voxels.push_back({-pt[0], -pt[1]});
682 voxels.push_back({-pt[0], -pt[1]}); << 694 }
683 } << 695
684 fCircleCollection.push_back(voxels); << 696 fCircleCollection.push_back(voxels);
685 fMaxCircleR++; << 697 fMaxCircleR++;
686 } << 698 }
687 699
688 return fCircleCollection[R]; << 700 return fCircleCollection[R];
689 } 701 }
690 702
691 ////////////////////////////////////////////// 703 ///////////////////////////////////////////////////////////////////////////////
692 // 704 //
693 // Make half circle by a modified midpoint cir 705 // Make half circle by a modified midpoint circle algorithm
694 // 706 //
695 // Modifies the algorithm so it doesn't have ' 707 // Modifies the algorithm so it doesn't have 'holes' when making many circles
696 // 708 //
697 // See B. Roget, J. Sitaraman, "Wall distance 709 // See B. Roget, J. Sitaraman, "Wall distance search algorithim using
698 // voxelized marching spheres," Journal of Com 710 // voxelized marching spheres," Journal of Computational Physics, Vol. 241,
699 // May 2013, pp. 76-94, https://doi.org/10.101 711 // May 2013, pp. 76-94, https://doi.org/10.1016/j.jcp.2013.01.035
700 // 712 //
701 std::vector<std::vector<int>> FastAerosol::Mak << 713 std::vector<std::vector<int>> FastAerosol::MakeHalfCircle(G4int R) {
702 { << 714 // makes an octant of a voxelized circle in the 1st quadrant starting at y-axis
703 // makes an octant of a voxelized circle in th << 715 fCircleType voxels = {{0, R}}; // hard code inclusion of (0,R)
704 fCircleType voxels = {{0, R}}; // hard code i << 716 int x = 1;
705 int x = 1; << 717 int y = R;
706 int y = R; << 718
707 int dx = 4-R; << 719 int dx = 4-R; // measure of whether (x+1,y) will be outside the sphere
708 // measure of whether (x+1,y) will be outside << 720 int dxup = 1; // measure of whether (x+1,y+1) will be outside the sphere
709 int dxup = 1; << 721 while (x<y) {
710 // measure of whether (x+1,y+1) will be outsi << 722 // mirror the points to be added to change octant->upper half plane
711 while (x<y) << 723 voxels.push_back({ x, y});
712 { << 724 voxels.push_back({ y, x});
713 // mirror the points to be added to change o << 725 voxels.push_back({-x, y});
714 voxels.push_back({ x, y}); << 726 voxels.push_back({-y, x});
715 voxels.push_back({ y, x}); << 727
716 voxels.push_back({-x, y}); << 728 // if next point outside the sphere, de-increment y
717 voxels.push_back({-y, x}); << 729 if (dx>0) {
718 // if next point outside the sphere, de-incr << 730 // if dxup<=0, the layer above just moves to the right
719 if (dx>0) << 731 // this would create a hole at (x+1,y) unless we add it
720 { << 732 dxup = dx + 2*y -2*R-1;
721 // if dxup<=0, the layer above just moves to << 733 if ((dxup<=0) && (x+1<y)) {
722 // this would create a hole at (x+1,y) unles << 734 voxels.push_back({ 1+x, y});
723 dxup = dx + 2*y -2*R-1; << 735 voxels.push_back({ y, 1+x});
724 if ((dxup<=0) && (x+1<y)) << 736 voxels.push_back({-1-x, y});
725 { << 737 voxels.push_back({ -y, 1+x});
726 voxels.push_back({ 1+x, y}); << 738 }
727 voxels.push_back({ y, 1+x}); <<
728 voxels.push_back({-1-x, y}); <<
729 voxels.push_back({ -y, 1+x}); <<
730 } <<
731 <<
732 dx += 4-2*y; <<
733 y--; <<
734 } <<
735 <<
736 dx += 1+2*x; <<
737 x++; <<
738 } <<
739 <<
740 // add points on y=+-x <<
741 if (x==y || (x==y+1 && dxup<=0)) <<
742 { <<
743 voxels.push_back({x,x}); <<
744 voxels.push_back({-x,x}); <<
745 } <<
746 739
747 return voxels; << 740 dx += 4-2*y;
>> 741 y--;
>> 742 }
>> 743
>> 744 dx += 1+2*x;
>> 745 x++;
>> 746 }
>> 747
>> 748 // add points on y=+-x
>> 749 if (x==y || (x==y+1 && dxup<=0)) {
>> 750 voxels.push_back({x,x});
>> 751 voxels.push_back({-x,x});
>> 752 }
>> 753
>> 754 return voxels;
748 } 755 }
749 756
750 ////////////////////////////////////////////// 757 ///////////////////////////////////////////////////////////////////////////////
751 // 758 //
752 // Get nearest droplet inside a voxel at (xGri 759 // Get nearest droplet inside a voxel at (xGrid, yGrid, zGrid)
753 // 760 //
754 void FastAerosol::GetNearestDropletInsideGrid( << 761 void FastAerosol::GetNearestDropletInsideGrid(G4double &minDistance, std::vector<G4ThreeVector> &candidates, std::vector<G4double> &distances, unsigned int xGrid, unsigned int yGrid, unsigned int zGrid, const G4ThreeVector &p) {
755 { << 762 unsigned int gi;
756 unsigned int gi; << 763 PopulateGrid(xGrid, yGrid, zGrid, gi);
757 PopulateGrid(xGrid, yGrid, zGrid, gi); <<
758 764
759 // initialize values << 765 // initialize values
760 G4double foundDistance; << 766 G4double foundDistance;
761 // std::vector<G4ThreeVector>::iterator bestP 767 // std::vector<G4ThreeVector>::iterator bestPt;
762 768
763 // find closest droplet << 769 // find closest droplet
764 for (auto it = fGrid[gi].begin(); it != fGrid << 770 for (auto it = fGrid[gi].begin(); it != fGrid[gi].end(); ++it) {
765 { << 771 foundDistance = std::sqrt((p-*it).mag2());
766 foundDistance = std::sqrt((p-*it).mag2()); << 772
767 << 773 if (foundDistance < minDistance+fdR) {
768 if (foundDistance < minDistance+fdR) << 774 minDistance = std::min(minDistance, foundDistance);
769 { << 775
770 minDistance = std::min(minDistance, foundDis << 776 candidates.push_back(*it);
771 candidates.push_back(*it); << 777 distances.push_back(foundDistance);
772 distances.push_back(foundDistance); << 778 }
773 } 779 }
774 } <<
775 } 780 }
776 781
777 ////////////////////////////////////////////// 782 ///////////////////////////////////////////////////////////////////////////////
778 // 783 //
779 // Get nearest droplet along vector 784 // Get nearest droplet along vector
780 // 785 //
781 // Maximum search radius is fVectorSearchRadiu 786 // Maximum search radius is fVectorSearchRadius
782 // Maximum distance travelled along vector is 787 // Maximum distance travelled along vector is maxSearch
783 // 788 //
784 // 1) Find the closest voxel along v that is i 789 // 1) Find the closest voxel along v that is inside the bulk
785 // 2) Search for droplets in a box width (2*fV 790 // 2) Search for droplets in a box width (2*fVectorSearchRadius+1) around this voxel
786 // that our ray would intersect 791 // that our ray would intersect
787 // 3) If so, save the droplet that is closest 792 // 3) If so, save the droplet that is closest along the ray to p and return true
788 // 4) If not, step 0.99*fGridPitch along v unt 793 // 4) If not, step 0.99*fGridPitch along v until in a new voxel
789 // 5) If outside bulk, jump until inside bulk 794 // 5) If outside bulk, jump until inside bulk and then start again from 2
790 // 6) If inside bulk, search the new voxels in 795 // 6) If inside bulk, search the new voxels in a box centered of width (2*fVectorSearchRadius+1)
791 // centered at our new point 796 // centered at our new point
792 // 7) If we find a point, do as in (3). If not 797 // 7) If we find a point, do as in (3). If not, do as in (4)
793 // 798 //
794 // This is searching box-shaped regions of vox 799 // This is searching box-shaped regions of voxels, an approach that is only valid
795 // for droplets which fit inside a voxel, whic 800 // for droplets which fit inside a voxel, which is one reason why the code checks
796 // to make sure that fR < fGridPitch at initia 801 // to make sure that fR < fGridPitch at initialization.
797 // 802 //
798 bool FastAerosol::GetNearestDroplet(const G4Th 803 bool FastAerosol::GetNearestDroplet(const G4ThreeVector &p, const G4ThreeVector &normalizedV, G4ThreeVector ¢er, G4double &minDistance, G4double maxSearch, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation)
799 { 804 {
800 // get the grid box this initial position is << 805 // get the grid box this initial position is inside
801 int xGrid, yGrid, zGrid; << 806 int xGrid, yGrid, zGrid;
802 GetGrid(p, xGrid, yGrid, zGrid); << 807 GetGrid(p, xGrid, yGrid, zGrid);
803 << 808
804 // initialize some values << 809 // initialize some values
805 G4ThreeVector currentP = p; << 810 G4ThreeVector currentP = p;
806 G4double distanceAlongVector = 0.0; << 811 G4double distanceAlongVector = 0.0;
807 G4double deltaDistanceAlongVector; << 812 G4double deltaDistanceAlongVector;
808 << 813
809 int newXGrid = 0; int newYGrid = 0; int newZG << 814 int newXGrid = 0; int newYGrid = 0; int newZGrid = 0;
810 int deltaX = 0; int deltaY = 0; int deltaZ = << 815 int deltaX = 0; int deltaY = 0; int deltaZ = 0;
811 int edgeX = 0 ; int edgeY = 0; int edgeZ = 0; << 816 int edgeX = 0 ; int edgeY = 0; int edgeZ = 0;
812 G4bool boxSearch = true; << 817 G4bool boxSearch = true;
813 818
814 G4double distanceToCloud; << 819 G4double distanceToCloud;
815 820
816 minDistance = kInfinity; << 821 minDistance = kInfinity;
817 822
818 // actual search << 823 // actual search
819 while(true) << 824 while(true) {
820 { << 825 deltaDistanceAlongVector = 0.0;
821 deltaDistanceAlongVector = 0.0; << 826
822 // Jump gaps << 827 // Jump gaps
823 distanceToCloud = DistanceToCloud(currentP, << 828 distanceToCloud = DistanceToCloud(currentP,normalizedV);
824 if (distanceToCloud == kInfinity) << 829 if (distanceToCloud == kInfinity)
825 {return(minDistance != kInfinity);} <<
826 else if (distanceToCloud > fGridPitch) <<
827 { 830 {
828 deltaDistanceAlongVector = distanceToCloud << 831 return(minDistance != kInfinity);
829 distanceAlongVector += deltaDistanceAlongV << 832 }
830 boxSearch = true;// if jumping gap, use bo << 833 else if (distanceToCloud > fGridPitch)
831 currentP += deltaDistanceAlongVector << 834 {
832 // skip gaps << 835 deltaDistanceAlongVector = distanceToCloud-fGridPitch;
>> 836 distanceAlongVector += deltaDistanceAlongVector;
>> 837 boxSearch = true; // if jumping gap, use box search
>> 838 currentP += deltaDistanceAlongVector*normalizedV; // skip gaps
833 } 839 }
834 840
835 // Search for droplets << 841 // Search for droplets
836 if (distanceAlongVector > maxSearch)// quit << 842 if (distanceAlongVector > maxSearch) // quit if will jump too far
837 { <<
838 return(minDistance != kInfinity); <<
839 } <<
840 else <<
841 { <<
842 if (boxSearch) // we jumped <<
843 { <<
844 GetGrid(currentP, xGrid, yGrid, zGrid); <<
845 GetNearestDropletInsideRegion(minDista <<
846 fVectorSearchRadius, fVect <<
847 } <<
848 else // didn't jump <<
849 { 843 {
850 // Searching endcaps << 844 return(minDistance != kInfinity);
851 // ================= << 845 }
852 if (deltaX != 0) << 846 else
853 { 847 {
854 GetNearestDropletInsideRegion(minDistance, << 848 if (boxSearch) // we jumped
855 edgeX, yGrid, zGrid, 0, fVectorSearchRadiu << 849 {
856 droplet, rotation);} << 850 GetGrid(currentP, xGrid, yGrid, zGrid);
857 << 851 GetNearestDropletInsideRegion(minDistance, center,
858 if (deltaY != 0) << 852 xGrid, yGrid, zGrid,
859 { << 853 fVectorSearchRadius, fVectorSearchRadius, fVectorSearchRadius,
860 GetNearestDropletInsideRegion(minDistanc << 854 p, normalizedV,
861 xGrid, edgeY, zGrid, fVectorSearchRadius << 855 droplet, rotation);
862 } << 856 }
863 << 857 else // didn't jump
864 if (deltaZ != 0) << 858 {
865 { << 859 // Searching endcaps
866 GetNearestDropletInsideRegion(minDistanc << 860 // =================
867 } << 861 if (deltaX != 0)
868 << 862 {
869 // Search bars << 863 GetNearestDropletInsideRegion(minDistance, center,
870 // =========== << 864 edgeX, yGrid, zGrid,
871 // (a bar joins two endcaps) << 865 0, fVectorSearchRadius, fVectorSearchRadius,
872 if (deltaX != 0 && deltaY != 0) << 866 p, normalizedV,
873 { GetNearestDropletInsideRegion(minDis << 867 droplet, rotation);
874 0, 0, fVectorSea << 868 }
875 droplet, rotation); << 869
>> 870 if (deltaY != 0)
>> 871 {
>> 872 GetNearestDropletInsideRegion(minDistance, center,
>> 873 xGrid, edgeY, zGrid,
>> 874 fVectorSearchRadius, 0, fVectorSearchRadius,
>> 875 p, normalizedV,
>> 876 droplet, rotation);
876 } 877 }
877 878
878 if (deltaX != 0 && deltaZ != 0) << 879 if (deltaZ != 0)
879 { GetNearestDropletInsideRegion(mi << 880 {
>> 881 GetNearestDropletInsideRegion(minDistance, center,
>> 882 xGrid, yGrid, edgeZ,
>> 883 fVectorSearchRadius, fVectorSearchRadius, 0,
>> 884 p, normalizedV,
>> 885 droplet, rotation);
>> 886 }
880 887
881 if (deltaY != 0 && deltaZ != 0) << 888 // Search bars
882 { << 889 // ===========
883 GetNearestDropletInsideRegion(minDistance, ce << 890 // (a bar joins two endcaps)
884 xGrid, edgeY, e << 891 if (deltaX != 0 && deltaY != 0)
885 fVectorSearchRa << 892 {
886 p, normalizedV, << 893 GetNearestDropletInsideRegion(minDistance, center,
887 } << 894 edgeX, edgeY, zGrid,
888 << 895 0, 0, fVectorSearchRadius,
889 // Search corner << 896 p, normalizedV,
890 // ============= << 897 droplet, rotation);
891 if (deltaX != 0 && deltaY != 0 && deltaZ != 0 << 898 }
892 << 899
893 // Update position << 900 if (deltaX != 0 && deltaZ != 0)
894 // ================================ << 901 {
895 xGrid = newXGrid; yGrid = newYGrid; zGrid = ne << 902 GetNearestDropletInsideRegion(minDistance, center,
896 } << 903 edgeX, yGrid, edgeZ,
>> 904 0, fVectorSearchRadius, 0,
>> 905 p, normalizedV,
>> 906 droplet, rotation);
>> 907 }
>> 908
>> 909 if (deltaY != 0 && deltaZ != 0)
>> 910 {
>> 911 GetNearestDropletInsideRegion(minDistance, center,
>> 912 xGrid, edgeY, edgeZ,
>> 913 fVectorSearchRadius, 0, 0,
>> 914 p, normalizedV,
>> 915 droplet, rotation);
>> 916 }
>> 917
>> 918 // Search corner
>> 919 // =============
>> 920 if (deltaX != 0 && deltaY != 0 && deltaZ != 0)
>> 921 {
>> 922 GetNearestDropletInsideRegion(minDistance, center,
>> 923 edgeX, edgeY, edgeZ,
>> 924 0, 0, 0,
>> 925 p, normalizedV,
>> 926 droplet, rotation);
>> 927 }
>> 928
>> 929 // Update position
>> 930 // ================================
>> 931 xGrid = newXGrid; yGrid = newYGrid; zGrid = newZGrid;
>> 932 }
897 933
898 // Check if we are done << 934 // Check if we are done
899 // ==================== << 935 // ====================
900 if (distanceAlongVector>minDistance+fdR) {retu << 936 if (distanceAlongVector>minDistance+fdR) {
>> 937 return(true);
>> 938 }
901 939
902 // walk along the grid <<
903 // advance by 0.99 fGridPitch (0.99 rather th <<
904 //caused by rounding errors for lines parallel <<
905 <<
906 while (true) { <<
907 deltaDistanceAlongVector = fGridPitch*0.99; <<
908 distanceAlongVector += deltaDistanceAlongVect <<
909 <<
910 // exit returning false if we have traveled m <<
911 if (distanceAlongVector > maxSearch) { return <<
912 <<
913 currentP += deltaDistanceAlongVector*normaliz <<
914 GetGrid(currentP, newXGrid, newYGrid, newZGri <<
915 <<
916 if ((newXGrid != xGrid) || (newYGrid != yGrid <<
917 { <<
918 deltaX = newXGrid - xGrid; <<
919 edgeX = xGrid+deltaX*(1+fVectorSearchRadius) <<
920 deltaY = newYGrid - yGrid; <<
921 edgeY = yGrid+deltaY*(1+fVectorSearchRadius) <<
922 deltaZ = newZGrid - zGrid; <<
923 edgeZ = zGrid+deltaZ*(1+fVectorSearchRadius) <<
924 <<
925 break; <<
926 } <<
927 } <<
928 boxSearch = false; <<
929 } <<
930 } <<
931 940
932 return false; << 941 // walk along the grid
>> 942 // advance by 0.99 fGridPitch (0.99 rather than 1.0 to avoid issues caused by
>> 943 // rounding errors for lines parallel to the grid and based on a grid boundary)
>> 944 while (true) {
>> 945 deltaDistanceAlongVector = fGridPitch*0.99;
>> 946 distanceAlongVector += deltaDistanceAlongVector;
>> 947
>> 948 // exit returning false if we have traveled more than MaxVectorFollowingDistance
>> 949 if (distanceAlongVector > maxSearch) { return(false); }
>> 950
>> 951 currentP += deltaDistanceAlongVector*normalizedV;
>> 952 GetGrid(currentP, newXGrid, newYGrid, newZGrid);
>> 953
>> 954 if ((newXGrid != xGrid) || (newYGrid != yGrid) || (newZGrid != zGrid)) {
>> 955 deltaX = newXGrid - xGrid; edgeX = xGrid+deltaX*(1+fVectorSearchRadius);
>> 956 deltaY = newYGrid - yGrid; edgeY = yGrid+deltaY*(1+fVectorSearchRadius);
>> 957 deltaZ = newZGrid - zGrid; edgeZ = zGrid+deltaZ*(1+fVectorSearchRadius);
>> 958
>> 959 break;
>> 960 }
>> 961 }
>> 962 boxSearch = false;
>> 963 }
>> 964 }
>> 965
>> 966 return false;
933 } 967 }
934 968
935 ////////////////////////////////////////////// 969 ///////////////////////////////////////////////////////////////////////////////
936 // 970 //
937 // Get nearest droplet (along a vector normali 971 // Get nearest droplet (along a vector normalizedV) inside a voxel centered at
938 // (xGridCenter, yGridCenter, zGridCenter) of 972 // (xGridCenter, yGridCenter, zGridCenter) of width (xWidth, yWidth, zWidth)
939 // 973 //
940 // This searches box-shaped regions. 974 // This searches box-shaped regions.
941 // 975 //
942 void FastAerosol::GetNearestDropletInsideRegio << 976 void FastAerosol::GetNearestDropletInsideRegion(G4double &minDistance, G4ThreeVector ¢er, int xGridCenter, int yGridCenter, int zGridCenter, int xWidth, int yWidth, int zWidth, const G4ThreeVector &p, const G4ThreeVector &normalizedV, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation) {
943 { << 977 for (int xGrid = (xGridCenter-xWidth); xGrid<=(xGridCenter+xWidth); xGrid++)
944 for (int xGrid = (xGridCenter-xWidth); xGrid< <<
945 { <<
946 if (0<=xGrid && xGrid<fNx) <<
947 { <<
948 for (int yGrid = (yGridCenter-yWidth); yG <<
949 { 978 {
950 if (0<=yGrid && yGrid<fNy) << 979 if (0<=xGrid && xGrid<fNx)
951 { << 980 {
952 for (int zGrid = (zGridCenter-zWidth); z << 981 for (int yGrid = (yGridCenter-yWidth); yGrid<=(yGridCenter+yWidth); yGrid++)
953 { << 982 {
954 if (0<=zGrid && zGrid<fNz) << 983 if (0<=yGrid && yGrid<fNy)
955 { << 984 {
>> 985 for (int zGrid = (zGridCenter-zWidth); zGrid<=(zGridCenter+zWidth); zGrid++)
>> 986 {
>> 987 if (0<=zGrid && zGrid<fNz)
>> 988 {
956 GetNearestDropletInsideGrid(minD 989 GetNearestDropletInsideGrid(minDistance, center, (unsigned)xGrid, (unsigned)yGrid, (unsigned)zGrid, p, normalizedV, droplet, rotation);
957 } << 990 }
958 } << 991 }
959 } << 992 }
>> 993 }
>> 994 }
960 } 995 }
961 } <<
962 } <<
963 } 996 }
964 997
965 ////////////////////////////////////////////// 998 ///////////////////////////////////////////////////////////////////////////////
966 // 999 //
967 // Get nearest droplet inside a voxel at (xGri 1000 // Get nearest droplet inside a voxel at (xGrid, yGrid, zGrid) along a vector
968 // 1001 //
969 // Return the closest one, as measured along t 1002 // Return the closest one, as measured along the line
970 // 1003 //
971 void FastAerosol::GetNearestDropletInsideGrid( << 1004 void FastAerosol::GetNearestDropletInsideGrid(G4double &minDistance, G4ThreeVector ¢er, unsigned int xGrid, unsigned int yGrid, unsigned int zGrid, const G4ThreeVector &p, const G4ThreeVector &normalizedV, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation) {
972 { << 1005 unsigned int gi;
973 unsigned int gi; << 1006 PopulateGrid(xGrid, yGrid, zGrid, gi);
974 PopulateGrid(xGrid, yGrid, zGrid, gi); <<
975 1007
976 G4double foundDistance; << 1008 G4double foundDistance;
977 1009
978 // find closest droplet << 1010 // find closest droplet
979 for (auto it = fGrid[gi].begin(); it != fGrid << 1011 for (auto it = fGrid[gi].begin(); it != fGrid[gi].end(); ++it) {
980 { << 1012 // could have the following check to see if the ray pierces the bounding sphere. Currently seems like unnecessary addition
981 // could have the following check to see if th << 1013 /*
982 /* << 1014 deltaP = *it-p;
983 deltaP = *it-p; << 1015 foundDistance = normalizedV.dot(deltaP);
984 foundDistance = normalizedV.dot(deltaP); <<
985 <<
986 if ((0<=foundDistance) && ((deltaP - normalize <<
987 */ <<
988 <<
989 G4RotationMatrix irotm = rotation(*it).invers <<
990 <<
991 G4ThreeVector relPos = irotm*(p - *it); <<
992 <<
993 if (droplet->Inside(relPos) == kInside) <<
994 { <<
995 minDistance = 0; <<
996 center = *it; <<
997 } <<
998 else <<
999 { <<
1000 foundDistance = droplet->DistanceToIn( rel <<
1001 1016
1002 if (foundDistance<minDistance) << 1017 if ((0<=foundDistance) && ((deltaP - normalizedV*foundDistance).mag2() < fR2)) {
1003 { << 1018
1004 minDistance = foundDistance; << 1019 }
1005 center = *it; << 1020 */
1006 } << 1021
1007 } << 1022 G4RotationMatrix irotm = rotation(*it).inverse();
1008 } << 1023
>> 1024 G4ThreeVector relPos = irotm*(p - *it);
>> 1025
>> 1026 if (droplet->Inside(relPos) == kInside)
>> 1027 {
>> 1028 minDistance = 0;
>> 1029 center = *it;
>> 1030 }
>> 1031 else
>> 1032 {
>> 1033 foundDistance = droplet->DistanceToIn( relPos , irotm*normalizedV);
>> 1034
>> 1035 if (foundDistance<minDistance)
>> 1036 {
>> 1037 minDistance = foundDistance;
>> 1038 center = *it;
>> 1039 }
>> 1040 }
>> 1041 }
1009 } 1042 }
1010 1043
1011 ///////////////////////////////////////////// 1044 ///////////////////////////////////////////////////////////////////////////////
1012 // 1045 //
1013 // Check if a droplet at the indicated point 1046 // Check if a droplet at the indicated point would collide with an existing droplet
1014 // 1047 //
1015 // Search neighboring voxels, too. Returns tr 1048 // Search neighboring voxels, too. Returns true if there is a collision
1016 // 1049 //
1017 bool FastAerosol::CheckCollision(G4double x, << 1050 bool FastAerosol::CheckCollision(G4double x, G4double y, G4double z) {
1018 { << 1051 G4ThreeVector p(x,y,z);
1019 G4ThreeVector p(x,y,z); <<
1020 1052
1021 std::pair<int, int> minMaxXGrid, minMaxYGrid << 1053 std::pair<int, int> minMaxXGrid, minMaxYGrid, minMaxZGrid;
1022 1054
1023 int xGrid, yGrid, zGrid; << 1055 int xGrid, yGrid, zGrid;
1024 GetGrid(p, xGrid, yGrid, zGrid); << 1056 GetGrid(p, xGrid, yGrid, zGrid);
1025 1057
1026 minMaxXGrid = GetMinMaxSide(xGrid, fNx); << 1058 minMaxXGrid = GetMinMaxSide(xGrid, fNx);
1027 minMaxYGrid = GetMinMaxSide(yGrid, fNy); << 1059 minMaxYGrid = GetMinMaxSide(yGrid, fNy);
1028 minMaxZGrid = GetMinMaxSide(zGrid, fNz); << 1060 minMaxZGrid = GetMinMaxSide(zGrid, fNz);
1029 << 1061
1030 for (int xi=minMaxXGrid.first; xi <= minMaxX << 1062 for (int xi=minMaxXGrid.first; xi <= minMaxXGrid.second; xi++) {
1031 for (int yi = minMaxYGrid.first; yi <= minM << 1063 for (int yi = minMaxYGrid.first; yi <= minMaxYGrid.second; yi++) {
1032 for (int zi = minMaxZGrid.first; zi <= mi << 1064 for (int zi = minMaxZGrid.first; zi <= minMaxZGrid.second; zi++) {
1033 if (CheckCollisionInsideGrid(x, y, z, (unsi << 1065 if (CheckCollisionInsideGrid(x, y, z, (unsigned)xi, (unsigned)yi, (unsigned)zi)) {
1034 fNumCollisions++; // log number of colli << 1066 fNumCollisions++; // log number of collisions for statistics print
1035 return(true); << 1067 return(true);
1036 } << 1068 }
1037 } << 1069 }
1038 } 1070 }
1039 } 1071 }
1040 return(false); 1072 return(false);
1041 } 1073 }
1042 1074
1043 ///////////////////////////////////////////// 1075 ///////////////////////////////////////////////////////////////////////////////
1044 // 1076 //
1045 // Check if a droplet at the indicated point 1077 // Check if a droplet at the indicated point would collide with an existing
1046 // droplet inside a specific grid 1078 // droplet inside a specific grid
1047 // 1079 //
1048 // Note that you don't need to lock the mutex 1080 // Note that you don't need to lock the mutex since this is only called by code
1049 // that already has the mutex (always called 1081 // that already has the mutex (always called by PopulateGrid).
1050 // 1082 //
1051 bool FastAerosol::CheckCollisionInsideGrid(G4 << 1083 bool FastAerosol::CheckCollisionInsideGrid(G4double x, G4double y, G4double z, unsigned int xi, unsigned int yi, unsigned int zi) {
1052 { << 1084 std::vector<G4ThreeVector> *thisGrid = &(fGrid[GetGridIndex(xi, yi, zi)]);
1053 std::vector<G4ThreeVector> *thisGrid = &(fGr << 1085 unsigned int numel = (unsigned int)thisGrid->size();
1054 unsigned int numel = (unsigned int)thisGrid- << 1086
1055 << 1087 for (unsigned int i=0; i < numel; i++) {
1056 for (unsigned int i=0; i < numel; i++) << 1088 if (CheckCollisionWithDroplet(x, y, z, (*thisGrid)[i]))
1057 { << 1089 return(true);
1058 if (CheckCollisionWithDroplet(x, y, z, (*th << 1090 }
1059 return(true); <<
1060 } <<
1061 1091
1062 return(false); << 1092 return(false);
1063 } 1093 }
1064 1094
1065 ///////////////////////////////////////////// 1095 ///////////////////////////////////////////////////////////////////////////////
1066 // 1096 //
1067 // Check for collsion with a specific droplet 1097 // Check for collsion with a specific droplet
1068 // 1098 //
1069 bool FastAerosol::CheckCollisionWithDroplet(G << 1099 bool FastAerosol::CheckCollisionWithDroplet(G4double x, G4double y, G4double z, G4ThreeVector p ) {
1070 { << 1100 return( std::pow(x-p.x(), 2.0) + std::pow(y-p.y(), 2.0) + std::pow(z-p.z(), 2.0) < fCollisionLimit2 );
1071 return( std::pow(x-p.x(), 2.0) + std::pow(y-p <<
1072 } 1101 }
1073 1102
1074 ///////////////////////////////////////////// 1103 ///////////////////////////////////////////////////////////////////////////////
1075 // 1104 //
1076 // Save droplet positions to a file for visua 1105 // Save droplet positions to a file for visualization and analysis
1077 // 1106 //
1078 void FastAerosol::SaveToFile(const char* file << 1107 void FastAerosol::SaveToFile(const char* filename) {
1079 { << 1108 G4cout << "Saving droplet positions to " << filename << "..." << G4endl;
1080 G4cout << "Saving droplet positions to " << << 1109 std::ofstream file;
1081 std::ofstream file; << 1110 file.open(filename);
1082 file.open(filename); << 1111
1083 << 1112 std::vector<G4ThreeVector> voxel;
1084 std::vector<G4ThreeVector> voxel; << 1113 G4ThreeVector pt;
1085 G4ThreeVector pt; << 1114
1086 << 1115 for (auto it1 = fGrid.begin(); it1 != fGrid.end(); ++it1) {
1087 for (auto it1 = fGrid.begin(); it1 != fGrid. << 1116 voxel = *it1;
1088 { << 1117
1089 voxel = *it1; << 1118 for (auto it2 = voxel.begin(); it2 != voxel.end(); ++it2) {
1090 << 1119 pt = *it2;
1091 for (auto it2 = voxel.begin(); it2 != voxe <<
1092 { <<
1093 pt = *it2; <<
1094 1120
1095 double x = pt.x(); << 1121 double x = pt.x();
1096 double y = pt.y(); << 1122 double y = pt.y();
1097 double z = pt.z(); << 1123 double z = pt.z();
1098 file << std::setprecision(15) << x/mm << << 1124
1099 << y/mm << "," << z/mm << "\n"; << 1125 file << std::setprecision(15) << x/mm << "," << y/mm << "," << z/mm << "\n";
1100 } 1126 }
1101 } 1127 }
1102 file.close(); << 1128
>> 1129 file.close();
1103 } 1130 }
1104 1131