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