fastAerosol/include/FastAerosol.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026
0027 // --------------------------------------------------------------------
0028 // GEANT 4 class header file
0029 //
0030 //
0031 // FastAerosol
0032 //
0033 // Class description:
0034 //
0035 //   A FastAerosol is a collection of points in a voxelized
0036 //   arbitrarily-shaped volume with methods implementing population of
0037 //   grids/voxels and for efficiently finding the nearest point.
0038 //
0039 // Author: A.Knaian (ara@nklabs.com), N.MacFadden (natemacfadden@gmail.com)
0040 // --------------------------------------------------------------------
0041
0042 #ifndef FastAerosol_h
0043 #define FastAerosol_h
0044
0045 #include "globals.hh"
0046 #include "Randomize.hh"
0047 #include "G4ThreeVector.hh"
0048 #include "G4VSolid.hh"
0049
0050 // rotations and number density distribution
0051 #include <functional>
0052 #include "G4RotationMatrix.hh"
0053 #include <atomic>
0054
0055 //using namespace std;
0056
0057 class FastAerosol
0058 {
0059 public:
0060  // Constructor; creates a random cloud of droplets
0061  FastAerosol(const G4String& pName, G4VSolid* pCloud,
0062            G4double pR, G4double pMinD,
0063            G4double pAvgNumDens, G4double pdR,
0064            std::function<G4double (G4ThreeVector)> pNumDensDistribution);
0065
0066 FastAerosol(const G4String& pName, G4VSolid* pCloud,
0067           G4double pR, G4double pMinD,
0068           G4double pNumDens, G4double pdR);
0069
0070 FastAerosol(const G4String& pName, G4VSolid* pCloud,
0071           G4double pR, G4double pMinD, G4double pNumDens);
0072
0073 ~FastAerosol()=default;
0074
0075 // Populate all grids. Otherwise, they are populated on-the-fly
0076 void PopulateAllGrids();
0077
0078 // Save locations of droplets to a file for visualization/analysis purposes
0079 void SaveToFile(const char *filename);
0080
0081 // Get absolutely nearest droplet - must be public as FastAerosolSolid uses it
0082 bool GetNearestDroplet(const G4ThreeVector &p, G4ThreeVector &center, G4double &closestDistance, G4double stepLim, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation);
0083
0084 // Get nearest droplet along a vector - must be public as FastAerosolSolid uses it
0085 bool GetNearestDroplet(const G4ThreeVector &p, const G4ThreeVector &v, G4ThreeVector &center, G4double &closestDistance, G4double stepLim, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation);
0086
0087 // ======
0088 // Inline
0089 // ======
0090 // Input quantities
0091 inline G4String GetName() const; //fasterosol name
0092 inline G4VSolid* GetBulk() const; // bulk shape
0093 inline G4double GetRadius() const; // droplet radius
0094 inline G4double GetAvgNumDens() const;  // droplet number density
0095 //inline G4double GetPitch() const;     // grid pitch
0096
0097 inline G4int GetNumDroplets() const;
0098 // in case the absolute number is more relevant than density
0099
0100 // Bulk quantities
0101 inline G4double GetXHalfLength() const;
0102 inline G4double GetYHalfLength() const;
0103 inline G4double GetZHalfLength() const;
0104 inline void GetBoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const;
0105 inline G4double GetCubicVolume() const;
0106
0107 inline G4double DistanceToCloud(const G4ThreeVector &p);
0108 inline G4double DistanceToCloud(const G4ThreeVector &p, const G4ThreeVector &v);
0109
0110 // Misc getters and setters
0111 inline long GetSeed();
0112 inline void SetSeed(long seed);
0113
0114 inline G4int GetNumPlacementTries();
0115 inline void SetNumPlacementTries(G4int numTries);
0116
0117 inline G4int GetPreSphereR();
0118 inline void SetPreSphereR(G4int fPreSphereRIn);
0119
0120 inline std::function<G4double (G4ThreeVector)> GetDistribution(); // the droplet number density distribution
0121
0122 inline G4double GetDropletsPerVoxel();
0123 inline void SetDropletsPerVoxel(G4double newDropletsPerVoxel);
0124
0125 // Printing diagnostic tool
0126 inline void PrintPopulationReport();
0127
0128 private:
0129  G4double kCarTolerance;
0130
0131 // Parameters, set in constructor
0132  G4String fName;
0133  G4VSolid* fCloud; // Solid volume of the cloud
0134  G4double fDx, fDy, fDz; // Half widths
0135  G4double fR;   // Bounding radius of each droplet
0136  G4double fR2;  // Bounding radius squared of each droplet
0137  G4double fdR;  // Uncertainty in DistanceToIn droplet when just using knowledge of droplet center
0138
0139 G4double fMinD; // Minimum distance allowed between faces of droplets when constructing random array of droplets
0140
0141 std::function<G4double (G4ThreeVector)> fDistribution;
0142 G4double fAvgNumDens;   // Average droplet number density
0143
0144 long int fNumDroplets = 0;
0145 // Number of droplets that have been created
0146
0147 G4double fGridPitch;
0148 // Pitch of collision detection grid.  Must be greater than diameter of droplets for correctness of collision detection.
0149
0150 // Ramdom engine
0151 CLHEP::HepJamesRandom fCloudEngine;
0152 long fSeed = 0; // Global random seed
0153
0154 G4double fDropletsPerVoxel = 4.0;
0155 // Expected number of droplets per voxel
0156
0157 // How far the voxel center must be inside the bulk
0158 //order for there to be no risk of placing a droplet outside
0159 G4double fEdgeDistance;
0160
0161 // Grid variables
0162 std::vector<std::vector<G4ThreeVector>> fGrid;
0163 // Grid of lists of inidices to grid points,
0164 //used for fast collsion checking
0165
0166 std::vector<G4double> fGridMean;
0167 // Array listing mean count for each voxel
0168
0169 std::atomic<bool> *fGridValid;
0170 // Array listing validity of each grid. uses atomic variables
0171
0172 G4int fNx, fNy, fNz; // Number of x, y, and z elements in fGrid
0173
0174 G4int fNxy; // Cached fNx*fNy
0175 long int fNumGridCells; // Cached fNx*fNy*fNz
0176
0177 G4double fCollisionLimit2;
0178 // Threshold distance squared when checking for collsion
0179
0180 G4int fNumNewPointTries = 100;
0181 // How many times we try to place droplets
0182
0183 G4double fMaxDropPercent = 1.0;
0184 // The maximal percentage of skipped droplets before crashing0
0185
0186 G4int fMaxDropCount;
0187 // The maximal number of skipped droplets before crashing
0188
0189 G4int fNumDropped = 0;
0190 // Number of skipped droplets due to collisions/out of bulk placement
0191
0192 G4int fNumCollisions = 0;
0193 // How many collisions occured when attempting to place
0194
0195 // Droplet search variables
0196 G4int fVectorSearchRadius;
0197 // maximum vector search radius
0198
0199 // Droplet placement functions
0200 // ===========================
0201 void InitializeGrid();
0202
0203 G4bool FindNewPoint(G4bool edgeVoxel, G4double dX, G4double dY, G4double dZ, G4double minX, G4double minY, G4double minZ, G4ThreeVector &foundVec);
0204
0205 G4double VoxelOverlap(G4ThreeVector voxelCenter, G4int nStat, G4double epsilon);
0206
0207 bool CheckCollision(G4double x, G4double y, G4double z);
0208 bool CheckCollisionInsideGrid(G4double x, G4double y, G4double z, unsigned int xi, unsigned int yi, unsigned int zi);
0209 bool CheckCollisionWithDroplet(G4double x, G4double y, G4double z, G4ThreeVector p);
0210
0211 // Droplet distance functions
0212 // ==========================
0213 void SearchSphere(G4int searchRad, G4double &minDistance, std::vector<G4ThreeVector> &candidates, std::vector<G4double> &distances2, G4int xGrid, G4int yGrid, G4int zGrid, const G4ThreeVector &p);
0214
0215 void GetNearestDropletInsideRegion(G4double &minDistance, G4ThreeVector &center, int xGrid, int yGrid, int zGrid, int xWidth, int yWidth, int zWidth, const G4ThreeVector &p, const G4ThreeVector &v, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation);
0216
0217 void GetNearestDropletInsideGrid(G4double &minDistance, std::vector<G4ThreeVector> &candidates, std::vector<G4double> &distances2, unsigned int xGrid, unsigned int yGrid, unsigned int zGrid, const G4ThreeVector &p);
0218
0219 void GetNearestDropletInsideGrid(G4double &minDistance, G4ThreeVector &center, unsigned int xGrid, unsigned int yGrid, unsigned int zGrid, const G4ThreeVector &p, const G4ThreeVector &v, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation);
0220
0221 // Voxelized sphere methods
0222 // ========================
0223 // a collection of points as in {{x1,y1},{x2,y2},...}
0224 typedef std::vector<std::vector<int>> fCircleType;
0225
0226 // a collection of points describing a spherical shell
0227 // with points (x,y,z)=(i-R,j-R,sphere[i][j][k])
0228 // that is, first index gives x-position, second index
0229 // gives y-position, and the value gives z-position
0230 //
0231 // this is done so that searching may be optimized:
0232 // if searching some x=i-R that is outside the
0233 // aerosol's bounding box, immediately increment i
0234 // (similar for y).
0235 typedef std::vector<std::vector<std::vector<int>>> fSphereType;
0236
0237 G4int fMaxCircleR;
0238 G4int fMaxSphereR;
0239 G4int fPreSphereR = 20;
0240 std::vector<fCircleType> fCircleCollection;
0241 std::vector<fSphereType> fSphereCollection;
0242 fSphereType MakeSphere(G4int R);
0243 fCircleType MakeCircle(G4int R);
0244 fCircleType MakeHalfCircle(G4int R);
0245
0246 void PopulateGrid(unsigned int xi, unsigned int yi, unsigned int zi, unsigned int& gi);
0247
0248 // ======
0249 // Inline
0250 // ======
0251 inline bool GetGrid(const G4ThreeVector &p, G4int &xGrid, G4int &yGrid, G4int &zGrid);
0252
0253 inline bool AnyIndexOutOfBounds(G4int xGrid, G4int yGrid, G4int zGrid);
0254
0255 inline unsigned int GetGridIndex(unsigned int xi, unsigned int yi, unsigned int zi);
0256
0257 inline G4ThreeVector GetIndexCoord(G4int index);
0258
0259 inline std::pair<G4int, G4int> GetMinMaxSide(G4int index, G4int numGrids);
0260 };
0261
0262 #include "FastAerosol.icc"
0263
0264 #endif