fastAerosol/include/FastAerosol.icc

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 // GEANT 4 inline definitions file
0028 //
0029 // FastAerosol.icc
0030 //
0031 // Implementation of inline methods of FastAerosol
0032 //
0033 // Author: A.Knaian (ara@nklabs.com), N.MacFadden (natemacfadden@gmail.com)
0034 // --------------------------------------------------------------------
0035
0036 inline
0037 G4String FastAerosol::GetName() const
0038 {
0039     return(fName);
0040 }
0041
0042 inline
0043 G4VSolid* FastAerosol::GetBulk() const
0044 {
0045     return(fCloud);
0046 }
0047
0048 inline
0049 G4double FastAerosol::GetRadius() const
0050 {
0051     return(fR);
0052 }
0053
0054 inline
0055 G4double FastAerosol::GetAvgNumDens() const
0056 {
0057     return(fAvgNumDens);
0058 }
0059
0060 inline
0061 G4int FastAerosol::GetNumDroplets() const
0062 {
0063     return((int)(fAvgNumDens*GetCubicVolume()));
0064 }
0065
0066
0067
0068 inline
0069 G4double FastAerosol::GetXHalfLength() const
0070 {
0071   return(fDx);
0072 }
0073
0074 inline
0075 G4double FastAerosol::GetYHalfLength() const
0076 {
0077   return(fDy);
0078 }
0079
0080 inline
0081 G4double FastAerosol::GetZHalfLength() const
0082 {
0083   return(fDz);
0084 }
0085
0086 inline
0087 void FastAerosol::GetBoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
0088 {
0089     pMin.setX(-fDx); pMax.setX(fDx);
0090     pMin.setY(-fDy); pMax.setY(fDy);
0091     pMin.setZ(-fDz); pMax.setZ(fDz);
0092 }
0093
0094 inline
0095 G4double FastAerosol::GetCubicVolume() const
0096 {
0097   return(fCloud->GetCubicVolume());
0098 }
0099
0100
0101 // Find the absolute distance to the cloud bulk from p
0102 inline
0103 G4double FastAerosol::DistanceToCloud(const G4ThreeVector &p) {
0104     if (fCloud->Inside(p)==kOutside)
0105     {
0106         return(fCloud->DistanceToIn(p));
0107     }
0108     else
0109     {
0110         return(0);
0111     }
0112 }
0113
0114 // Find the distance to the cloud bulk from p along a vector v
0115 inline
0116 G4double FastAerosol::DistanceToCloud(const G4ThreeVector &p, const G4ThreeVector &v) {
0117     if (fCloud->Inside(p)==kInside)
0118     {
0119         return(0);
0120     }
0121     else
0122     {
0123         return(fCloud->DistanceToIn(p,v));
0124     }
0125 }
0126
0127
0128 // Get and set the base of the random seed used to set droplet positions
0129 // For a given seed and a given geometry, the droplet positions are the same
0130 inline
0131 long FastAerosol::GetSeed() {
0132     return(fSeed);
0133 }
0134
0135 inline
0136 void FastAerosol::SetSeed(long seedIn) {
0137     fSeed = seedIn;
0138 }
0139
0140
0141 // Get and set the maximum radius of the voxelized spheres to pre-populate
0142 inline
0143 G4int FastAerosol::GetPreSphereR() {
0144     return(fPreSphereR);
0145 }
0146
0147 inline
0148 void FastAerosol::SetPreSphereR(G4int preSphereRIn) {
0149     fPreSphereR = preSphereRIn;
0150 }
0151
0152
0153 // Get the droplet distribution function
0154 inline
0155 std::function<G4double (G4ThreeVector)> FastAerosol::GetDistribution() {
0156     return(fDistribution);
0157 }
0158
0159
0160 // Get and set the maximum number of droplet placement tries in the solid
0161 // Skip placement if attempting to place a single droplet more than this
0162 inline
0163 G4int FastAerosol::GetNumPlacementTries()
0164 {
0165     return(fNumNewPointTries);
0166 }
0167
0168 inline
0169 void FastAerosol::SetNumPlacementTries(G4int numTries)
0170 {
0171     fNumNewPointTries = numTries;
0172 }
0173
0174 // Get and set the expected number of droplets per voxel (on average)
0175 inline
0176 G4double FastAerosol::GetDropletsPerVoxel()
0177 {
0178   return(fDropletsPerVoxel);
0179 }
0180
0181 inline
0182 void FastAerosol::SetDropletsPerVoxel(G4double newDropletsPerVoxel)
0183 {
0184   if (newDropletsPerVoxel >= std::pow(4.0*std::sqrt(2),-1.0))
0185   {
0186     fDropletsPerVoxel = newDropletsPerVoxel;
0187     InitializeGrid();
0188   }
0189   else
0190   {
0191     std::ostringstream message;
0192         message << "Invalid droplets/voxel for cloud: " << GetName() << "!" << G4endl
0193                 << "        For grid pitch to be larger than radius (currently assumed),"
0194                 << "        droplets/voxel must be greater than or equal to 1/(4*sqrt(2))"
0195                 << "        newDropletsPerVoxel = " << newDropletsPerVoxel;
0196         G4Exception("FastAerosol::SetDropletsPerVoxel()", "GeomSolids0002",
0197                     FatalErrorInArgument, message);
0198   }
0199 }
0200
0201
0202 inline
0203 void FastAerosol::PrintPopulationReport() {
0204     G4cout << "Total grids: " << fNumGridCells << G4endl;
0205     G4cout << "Droplets created: " << fNumDroplets << G4endl;
0206     G4cout << "Average Number density: " << fAvgNumDens << G4endl;
0207     G4cout << "Droplets expected: " << GetNumDroplets() << G4endl;
0208 }
0209
0210
0211 // =======
0212 // Private
0213 // =======
0214 // Find the grid associated with a point.  Return true if that is a valid grid,
0215 // false if it is out of bounds.
0216 inline
0217 bool FastAerosol::GetGrid(const G4ThreeVector &p, int &xGrid, int &yGrid, int &zGrid) {
0218     xGrid = (int)floorl((p.x() + fDx) / fGridPitch);
0219     yGrid = (int)floorl((p.y() + fDy) / fGridPitch);
0220     zGrid = (int)floorl((p.z() + fDz) / fGridPitch);
0221     return(!AnyIndexOutOfBounds(xGrid, yGrid, zGrid));
0222 }
0223
0224 // Return true if any grid index given is out of bounds, false otherwise
0225 inline
0226 bool FastAerosol::AnyIndexOutOfBounds(G4int xGrid, G4int yGrid, G4int zGrid) {
0227     return ((xGrid < 0) || (xGrid >= fNx) ||
0228             (yGrid < 0) || (yGrid >= fNy) ||
0229             (zGrid < 0) || (zGrid >= fNz));
0230 }
0231
0232 // Create index for grid
0233 inline
0234 unsigned int FastAerosol::GetGridIndex(unsigned int xi, unsigned int yi, unsigned int zi) {
0235     return(zi*fNxy + yi*fNx + xi);
0236 }
0237
0238 // Create get coordinates of index
0239 inline
0240 G4ThreeVector FastAerosol::GetIndexCoord(G4int index) {
0241     G4int x = index % fNx;
0242     G4int y = ( (index -x)/fNx ) % fNy;
0243     G4int z = (index -x -fNx*y )/(fNx*fNy);
0244     return(G4ThreeVector(x,y,z));
0245 }
0246
0247 // find lower and upper grid boundary
0248 inline
0249 std::pair<G4int, G4int> FastAerosol::GetMinMaxSide(G4int i, G4int numGrids) {
0250     return(std::make_pair((i == 0) ? 0 : (i-1),
0251                            (i == (numGrids-1)) ? i : (i+1)));
0252 }