Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/fastAerosol/src/FastAerosolSolid.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/fastAerosol/src/FastAerosolSolid.cc (Version 11.3.0) and /examples/advanced/fastAerosol/src/FastAerosolSolid.cc (Version 11.2.1)


  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 FastAerosolSolid class       28 // Implementation for FastAerosolSolid 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 "FastAerosolSolid.hh"                     32 #include "FastAerosolSolid.hh"
 33                                                    33 
 34 #include "G4SystemOfUnits.hh"                      34 #include "G4SystemOfUnits.hh"
 35                                                    35 
 36 // calculate extent                                36 // calculate extent
 37 #include "G4BoundingEnvelope.hh"                   37 #include "G4BoundingEnvelope.hh"
 38 #include "G4AffineTransform.hh"                    38 #include "G4AffineTransform.hh"
 39 #include "G4VoxelLimits.hh"                        39 #include "G4VoxelLimits.hh"
 40                                                    40 
 41 // visualization                                   41 // visualization
 42 #include "G4VGraphicsScene.hh"                     42 #include "G4VGraphicsScene.hh"
 43 #include "G4VisExtent.hh"                          43 #include "G4VisExtent.hh"
 44                                                    44 
 45 // polyhedron                                      45 // polyhedron
 46 #include "G4AutoLock.hh"                           46 #include "G4AutoLock.hh"
 47 #include "G4Polyhedron.hh"                         47 #include "G4Polyhedron.hh"
 48 #include "HepPolyhedronProcessor.h"                48 #include "HepPolyhedronProcessor.h"
 49                                                    49 
 50 namespace                                          50 namespace
 51 {                                                  51 {
 52  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER     52  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
 53 }                                                  53 }
 54                                                    54 
 55                                                    55 
 56 //////////////////////////////////////////////     56 ///////////////////////////////////////////////////////////////////////////////
 57 //                                                 57 //
 58 // Constructor                                     58 // Constructor
 59 //                                                 59 //
 60 FastAerosolSolid::FastAerosolSolid(const G4Str     60 FastAerosolSolid::FastAerosolSolid(const G4String& pName,                FastAerosol* pCloud,                    G4VSolid* pDroplet,
 61                      std::function<G4RotationM     61                      std::function<G4RotationMatrix (G4ThreeVector)> pRotation)
 62   : G4VSolid(pName), fCloud(pCloud), fDroplet(     62   : G4VSolid(pName), fCloud(pCloud), fDroplet(pDroplet), fRotation(pRotation), fRebuildPolyhedron(false), fpPolyhedron(nullptr)
 63 {                                                  63 {
 64   // Get cloud size from fCloud                    64   // Get cloud size from fCloud
 65   G4ThreeVector cloudPMin, cloudPMax;              65   G4ThreeVector cloudPMin, cloudPMax;
 66   fCloud->GetBoundingLimits(cloudPMin, cloudPM     66   fCloud->GetBoundingLimits(cloudPMin, cloudPMax);
 67                                                    67 
 68   fVisDx = cloudPMax.x();                          68   fVisDx = cloudPMax.x();
 69   fVisDy = cloudPMax.y();                          69   fVisDy = cloudPMax.y();
 70   fVisDz = cloudPMax.z();                          70   fVisDz = cloudPMax.z();
 71                                                    71   
 72   // Check and set droplet radius                  72   // Check and set droplet radius
 73   G4double pR = fCloud->GetRadius();               73   G4double pR = fCloud->GetRadius();
 74   // would be nice to add a check to make sure     74   // would be nice to add a check to make sure pDroplet fits in sphere of radius pR
 75   fR = pR;                                         75   fR = pR;
 76                                                    76 
 77   fBulk = fCloud->GetBulk();                       77   fBulk = fCloud->GetBulk(); 
 78                                                    78 
 79   farFromCloudDist = fCloud->GetPreSphereR()*f     79   farFromCloudDist = fCloud->GetPreSphereR()*fR;
 80 }                                                  80 }
 81                                                    81 
 82                                                    82 
 83 //////////////////////////////////////////////     83 ///////////////////////////////////////////////////////////////////////////////
 84 //                                                 84 //
 85 // Alternative constructor (constant rotation      85 // Alternative constructor (constant rotation function)
 86 //                                                 86 //
 87 FastAerosolSolid::FastAerosolSolid(const G4Str     87 FastAerosolSolid::FastAerosolSolid(const G4String& pName,
 88                      FastAerosol* pCloud,          88                      FastAerosol* pCloud,
 89                      G4VSolid* pDroplet):          89                      G4VSolid* pDroplet):
 90   FastAerosolSolid(pName, pCloud, pDroplet,        90   FastAerosolSolid(pName, pCloud, pDroplet,
 91            [](G4ThreeVector) {return G4Rotatio     91            [](G4ThreeVector) {return G4RotationMatrix();})
 92 {}                                                 92 {}
 93                                                    93 
 94 //////////////////////////////////////////////     94 ///////////////////////////////////////////////////////////////////////////////
 95 //                                                 95 //
 96 // Fake default constructor - sets only member     96 // Fake default constructor - sets only member data and allocates memory
 97 //                            for usage restri     97 //                            for usage restricted to object persistency.
 98 //                                                 98 //
 99 FastAerosolSolid::FastAerosolSolid( __void__&      99 FastAerosolSolid::FastAerosolSolid( __void__& a )
100   : G4VSolid(a), fCloud(nullptr), fDroplet(nul    100   : G4VSolid(a), fCloud(nullptr), fDroplet(nullptr),
101   fBulk(nullptr), fR(0.),                         101   fBulk(nullptr), fR(0.),
102   fVisDx(0.), fVisDy(0.), fVisDz(0.),             102   fVisDx(0.), fVisDy(0.), fVisDz(0.),
103   fCubicVolume(0.), fSurfaceArea(0.),             103   fCubicVolume(0.), fSurfaceArea(0.),
104   farFromCloudDist(0.),                           104   farFromCloudDist(0.),
105   fRotation([](G4ThreeVector) {return G4Rotati    105   fRotation([](G4ThreeVector) {return G4RotationMatrix();}),
106   fRebuildPolyhedron(false), fpPolyhedron(null    106   fRebuildPolyhedron(false), fpPolyhedron(nullptr)
107 {                                                 107 {
108 }                                                 108 }
109                                                   109 
110 //////////////////////////////////////////////    110 ///////////////////////////////////////////////////////////////////////////////
111 //                                                111 //
112 // Copy constructor                               112 // Copy constructor
113 //                                                113 //
114 FastAerosolSolid::FastAerosolSolid(const FastA    114 FastAerosolSolid::FastAerosolSolid(const FastAerosolSolid &rhs)
115   : G4VSolid(rhs), fCloud(rhs.fCloud), fDrople    115   : G4VSolid(rhs), fCloud(rhs.fCloud), fDroplet(rhs.fDroplet),
116   fBulk(rhs.fBulk), fR(rhs.fR),                   116   fBulk(rhs.fBulk), fR(rhs.fR),
117   fVisDx(rhs.fVisDx), fVisDy(rhs.fVisDy), fVis    117   fVisDx(rhs.fVisDx), fVisDy(rhs.fVisDy), fVisDz(rhs.fVisDz),
118   fCubicVolume(rhs.fCubicVolume), fSurfaceArea    118   fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
119   farFromCloudDist(rhs.farFromCloudDist),         119   farFromCloudDist(rhs.farFromCloudDist),
120   fRotation(rhs.fRotation),                       120   fRotation(rhs.fRotation),
121   fRebuildPolyhedron(rhs.fRebuildPolyhedron),     121   fRebuildPolyhedron(rhs.fRebuildPolyhedron), fpPolyhedron(rhs.fpPolyhedron)
122 {                                                 122 {
123 }                                                 123 }
124                                                   124 
125 //////////////////////////////////////////////    125 //////////////////////////////////////////////////////////////////////////
126 //                                                126 //
127 // Assignment operator                            127 // Assignment operator
128 //                                                128 //
129 FastAerosolSolid &FastAerosolSolid::operator =    129 FastAerosolSolid &FastAerosolSolid::operator = (const FastAerosolSolid &rhs)
130 {                                                 130 {
131   // Check assignment to self                     131   // Check assignment to self
132   //                                              132   //
133   if (this == &rhs)                               133   if (this == &rhs)
134   {                                               134   {
135     return *this;                                 135     return *this;
136   }                                               136   }
137                                                   137 
138   // Copy base class data                         138   // Copy base class data
139   //                                              139   //
140   G4VSolid::operator=(rhs);                       140   G4VSolid::operator=(rhs);
141                                                   141 
142   // Copy data                                    142   // Copy data
143   //                                              143   //
144   fCloud = rhs.fCloud;                            144   fCloud = rhs.fCloud;
145   fDroplet = rhs.fDroplet;                        145   fDroplet = rhs.fDroplet;
146   fBulk = rhs.fBulk;                              146   fBulk = rhs.fBulk;
147   fR = rhs.fR;                                    147   fR = rhs.fR;
148   fVisDx = rhs.fVisDx;                            148   fVisDx = rhs.fVisDx;
149   fVisDy = rhs.fVisDy;                            149   fVisDy = rhs.fVisDy;
150   fVisDz = rhs.fVisDz;                            150   fVisDz = rhs.fVisDz;
151   fCubicVolume = rhs.fCubicVolume;                151   fCubicVolume = rhs.fCubicVolume;
152   fSurfaceArea = rhs.fSurfaceArea;                152   fSurfaceArea = rhs.fSurfaceArea;
153   farFromCloudDist = rhs.farFromCloudDist;        153   farFromCloudDist = rhs.farFromCloudDist;
154   fRotation = rhs.fRotation;                      154   fRotation = rhs.fRotation;
155   fRebuildPolyhedron = rhs.fRebuildPolyhedron;    155   fRebuildPolyhedron = rhs.fRebuildPolyhedron;
156   fpPolyhedron = rhs.fpPolyhedron;                156   fpPolyhedron = rhs.fpPolyhedron;
157                                                   157 
158                                                   158 
159   return *this;                                   159   return *this;
160 }                                                 160 }
161                                                   161 
162 //////////////////////////////////////////////    162 ///////////////////////////////////////////////////////////////////////////////
163 //                                                163 //
164 // Calculate extent under transform and specif    164 // Calculate extent under transform and specified limit
165 //                                                165 //
166 G4bool FastAerosolSolid::CalculateExtent(const    166 G4bool FastAerosolSolid::CalculateExtent(const EAxis pAxis,
167                       const G4VoxelLimits &pVo    167                       const G4VoxelLimits &pVoxelLimit,
168                       const G4AffineTransform     168                       const G4AffineTransform &pTransform,
169                       G4double &pMin, G4double    169                       G4double &pMin, G4double &pMax) const
170 {                                                 170 {
171   // Get smallest box to fully contain the clo    171   // Get smallest box to fully contain the cloud of objects, not just the centers
172   //                                              172   //
173   G4ThreeVector bmin, bmax;                       173   G4ThreeVector bmin, bmax;
174   fCloud->GetBoundingLimits(bmin, bmax);          174   fCloud->GetBoundingLimits(bmin, bmax);
175                                                   175 
176   // Find extent                                  176   // Find extent
177   //                                              177   //
178   G4BoundingEnvelope bbox(bmin, bmax);            178   G4BoundingEnvelope bbox(bmin, bmax);
179   return bbox.CalculateExtent(pAxis, pVoxelLim    179   return bbox.CalculateExtent(pAxis, pVoxelLimit, pTransform, pMin, pMax);
180 }                                                 180 }
181                                                   181 
182 //////////////////////////////////////////////    182 ///////////////////////////////////////////////////////////////////////////////
183 //                                                183 //
184 // Return whether point inside/outside/on surf    184 // Return whether point inside/outside/on surface
185 //                                                185 //
186 // This function assumes the cloud has at leas    186 // This function assumes the cloud has at least 1 droplet
187 //                                                187 //
188 EInside FastAerosolSolid::Inside(const G4Three    188 EInside FastAerosolSolid::Inside(const G4ThreeVector &p) const
189 {                                                 189 {
190   G4ThreeVector center;                           190   G4ThreeVector center;
191   G4double closestDistance;                       191   G4double closestDistance;
192                                                   192 
193   fCloud->GetNearestDroplet(p, center, closest    193   fCloud->GetNearestDroplet(p, center, closestDistance, fR, fDroplet, fRotation);
194                                                   194 
195   if (closestDistance==0.0)                       195   if (closestDistance==0.0)
196   {                                               196   {
197     G4RotationMatrix irotm = fRotation(center)    197     G4RotationMatrix irotm = fRotation(center).inverse();
198                                                   198 
199     return fDroplet->Inside( irotm*(p - center    199     return fDroplet->Inside( irotm*(p - center) );
200   }                                               200   }
201   else                                            201   else
202   {                                               202   {
203     return kOutside;                              203     return kOutside;
204   }                                               204   }
205 }                                                 205 }
206                                                   206 
207 //////////////////////////////////////////////    207 /////////////////////////////////////////////////////////////////////
208 //                                                208 //
209 // Return unit normal of surface closest to p     209 // Return unit normal of surface closest to p
210 //                                                210 //
211 // This function assumes the cloud has at leas    211 // This function assumes the cloud has at least 1 droplet
212 //                                                212 //
213 G4ThreeVector FastAerosolSolid::SurfaceNormal(    213 G4ThreeVector FastAerosolSolid::SurfaceNormal(const G4ThreeVector &p) const
214 {                                                 214 {
215   G4ThreeVector center;                           215   G4ThreeVector center;
216   G4double closestDistance;                       216   G4double closestDistance;
217                                                   217 
218   fCloud->GetNearestDroplet(p, center, closest    218   fCloud->GetNearestDroplet(p, center, closestDistance, DBL_MAX, fDroplet, fRotation);
219                                                   219 
220   G4RotationMatrix rotm = fRotation(center);      220   G4RotationMatrix rotm = fRotation(center);
221                                                   221 
222   return rotm*( fDroplet->SurfaceNormal( rotm.    222   return rotm*( fDroplet->SurfaceNormal( rotm.inverse()*(p - center) ) );
223 }                                                 223 }
224                                                   224 
225 //////////////////////////////////////////////    225 ///////////////////////////////////////////////////////////////////////////////
226 //                                                226 //
227 // Calculate distance to shape from outside, a    227 // Calculate distance to shape from outside, along normalised vector
228 //                                                228 //
229 // This CANNOT be an underestimate                229 // This CANNOT be an underestimate
230 //                                                230 //
231 G4double FastAerosolSolid::DistanceToIn(const     231 G4double FastAerosolSolid::DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
232 {                                                 232 {
233   G4ThreeVector center;                           233   G4ThreeVector center;
234   G4double closestDistance;                       234   G4double closestDistance;
235                                                   235 
236   if (fCloud->GetNearestDroplet(p, v, center,     236   if (fCloud->GetNearestDroplet(p, v, center, closestDistance, fStepLim, fDroplet, fRotation))  // if we found a droplet within fStepLim of query
237   {                                               237   {
238     return closestDistance;                       238     return closestDistance;
239   }                                               239   }
240   else if (fCloud->DistanceToCloud(p,v)<DBL_MA    240   else if (fCloud->DistanceToCloud(p,v)<DBL_MAX)      // if there is cloud in front of us
241   {                                               241   {
242     return 1.1*fStepLim;                          242     return 1.1*fStepLim;
243   }                                               243   }
244   else                          // flying away    244   else                          // flying away from cloud
245   {                                               245   {
246     return kInfinity;                             246     return kInfinity;
247   }                                               247   }
248 }                                                 248 }
249                                                   249 
250 //////////////////////////////////////////////    250 //////////////////////////////////////////////////////////////////////
251 //                                                251 //
252 // Calculate distance (<= actual) to closest s    252 // Calculate distance (<= actual) to closest surface of shape from outside
253 //                                                253 //
254 // This function assumes the cloud has at leas    254 // This function assumes the cloud has at least 1 droplet
255 //                                                255 //
256 // This can be an underestimate                   256 // This can be an underestimate
257 //                                                257 //
258 G4double FastAerosolSolid::DistanceToIn(const     258 G4double FastAerosolSolid::DistanceToIn(const G4ThreeVector &p) const
259 {                                                 259 {
260   G4ThreeVector center;                           260   G4ThreeVector center;
261   G4double closestDistance;                       261   G4double closestDistance;
262                                                   262 
263   G4double distanceToCloud = fBulk->DistanceTo    263   G4double distanceToCloud = fBulk->DistanceToIn(p);
264                                                   264 
265   if (fBulk->Inside(p)==kOutside && distanceTo    265   if (fBulk->Inside(p)==kOutside && distanceToCloud>=farFromCloudDist)
266   {                                               266   {
267     return distanceToCloud;                       267     return distanceToCloud;
268   }                                               268   }
269   else if (fCloud->GetNearestDroplet(p, center    269   else if (fCloud->GetNearestDroplet(p, center, closestDistance, fStepLim, fDroplet, fRotation))    // if we found a droplet within fStepLim of query
270   {                                               270   {
271     return closestDistance;                       271     return closestDistance;
272   }                                               272   }
273   else                                            273   else
274   {                                               274   {
275     return 1.1*fStepLim;                          275     return 1.1*fStepLim;
276   }                                               276   }
277 }                                                 277 }
278                                                   278 
279 //////////////////////////////////////////////    279 //////////////////////////////////////////////////////////////////////
280 //                                                280 //
281 // Calculate distance (<= actual) to closest s    281 // Calculate distance (<= actual) to closest surface of shape from inside
282 //                                                282 //
283 // Despite being a vector distance, we find th    283 // Despite being a vector distance, we find the absolutely closest
284 // droplet to our point since we assume that p    284 // droplet to our point since we assume that p is in a droplet and p
285 // could be past the center                       285 // could be past the center
286 //                                                286 //
287 // This CANNOT be an underestimate                287 // This CANNOT be an underestimate
288 //                                                288 //
289 G4double FastAerosolSolid::DistanceToOut(const    289 G4double FastAerosolSolid::DistanceToOut(const G4ThreeVector &p,
290                      const G4ThreeVector &v,      290                      const G4ThreeVector &v,
291                      const G4bool calcNorm,       291                      const G4bool calcNorm,
292                      G4bool *validNorm,           292                      G4bool *validNorm,
293                      G4ThreeVector *n) const      293                      G4ThreeVector *n) const
294 {                                                 294 {
295   G4ThreeVector center;                           295   G4ThreeVector center;
296   G4double distanceToIn; // should be 0           296   G4double distanceToIn; // should be 0
297                                                   297 
298   fCloud->GetNearestDroplet(p, center, distanc    298   fCloud->GetNearestDroplet(p, center, distanceToIn, fR, fDroplet, fRotation);  // if we call this function, must be inside and thus must have a droplet within fR
299                                                   299 
300   G4RotationMatrix rotm = fRotation(center);      300   G4RotationMatrix rotm = fRotation(center);
301   G4RotationMatrix irotm = rotm.inverse();        301   G4RotationMatrix irotm = rotm.inverse();
302                                                   302 
303   G4ThreeVector relPos = irotm*(p-center);        303   G4ThreeVector relPos = irotm*(p-center);
304                                                   304 
305   if (fDroplet->Inside(relPos) == kOutside) //    305   if (fDroplet->Inside(relPos) == kOutside) // something went wrong... we should be inside
306   {                                               306   {
307     std::ostringstream message;                   307     std::ostringstream message;
308     message << std::setprecision(15) << "The p    308     message << std::setprecision(15) << "The particle at point p = " << p/mm << "mm"
309         << std::setprecision(15) << " called D    309         << std::setprecision(15) << " called DistanceToOut(p,v) and found the closest droplet to be at center = " << center/mm << "mm"
310         << " but p is outside the droplet!";      310         << " but p is outside the droplet!";
311     G4Exception("FastAerosolSolid::DistanceToO    311     G4Exception("FastAerosolSolid::DistanceToOut()", "GeomSolids0002",
312           FatalErrorInArgument, message);         312           FatalErrorInArgument, message);
313   }                                               313   }
314                                                   314 
315   G4double dist = fDroplet->DistanceToOut(relP    315   G4double dist = fDroplet->DistanceToOut(relPos, irotm*v, calcNorm, validNorm, n);
316   *n = rotm*(*n);                                 316   *n = rotm*(*n);
317   *validNorm = false; // even if droplet is co    317   *validNorm = false; // even if droplet is convex, the aerosol isn't
318                                                   318 
319   return dist;                                    319   return dist;
320 }                                                 320 }
321                                                   321 
322 //////////////////////////////////////////////    322 /////////////////////////////////////////////////////////////////////////
323 //                                                323 //
324 // Calculate distance (<=actual) to closest su    324 // Calculate distance (<=actual) to closest surface of shape from inside
325 //                                                325 //
326 // This can be an underestimate                   326 // This can be an underestimate
327 //                                                327 //
328 G4double FastAerosolSolid::DistanceToOut(const    328 G4double FastAerosolSolid::DistanceToOut(const G4ThreeVector &p) const
329 {                                                 329 {
330   G4ThreeVector center;                           330   G4ThreeVector center;
331   G4double distanceToIn; // should be 0           331   G4double distanceToIn; // should be 0
332                                                   332 
333   fCloud->GetNearestDroplet(p, center, distanc    333   fCloud->GetNearestDroplet(p, center, distanceToIn, fR, fDroplet, fRotation);  // if we call this function, must be inside and thus must have a droplet within fR
334                                                   334 
335   G4RotationMatrix irotm = fRotation(center).i    335   G4RotationMatrix irotm = fRotation(center).inverse();
336   G4ThreeVector relPos = irotm*(p-center);        336   G4ThreeVector relPos = irotm*(p-center);
337                                                   337 
338   if (fDroplet->Inside(relPos) == kOutside) //    338   if (fDroplet->Inside(relPos) == kOutside) // something went wrong... we should be inside
339   {                                               339   {
340     std::ostringstream message;                   340     std::ostringstream message;
341     message << "The particle at point p = " <<    341     message << "The particle at point p = " << p/mm << "mm"
342         << " called DistanceToOut(p) and found    342         << " called DistanceToOut(p) and found the closest droplet to be at center = " << center/mm << "mm"
343         << " but p is outside the droplet!";      343         << " but p is outside the droplet!";
344     G4Exception("FastAerosolSolid::DistanceToO    344     G4Exception("FastAerosolSolid::DistanceToOut()", "GeomSolids0002", 
345           FatalErrorInArgument, message);         345           FatalErrorInArgument, message);
346   }                                               346   }
347                                                   347 
348   return fDroplet->DistanceToOut(relPos);         348   return fDroplet->DistanceToOut(relPos);
349 }                                                 349 }
350                                                   350 
351 //////////////////////////////////////////////    351 //////////////////////////////////////////////////////////////////////////
352 //                                                352 //
353 // G4EntityType                                   353 // G4EntityType
354 //                                                354 //
355 G4GeometryType FastAerosolSolid::GetEntityType    355 G4GeometryType FastAerosolSolid::GetEntityType() const
356 {                                                 356 {
357   return G4String("FastAerosolSolid");            357   return G4String("FastAerosolSolid");
358 }                                                 358 }
359                                                   359 
360 //////////////////////////////////////////////    360 //////////////////////////////////////////////////////////////////////////
361 //                                                361 //
362 // G4EntityType                                   362 // G4EntityType
363 //                                                363 //
364 G4VSolid* FastAerosolSolid::Clone() const         364 G4VSolid* FastAerosolSolid::Clone() const
365 {                                                 365 {
366   return new FastAerosolSolid(*this);             366   return new FastAerosolSolid(*this);
367 }                                                 367 }
368                                                   368 
369 //////////////////////////////////////////////    369 //////////////////////////////////////////////////////////////////////////
370 //                                                370 //
371 // Stream object contents to an output stream     371 // Stream object contents to an output stream
372 //                                                372 //
373 std::ostream &FastAerosolSolid::StreamInfo(std    373 std::ostream &FastAerosolSolid::StreamInfo(std::ostream &os) const
374 {                                                 374 {
375   os << "-------------------------------------    375   os << "-----------------------------------------------------------\n"
376      << "    *** Dump for solid - " << GetName    376      << "    *** Dump for solid - " << GetName() << " ***\n"
377      << "    =================================    377      << "    ===================================================\n"
378      << " Solid type: FastAerosolSolid\n"         378      << " Solid type: FastAerosolSolid\n"
379      << " Parameters: \n"                         379      << " Parameters: \n"
380      << "    numDroplets: " << fCloud->GetNumD    380      << "    numDroplets: " << fCloud->GetNumDroplets() << "\n"
381      << "    fDroplet type: " << fDroplet->Get    381      << "    fDroplet type: " << fDroplet->GetName() << "\n"
382      << "    fDroplet parameters: \n";            382      << "    fDroplet parameters: \n";
383      fDroplet->StreamInfo(os);                    383      fDroplet->StreamInfo(os);
384     os  << "----------------------------------    384     os  << "-----------------------------------------------------------\n";
385   return os;                                      385   return os;
386 }                                                 386 }
387                                                   387 
388 //////////////////////////////////////////////    388 ////////////////////////////////////////////////////////////////////////////////
389 //                                                389 //
390 // GetPointOnSurface                              390 // GetPointOnSurface
391 //                                                391 //
392 // Currently hardcoded to look at all droplets    392 // Currently hardcoded to look at all droplets, not just the populated ones
393 //                                                393 //
394 G4ThreeVector FastAerosolSolid::GetPointOnSurf    394 G4ThreeVector FastAerosolSolid::GetPointOnSurface() const
395 {                                                 395 {
396   G4ThreeVector center;                           396   G4ThreeVector center;
397   G4double closestDistance;                       397   G4double closestDistance;
398                                                   398 
399   G4double fDx = fCloud->GetXHalfLength();        399   G4double fDx = fCloud->GetXHalfLength();
400   G4double fDy = fCloud->GetYHalfLength();        400   G4double fDy = fCloud->GetYHalfLength();
401   G4double fDz = fCloud->GetZHalfLength();        401   G4double fDz = fCloud->GetZHalfLength();
402                                                   402 
403   G4ThreeVector p(2.0*fDx*G4UniformRand(),2.0*    403   G4ThreeVector p(2.0*fDx*G4UniformRand(),2.0*fDy*G4UniformRand(),2.0*fDz*G4UniformRand());
404   p -= G4ThreeVector(fDx, fDy, fDz);              404   p -= G4ThreeVector(fDx, fDy, fDz);
405                                                   405 
406   fCloud->GetNearestDroplet(p, center, closest    406   fCloud->GetNearestDroplet(p, center, closestDistance, DBL_MAX, fDroplet, fRotation);
407                                                   407 
408   return(center + fRotation(center)*fDroplet->    408   return(center + fRotation(center)*fDroplet->GetPointOnSurface());
409 }                                                 409 }
410                                                   410 
411                                                   411 
412 //////////////////////////////////////////////    412 /////////////////////////////////////////////////////////////////////////////
413 //                                                413 //
414 // Methods for visualisation                      414 // Methods for visualisation
415 //                                                415 //
416 void FastAerosolSolid::DescribeYourselfTo (G4V    416 void FastAerosolSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const
417 {                                                 417 {
418   scene.AddSolid(*this);                          418   scene.AddSolid(*this);
419 }                                                 419 }
420                                                   420 
421 G4VisExtent FastAerosolSolid::GetExtent() cons    421 G4VisExtent FastAerosolSolid::GetExtent() const
422 {                                                 422 {
423   return G4VisExtent (-fVisDx, fVisDx, -fVisDy    423   return G4VisExtent (-fVisDx, fVisDx, -fVisDy, fVisDy, -fVisDz, fVisDz);
424 }                                                 424 }
425                                                   425 
426 G4Polyhedron* FastAerosolSolid::CreatePolyhedr    426 G4Polyhedron* FastAerosolSolid::CreatePolyhedron () const
427 {                                                 427 {
428   return fBulk->CreatePolyhedron();               428   return fBulk->CreatePolyhedron();
429 }                                                 429 }
430                                                   430 
431                                                   431 
432 // copied from G4Ellipsoid                        432 // copied from G4Ellipsoid
433 G4Polyhedron* FastAerosolSolid::GetPolyhedron     433 G4Polyhedron* FastAerosolSolid::GetPolyhedron () const
434 {                                                 434 {
435   if (!fpPolyhedron ||                            435   if (!fpPolyhedron ||
436     fRebuildPolyhedron ||                         436     fRebuildPolyhedron ||
437     fpPolyhedron->GetNumberOfRotationStepsAtTi    437     fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
438     fpPolyhedron->GetNumberOfRotationSteps())     438     fpPolyhedron->GetNumberOfRotationSteps())
439   {                                               439   {
440     G4AutoLock l(&polyhedronMutex);               440     G4AutoLock l(&polyhedronMutex);
441     delete fpPolyhedron;                          441     delete fpPolyhedron;
442     fpPolyhedron = CreatePolyhedron();            442     fpPolyhedron = CreatePolyhedron();
443     fRebuildPolyhedron = false;                   443     fRebuildPolyhedron = false;
444     l.unlock();                                   444     l.unlock();
445   }                                               445   }
446   return fpPolyhedron;                            446   return fpPolyhedron;
447 }                                                 447 }
448                                                   448