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