| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26
27 // -------------------------------------------
28 // Implementation for FastAerosol class
29 // Author: A.Knaian (ara@nklabs.com), N.MacFad
30 // -------------------------------------------
31
32 #include "FastAerosol.hh"
33
34 #include "G4SystemOfUnits.hh"
35 #include "G4GeometryTolerance.hh"
36
37 #include <numeric> // for summing vectors wit
38
39 // multithreaded safety
40 //#include <atomic>
41 #include "G4AutoLock.hh"
42 namespace
43 {
44 G4Mutex gridMutex = G4MUTEX_INITIALIZER;
45 G4Mutex sphereMutex = G4MUTEX_INITIALIZER;
46 }
47
48 //////////////////////////////////////////////
49 //
50 // Constructor - check inputs
51 // - initialize grid
52 // - cache values
53 //
54 FastAerosol::FastAerosol(const G4String& pName
55 G4double pR, G4double pMinD, G4do
56 G4double pdR,
57 std::function<G4double (G4ThreeV
58 : fName(pName), fCloud(pCloud), fMinD(pMinD)
59 {
60 kCarTolerance = G4GeometryTolerance::GetInsta
61
62 // get and set bounding box dimensions
63 G4ThreeVector minBounding, maxBounding;
64 fCloud->BoundingLimits(minBounding, maxBoundi
65 G4ThreeVector halfSizeVec = 0.5*(maxBounding
66
67 G4double pX = halfSizeVec[0];
68 G4double pY = halfSizeVec[1];
69 G4double pZ = halfSizeVec[2];
70
71 if (pX < 2*kCarTolerance ||
72 pY < 2*kCarTolerance ||
73 pZ < 2*kCarTolerance) // limit to thicknes
74 {
75 std::ostringstream message;
76 message << "Dimensions too small for cloud:
77 << GetName() << "!" << G4endl
78 << " fDx, fDy, fDz = " << pX << ", "
79 G4Exception("FastAerosol::FastAerosol()",
80 FatalException, message);
81 }
82 fDx = pX;
83 fDy = pY;
84 fDz = pZ;
85
86
87 // check and set droplet radius
88 if (pR < 0.0)
89 {
90 std::ostringstream message;
91 message << "Invalid radius for cloud: " <<
92 << " Radius must be positive."
93 << " Inputs: pR = " << pR;
94 G4Exception("FastAerosol::FastAerosol()",
95 FatalErrorInArgument, message);
96 }
97 fR = pR;
98 fR2 = fR*fR;
99
100
101 // check and set droplet radius safety
102 if (pdR < 0.0 || pdR > fR)
103 {
104 std::ostringstream message;
105 message << "Invalid sphericality measure f
106 << " Radius uncertainty must be
107 << " Inputs: pdR = " << pdR
108 << " Inputs: pR = " << pR;
109 G4Exception("FastAerosol::FastAerosol()",
110 FatalErrorInArgument, message);
111 }
112 fdR = pdR;
113
114
115 // check and set number density
116 if (pAvgNumDens <= 0)
117 {
118 std::ostringstream message;
119 message << "Invalid average number density
120 << " pAvgNumDens = " << pAvgNumDen
121 G4Exception("FastAerosol::FastAerosol()",
122 FatalException, message);
123 }
124 fAvgNumDens = pAvgNumDens;
125
126
127 // Set collision limit for collsion between
128 // no droplets will be placed closer than th
129 fCollisionLimit2 = (2*fR + fMinD)*(2*fR + fM
130
131 // set maximum number of droplet that we are
132 fMaxDropCount = (G4int)floor(fAvgNumDens*(fC
133
134 // initialize grid variables
135 InitializeGrid();
136
137
138 // begin cache of voxelized circles and sphe
139 // see header for more details on these data
140 G4AutoLock lockSphere(&sphereMutex); // loc
141
142 fCircleCollection.push_back({{0, 0}}); // t
143 fSphereCollection.push_back({{{0}}}); // the
144
145 fMaxCircleR = 0;
146 fMaxSphereR = 0;
147
148 MakeSphere(fPreSphereR);
149
150 lockSphere.unlock(); // unlock
151
152 // vector search radius. In terms of voxel w
153 // you need to search a larger area if fR is
154 fVectorSearchRadius = (G4int)ceill(fR/fGridP
155 }
156
157 //////////////////////////////////////////////
158 //
159 // Alternative Constructor
160 //
161 // Same as standard constructor except with a
162 //
163 FastAerosol::FastAerosol(const G4String& pName
164 FastAerosol(pName, pCloud, pR, pMinD, pNumDe
165 [](G4ThreeVector) {return 1.0;})
166 {}
167
168 //////////////////////////////////////////////
169 //
170 // Alternative Constructor
171 //
172 // Same as standard constructor except with a
173 // assuming no uncertainty in sphericality
174 //
175 FastAerosol::FastAerosol(const G4String& pName
176 {}
177
178 //////////////////////////////////////////////
179 //
180 // Initialize grid
181 //
182 // Sets grids to initial values and calculates
183 // for each voxel
184 //
185 void FastAerosol::InitializeGrid()
186 {
187 // set pitch so, on average, fDropletsPerVoxe
188 fGridPitch = std::pow(fDropletsPerVoxel/fAvgN
189
190 // if a voxel has center farther than this di
191 // we know it is fully contained in the bulk
192 fEdgeDistance = fGridPitch*std::sqrt(3.0)/2.0
193
194 // set number of grid cells
195 fNx = (G4int)ceill(2*fDx / fGridPitch);
196 fNy = (G4int)ceill(2*fDy / fGridPitch);
197 fNz = (G4int)ceill(2*fDz / fGridPitch);
198 fNumGridCells = (long) fNx*fNy*fNz;
199 fNxy = fNx*fNy;
200
201 // try... catch... because vectors can only b
202 // thus you can't set grid too fine
203 try
204 {
205 std::vector<G4ThreeVector> emptyVoxel{};
206 fGrid.resize(fNumGridCells, emptyVoxel);
207 fGridMean.resize(fNumGridCells, 0);
208 fGridValid = new std::atomic<bool>[fNumGrid
209 }
210 catch ( const std::bad_alloc&)
211 {
212 std::ostringstream message;
213 message << "Out of memory! Grid pitch too s
214 << GetName() << "!" << G4endl
215 << " Asked for fNumGridCells = "
216 << fNumGridCells << G4endl
217 << " each with element of size "
218 << sizeof(std::vector<G4ThreeVector>) <
219 << " each with element of size " <
220 << G4endl
221 << " but the max size is " << fGri
222 G4Exception("FastAerosol::FastAerosol()", "
223 FatalErrorInArgument, message); }
224
225 // initialize grid cache
226 G4ThreeVector voxelCenter;
227 G4bool valid;
228
229 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 {
255 volScaling = VoxelOverlap(voxelCenter, 100,
256 }
257 // calculates number of droplets based off
258 fGridMean[i] = std::max(0.0, volScaling*fDis
259 }
260 }
261
262 // must scale fGridMean[i] so that the desired
263 // this is because no restrictions are applied
264 G4double tempScaling = (fCloud->GetCubicVolume
265 for (int i=0; i<fNumGridCells; i++) {fGridMe
266 }
267
268 //////////////////////////////////////////////
269 //
270 // Calculate the overlap of the voxel with the
271 //
272 // The method is largely copied from G4VSolid:
273 //
274 G4double FastAerosol::VoxelOverlap(G4ThreeVect
275 {
276 G4double cenX = voxelCenter.x();
277 G4double cenY = voxelCenter.y();
278 G4double cenZ = voxelCenter.z();
279
280 G4int iInside=0;
281 G4double px,py,pz;
282 G4ThreeVector p;
283 G4bool in;
284
285 if(nStat < 100) nStat = 100;
286 if(epsilon > 0.01) epsilon = 0.01;
287
288 for(G4int i = 0; i < nStat; i++ )
289 {
290 px = cenX+(fGridPitch-epsilon)*(G4UniformRan
291 py = cenY+(fGridPitch-epsilon)*(G4UniformRan
292 pz = cenZ+(fGridPitch-epsilon)*(G4UniformRan
293 p = G4ThreeVector(px,py,pz);
294 in = (fCloud->Inside(p) == kInside) && (fClo
295 if(in) iInside++;
296 }
297
298 return (double)iInside/nStat;
299 }
300
301 //////////////////////////////////////////////
302 //
303 // Populate all voxels
304 //
305 // Allows for partially populated clouds.
306 //
307 void FastAerosol::PopulateAllGrids()
308 {
309 unsigned int gi;
310
311 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 {
317 // PopulateGrid takes unsigned arguments
318 // fNx, fNy, and fNz are signed so that s
319 // thus we iterate with xi, yi, and zi sig
320 PopulateGrid((unsigned)xi,(unsigned)yi,(un
321 }
322 }
323 }
324 }
325
326 //////////////////////////////////////////////
327 //
328 // Populate a single grid
329 //
330 // If grid is already populated, does nothing.
331 //
332 void FastAerosol::PopulateGrid(unsigned int xi
333 {
334 gi = GetGridIndex(xi,yi,zi);
335
336 // Check if this grid needs update
337 bool tmpValid = fGridValid[gi].load(std::memo
338 // have to do this weirdly because we are usi
339 std::atomic_thread_fence(std::memory_order_ac
340
341 if (!tmpValid) // if true then either outside
342 {
343 G4AutoLock lockGrid(&gridMutex);
344 // Check again now that we have the lock - in
345 tmpValid = fGridValid[gi].load(std::memory_o
346
347 if (!tmpValid)
348 {
349 // uniquely set the seed to randomly place dro
350 // changing global seed gaurantees a totally n
351 fCloudEngine.setSeed((long)gi + (long)fNum
352
353 // find if the voxel is near the bulk edge. In
354 // if not an edge voxel, we can place droplets
355 G4ThreeVector voxelCenter = fGridPitch*G4Thr
356 G4bool edgeVoxel = ( kInside != fCloud->Insid
357
358 // number of droplets to place
359 unsigned int numDropletsToPlace = CLHEP::Ran
360
361 // actually add the points
362 G4ThreeVector point;
363
364 while (numDropletsToPlace > 0)
365 {
366 // find a new point not overlapping with th
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
379
380 fGridValid[gi].store(true, std::memory_order_r
381 }
382
383 lockGrid.unlock();
384 }
385 }
386
387 //////////////////////////////////////////////
388 //
389 // Find a new droplet position in a box of ful
390 // with minimum corner at (minX, minY, minZ).
391 //
392 // -------------------------------------------
393 //
394 // Ex: FindNewPoint(fGridPitch, fGridPitch, fG
395 // (G4double)xi*fGridPitch-f
396 // (G4double)zi*fGridPitch-f
397 //
398 // 1) CLHEP::RandFlat::shoot(&fCloudEngine)*fG
399 // [0, fGridPitch]
400 //
401 // 1.5) Note that the minimum x value of the (
402 // xi*fGridPitch-fDx
403 //
404 // 2) xi*fGridPitch-fDx + (1) adds the minimum
405 // voxel
406 //
407 // -------------------------------------------
408 //
409 // Ex: FindNewPoint(2.0*fDx, 2.0*fDy, 2.0*fDz,
410 //
411 // 1) CLHEP::RandFlat::shoot(&fCloudEngine)*2.
412 // range [0, 2.0*fDx]
413 //
414 // 2) add -fDx to change range into [-fDx, fDx
415 //
416 G4bool FastAerosol::FindNewPoint(G4bool edgeVo
417 {
418 G4int tries = 0;
419 // counter of tries. Give up after fNumNewPoi
420
421 G4double x, y, z;
422
423 G4ThreeVector point;
424 G4bool placedOutside = false;
425 // variable whether or not the point was plac
426
427 // Generate a droplet and check if it overlap
428 do {
429 tries++;
430 if (tries > fNumNewPointTries)
431 // skip if we tried more than fNumNewPoi
432 {
433 fNumDropped++;
434 if (fNumDropped < fMaxDropCount)
435 // return error if we skipped more than fMa
436 {
437 return false;
438 }
439
440 std::ostringstream message;
441 message << "Threw out too many droplets for
442 << GetName() << G4endl
443 << " Tried to place individual drop
444 << fNumNewPointTries << " times." << G4end
445 << " This failed for " << fMaxDropC
446 << " droplets.";
447 G4Exception("FastAerosol::FindNewPoint()",
448 FatalErrorInArgument, message);
449 }
450
451 x = minX + CLHEP::RandFlat::shoot(&fCloudE
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 }
460 }
461 while (CheckCollision(x,y,z) || placedOutsid
462
463 foundPt = G4ThreeVector(x,y,z);
464 return true;
465 }
466
467 //////////////////////////////////////////////
468 //
469 // Get absolutely nearest droplet
470 //
471 // Maximum search radius is maxSearch
472 //
473 // Start in grid of point p, check if any sper
474 // and then go out one more (spherically-shape
475 // and so on, until the whole volume has been
476 // in general as you go out, you always need t
477 // the one you have (at the closer level) is t
478 //
479 bool FastAerosol::GetNearestDroplet(const G4Th
480 {
481 G4double cloudDistance = fCloud->DistanceToIn
482
483 // starting radius/diam of voxel layer
484 G4int searchRad = (G4int)floor(0.5+cloudDista
485
486 // find the voxel containing p
487 int xGrid, yGrid, zGrid;
488 GetGrid(p, xGrid, yGrid, zGrid);
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
505 SearchSphere(searchRad, minDistance, candid
506 maxCheckDistance = minDistance+fdR;
507
508 // theory says that, to safely have searched
509 // *** unsure if fully accurate. Calculati
510 unsafe = searchRad < std::ceil(0.25+maxCheckD
511
512 searchRad++;
513 }
514
515 // delete any collected droplets that are too
516 unsigned int index = 0;
517
518 G4ThreeVector tempCenter;
519 G4double tempDistance;
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
544 return true;
545 }
546
547 //////////////////////////////////////////////
548 //
549 // Search sphere of radius R for droplets
550 //
551 // Saves any droplets to "center" if they are
552 // (updating this minimum found distance at ea
553 //
554 // Returns if minDistance<infinity (i.e., if w
555 //
556 void FastAerosol::SearchSphere(G4int searchRad
557 {
558 // search variables
559 G4int xSearch, ySearch, zSearch;
560 // x, y, and z coordinates of the currently s
561 // in the shell voxel layer variables
562 fSphereType shell; // the shell that we are s
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
597 if (0<=ySearch && ySearch<fNy)
598 {
599 for (auto it = shell[i][j].begin(); it !=
600 zSearch = *it+zGrid;
601
602 if (0<=zSearch && zSearch<fNz)
603 { GetNearestDropletInsideGrid(
604 }
605 }
606 }
607 }
608 }
609 }
610 }
611
612 //////////////////////////////////////////////
613 //
614 // Make voxelized spheres up to radius R for f
615 //
616 // These spheres are used for searching for dr
617 //
618 // To create them, we first create a helper ha
619 // voxel/point in zr represents a layer of our
620 // zr represents the radius of the layer and t
621 // layer's displacement from the origin along
622 // z in our convention).
623 //
624 std::vector<std::vector<std::vector<int>>> Fas
625 // inductively make smaller spheres since fSph
626 if (fMaxSphereR<(R-1)) { MakeSphere(R-1);}
627
628 std::vector<int> z;
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
646 for (auto it2 = circ.begin(); it2 != circ.en
647 {
648 std::vector<int> pt2 = *it2;
649 x[pt2[0]+R][pt2[1]+R].push_back(pt1[0]);
650 }
651 }
652
653 fSphereCollection.push_back(x);
654 fMaxSphereR = R;
655 return x;
656 }
657
658 //////////////////////////////////////////////
659 //
660 // Make voxelized circles up to radius R for f
661 //
662 // These circles are used to create voxelized
663 // used for searching for droplets). This just
664 // a half-circle, adding the points (R,0) and
665 // circle along the x-axis.
666 //
667 std::vector<std::vector<int>> FastAerosol::Mak
668 {
669 // inductively make smaller circles since fCir
670 if (fMaxCircleR<R)
671 {
672 if (fMaxCircleR<(R-1)) {MakeCircle(R-1);}
673 fCircleType voxels = {{R, 0}, {-R, 0}};
674 // add (R,0) and (-R,0) since half circle e
675 fCircleType halfVoxels = MakeHalfCircle(R);
676
677 // join the voxels, halfVoxels, and -halfVo
678 for (auto it = halfVoxels.begin(); it != ha
679 {
680 std::vector<int> pt = *it;
681 voxels.push_back(pt);
682 voxels.push_back({-pt[0], -pt[1]});
683 }
684 fCircleCollection.push_back(voxels);
685 fMaxCircleR++;
686 }
687
688 return fCircleCollection[R];
689 }
690
691 //////////////////////////////////////////////
692 //
693 // Make half circle by a modified midpoint cir
694 //
695 // Modifies the algorithm so it doesn't have '
696 //
697 // See B. Roget, J. Sitaraman, "Wall distance
698 // voxelized marching spheres," Journal of Com
699 // May 2013, pp. 76-94, https://doi.org/10.101
700 //
701 std::vector<std::vector<int>> FastAerosol::Mak
702 {
703 // makes an octant of a voxelized circle in th
704 fCircleType voxels = {{0, R}}; // hard code i
705 int x = 1;
706 int y = R;
707 int dx = 4-R;
708 // measure of whether (x+1,y) will be outside
709 int dxup = 1;
710 // measure of whether (x+1,y+1) will be outsi
711 while (x<y)
712 {
713 // mirror the points to be added to change o
714 voxels.push_back({ x, y});
715 voxels.push_back({ y, x});
716 voxels.push_back({-x, y});
717 voxels.push_back({-y, x});
718 // if next point outside the sphere, de-incr
719 if (dx>0)
720 {
721 // if dxup<=0, the layer above just moves to
722 // this would create a hole at (x+1,y) unles
723 dxup = dx + 2*y -2*R-1;
724 if ((dxup<=0) && (x+1<y))
725 {
726 voxels.push_back({ 1+x, y});
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
747 return voxels;
748 }
749
750 //////////////////////////////////////////////
751 //
752 // Get nearest droplet inside a voxel at (xGri
753 //
754 void FastAerosol::GetNearestDropletInsideGrid(
755 {
756 unsigned int gi;
757 PopulateGrid(xGrid, yGrid, zGrid, gi);
758
759 // initialize values
760 G4double foundDistance;
761 // std::vector<G4ThreeVector>::iterator bestP
762
763 // find closest droplet
764 for (auto it = fGrid[gi].begin(); it != fGrid
765 {
766 foundDistance = std::sqrt((p-*it).mag2());
767
768 if (foundDistance < minDistance+fdR)
769 {
770 minDistance = std::min(minDistance, foundDis
771 candidates.push_back(*it);
772 distances.push_back(foundDistance);
773 }
774 }
775 }
776
777 //////////////////////////////////////////////
778 //
779 // Get nearest droplet along vector
780 //
781 // Maximum search radius is fVectorSearchRadiu
782 // Maximum distance travelled along vector is
783 //
784 // 1) Find the closest voxel along v that is i
785 // 2) Search for droplets in a box width (2*fV
786 // that our ray would intersect
787 // 3) If so, save the droplet that is closest
788 // 4) If not, step 0.99*fGridPitch along v unt
789 // 5) If outside bulk, jump until inside bulk
790 // 6) If inside bulk, search the new voxels in
791 // centered at our new point
792 // 7) If we find a point, do as in (3). If not
793 //
794 // This is searching box-shaped regions of vox
795 // for droplets which fit inside a voxel, whic
796 // to make sure that fR < fGridPitch at initia
797 //
798 bool FastAerosol::GetNearestDroplet(const G4Th
799 {
800 // get the grid box this initial position is
801 int xGrid, yGrid, zGrid;
802 GetGrid(p, xGrid, yGrid, zGrid);
803
804 // initialize some values
805 G4ThreeVector currentP = p;
806 G4double distanceAlongVector = 0.0;
807 G4double deltaDistanceAlongVector;
808
809 int newXGrid = 0; int newYGrid = 0; int newZG
810 int deltaX = 0; int deltaY = 0; int deltaZ =
811 int edgeX = 0 ; int edgeY = 0; int edgeZ = 0;
812 G4bool boxSearch = true;
813
814 G4double distanceToCloud;
815
816 minDistance = kInfinity;
817
818 // actual search
819 while(true)
820 {
821 deltaDistanceAlongVector = 0.0;
822 // Jump gaps
823 distanceToCloud = DistanceToCloud(currentP,
824 if (distanceToCloud == kInfinity)
825 {return(minDistance != kInfinity);}
826 else if (distanceToCloud > fGridPitch)
827 {
828 deltaDistanceAlongVector = distanceToCloud
829 distanceAlongVector += deltaDistanceAlongV
830 boxSearch = true;// if jumping gap, use bo
831 currentP += deltaDistanceAlongVector
832 // skip gaps
833 }
834
835 // Search for droplets
836 if (distanceAlongVector > maxSearch)// quit
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 {
850 // Searching endcaps
851 // =================
852 if (deltaX != 0)
853 {
854 GetNearestDropletInsideRegion(minDistance,
855 edgeX, yGrid, zGrid, 0, fVectorSearchRadiu
856 droplet, rotation);}
857
858 if (deltaY != 0)
859 {
860 GetNearestDropletInsideRegion(minDistanc
861 xGrid, edgeY, zGrid, fVectorSearchRadius
862 }
863
864 if (deltaZ != 0)
865 {
866 GetNearestDropletInsideRegion(minDistanc
867 }
868
869 // Search bars
870 // ===========
871 // (a bar joins two endcaps)
872 if (deltaX != 0 && deltaY != 0)
873 { GetNearestDropletInsideRegion(minDis
874 0, 0, fVectorSea
875 droplet, rotation);
876 }
877
878 if (deltaX != 0 && deltaZ != 0)
879 { GetNearestDropletInsideRegion(mi
880
881 if (deltaY != 0 && deltaZ != 0)
882 {
883 GetNearestDropletInsideRegion(minDistance, ce
884 xGrid, edgeY, e
885 fVectorSearchRa
886 p, normalizedV,
887 }
888
889 // Search corner
890 // =============
891 if (deltaX != 0 && deltaY != 0 && deltaZ != 0
892
893 // Update position
894 // ================================
895 xGrid = newXGrid; yGrid = newYGrid; zGrid = ne
896 }
897
898 // Check if we are done
899 // ====================
900 if (distanceAlongVector>minDistance+fdR) {retu
901
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
932 return false;
933 }
934
935 //////////////////////////////////////////////
936 //
937 // Get nearest droplet (along a vector normali
938 // (xGridCenter, yGridCenter, zGridCenter) of
939 //
940 // This searches box-shaped regions.
941 //
942 void FastAerosol::GetNearestDropletInsideRegio
943 {
944 for (int xGrid = (xGridCenter-xWidth); xGrid<
945 {
946 if (0<=xGrid && xGrid<fNx)
947 {
948 for (int yGrid = (yGridCenter-yWidth); yG
949 {
950 if (0<=yGrid && yGrid<fNy)
951 {
952 for (int zGrid = (zGridCenter-zWidth); z
953 {
954 if (0<=zGrid && zGrid<fNz)
955 {
956 GetNearestDropletInsideGrid(minD
957 }
958 }
959 }
960 }
961 }
962 }
963 }
964
965 //////////////////////////////////////////////
966 //
967 // Get nearest droplet inside a voxel at (xGri
968 //
969 // Return the closest one, as measured along t
970 //
971 void FastAerosol::GetNearestDropletInsideGrid(
972 {
973 unsigned int gi;
974 PopulateGrid(xGrid, yGrid, zGrid, gi);
975
976 G4double foundDistance;
977
978 // find closest droplet
979 for (auto it = fGrid[gi].begin(); it != fGrid
980 {
981 // could have the following check to see if th
982 /*
983 deltaP = *it-p;
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
1002 if (foundDistance<minDistance)
1003 {
1004 minDistance = foundDistance;
1005 center = *it;
1006 }
1007 }
1008 }
1009 }
1010
1011 /////////////////////////////////////////////
1012 //
1013 // Check if a droplet at the indicated point
1014 //
1015 // Search neighboring voxels, too. Returns tr
1016 //
1017 bool FastAerosol::CheckCollision(G4double x,
1018 {
1019 G4ThreeVector p(x,y,z);
1020
1021 std::pair<int, int> minMaxXGrid, minMaxYGrid
1022
1023 int xGrid, yGrid, zGrid;
1024 GetGrid(p, xGrid, yGrid, zGrid);
1025
1026 minMaxXGrid = GetMinMaxSide(xGrid, fNx);
1027 minMaxYGrid = GetMinMaxSide(yGrid, fNy);
1028 minMaxZGrid = GetMinMaxSide(zGrid, fNz);
1029
1030 for (int xi=minMaxXGrid.first; xi <= minMaxX
1031 for (int yi = minMaxYGrid.first; yi <= minM
1032 for (int zi = minMaxZGrid.first; zi <= mi
1033 if (CheckCollisionInsideGrid(x, y, z, (unsi
1034 fNumCollisions++; // log number of colli
1035 return(true);
1036 }
1037 }
1038 }
1039 }
1040 return(false);
1041 }
1042
1043 /////////////////////////////////////////////
1044 //
1045 // Check if a droplet at the indicated point
1046 // droplet inside a specific grid
1047 //
1048 // Note that you don't need to lock the mutex
1049 // that already has the mutex (always called
1050 //
1051 bool FastAerosol::CheckCollisionInsideGrid(G4
1052 {
1053 std::vector<G4ThreeVector> *thisGrid = &(fGr
1054 unsigned int numel = (unsigned int)thisGrid-
1055
1056 for (unsigned int i=0; i < numel; i++)
1057 {
1058 if (CheckCollisionWithDroplet(x, y, z, (*th
1059 return(true);
1060 }
1061
1062 return(false);
1063 }
1064
1065 /////////////////////////////////////////////
1066 //
1067 // Check for collsion with a specific droplet
1068 //
1069 bool FastAerosol::CheckCollisionWithDroplet(G
1070 {
1071 return( std::pow(x-p.x(), 2.0) + std::pow(y-p
1072 }
1073
1074 /////////////////////////////////////////////
1075 //
1076 // Save droplet positions to a file for visua
1077 //
1078 void FastAerosol::SaveToFile(const char* file
1079 {
1080 G4cout << "Saving droplet positions to " <<
1081 std::ofstream file;
1082 file.open(filename);
1083
1084 std::vector<G4ThreeVector> voxel;
1085 G4ThreeVector pt;
1086
1087 for (auto it1 = fGrid.begin(); it1 != fGrid.
1088 {
1089 voxel = *it1;
1090
1091 for (auto it2 = voxel.begin(); it2 != voxe
1092 {
1093 pt = *it2;
1094
1095 double x = pt.x();
1096 double y = pt.y();
1097 double z = pt.z();
1098 file << std::setprecision(15) << x/mm <<
1099 << y/mm << "," << z/mm << "\n";
1100 }
1101 }
1102 file.close();
1103 }
1104