| 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 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