| 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 // G4VSolid implementation for solid base clas <<
27 // 26 //
28 // 10.10.18 E.Tcherniaev, more robust Estimate << 27 // $Id: G4VSolid.cc,v 1.39 2008/09/23 13:07:41 gcosmo Exp $
29 // 30.06.95 P.Kent, Created. << 28 // GEANT4 tag $Name: geant4-09-02-patch-04 $
>> 29 //
>> 30 // class G4VSolid
>> 31 //
>> 32 // Implementation for solid base class
>> 33 //
>> 34 // History:
>> 35 //
>> 36 // 06.12.02 V.Grichine, restored original conditions in ClipPolygon()
>> 37 // 10.05.02 V.Grichine, ClipPolygon(): clip only other axis and limited voxels
>> 38 // 15.04.02 V.Grichine, bug fixed in ClipPolygon(): clip only one axis
>> 39 // 13.03.02 V.Grichine, cosmetics of voxel limit functions
>> 40 // 15.11.00 D.Williams, V.Grichine, fix in CalculateClippedPolygonExtent()
>> 41 // 10.07.95 P.Kent, Added == operator, solid Store entry
>> 42 // 30.06.95 P.Kent, Created.
30 // ------------------------------------------- 43 // --------------------------------------------------------------------
31 44
32 #include "G4VSolid.hh" 45 #include "G4VSolid.hh"
33 #include "G4SolidStore.hh" 46 #include "G4SolidStore.hh"
34 #include "globals.hh" 47 #include "globals.hh"
35 #include "G4QuickRand.hh" << 48 #include "Randomize.hh"
36 #include "G4GeometryTolerance.hh" 49 #include "G4GeometryTolerance.hh"
37 50
38 #include "G4VoxelLimits.hh" 51 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh" 52 #include "G4AffineTransform.hh"
40 #include "G4VisExtent.hh" 53 #include "G4VisExtent.hh"
41 54
42 ////////////////////////////////////////////// 55 //////////////////////////////////////////////////////////////////////////
43 // 56 //
44 // Streaming operator dumping solid contents <<
45 <<
46 std::ostream& operator<< ( std::ostream& os, c <<
47 { <<
48 return e.StreamInfo(os); <<
49 } <<
50 <<
51 ////////////////////////////////////////////// <<
52 // <<
53 // Constructor 57 // Constructor
54 // - Copies name 58 // - Copies name
55 // - Add ourselves to solid Store 59 // - Add ourselves to solid Store
56 60
57 G4VSolid::G4VSolid(const G4String& name) 61 G4VSolid::G4VSolid(const G4String& name)
58 : fshapeName(name) 62 : fshapeName(name)
59 { 63 {
60 kCarTolerance = G4GeometryTolerance::GetIn 64 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
61 65
62 // Register to store 66 // Register to store
63 // 67 //
64 G4SolidStore::GetInstance()->Register(this 68 G4SolidStore::GetInstance()->Register(this);
65 } 69 }
66 70
67 ////////////////////////////////////////////// 71 //////////////////////////////////////////////////////////////////////////
68 // 72 //
69 // Copy constructor 73 // Copy constructor
70 // 74 //
71 75
72 G4VSolid::G4VSolid(const G4VSolid& rhs) 76 G4VSolid::G4VSolid(const G4VSolid& rhs)
73 : kCarTolerance(rhs.kCarTolerance), fshapeNa 77 : kCarTolerance(rhs.kCarTolerance), fshapeName(rhs.fshapeName)
74 { 78 {
75 // Register to store 79 // Register to store
76 // 80 //
77 G4SolidStore::GetInstance()->Register(this 81 G4SolidStore::GetInstance()->Register(this);
78 } 82 }
79 83
80 ////////////////////////////////////////////// 84 //////////////////////////////////////////////////////////////////////////
81 // 85 //
82 // Fake default constructor - sets only member 86 // Fake default constructor - sets only member data and allocates memory
83 // for usage restri 87 // for usage restricted to object persistency.
84 // 88 //
85 G4VSolid::G4VSolid( __void__& ) 89 G4VSolid::G4VSolid( __void__& )
86 : fshapeName("") 90 : fshapeName("")
87 { 91 {
88 // Register to store 92 // Register to store
89 // 93 //
90 G4SolidStore::GetInstance()->Register(this 94 G4SolidStore::GetInstance()->Register(this);
91 } 95 }
92 96
93 ////////////////////////////////////////////// 97 //////////////////////////////////////////////////////////////////////////
94 // 98 //
95 // Destructor (virtual) 99 // Destructor (virtual)
96 // - Remove ourselves from solid Store 100 // - Remove ourselves from solid Store
97 101
98 G4VSolid::~G4VSolid() 102 G4VSolid::~G4VSolid()
99 { 103 {
100 G4SolidStore::GetInstance()->DeRegister(th 104 G4SolidStore::GetInstance()->DeRegister(this);
101 } 105 }
102 106
103 ////////////////////////////////////////////// 107 //////////////////////////////////////////////////////////////////////////
104 // 108 //
105 // Assignment operator 109 // Assignment operator
106 110
107 G4VSolid& G4VSolid::operator = (const G4VSolid << 111 G4VSolid& G4VSolid::operator = (const G4VSolid& rhs)
108 { 112 {
109 // Check assignment to self 113 // Check assignment to self
110 // 114 //
111 if (this == &rhs) { return *this; } 115 if (this == &rhs) { return *this; }
112 116
113 // Copy data 117 // Copy data
114 // 118 //
115 kCarTolerance = rhs.kCarTolerance; 119 kCarTolerance = rhs.kCarTolerance;
116 fshapeName = rhs.fshapeName; 120 fshapeName = rhs.fshapeName;
117 121
118 return *this; 122 return *this;
119 } << 123 }
120 <<
121 <<
122 124
123 ////////////////////////////////////////////// 125 //////////////////////////////////////////////////////////////////////////
124 // 126 //
125 // Set solid name and notify store of the chan << 127 // Streaming operator dumping solid contents
126 128
127 void G4VSolid::SetName(const G4String& name) << 129 std::ostream& operator<< ( std::ostream& os, const G4VSolid& e )
128 { 130 {
129 fshapeName = name; << 131 return e.StreamInfo(os);
130 G4SolidStore::GetInstance()->SetMapValid(fal <<
131 } 132 }
132 133
133 ////////////////////////////////////////////// 134 //////////////////////////////////////////////////////////////////////////
134 // 135 //
135 // Throw exception if ComputeDimensions called 136 // Throw exception if ComputeDimensions called for illegal derived class
136 137
137 void G4VSolid::ComputeDimensions(G4VPVParamete 138 void G4VSolid::ComputeDimensions(G4VPVParameterisation*,
138 const G4int, 139 const G4int,
139 const G4VPhys 140 const G4VPhysicalVolume*)
140 { 141 {
141 std::ostringstream message; << 142 G4cerr << "ERROR - Illegal call to G4VSolid::ComputeDimensions()" << G4endl
142 message << "Illegal call to G4VSolid::Comp << 143 << " Method not overloaded by derived class !" << G4endl;
143 << "Method not overloaded by deriv << 144 G4Exception("G4VSolid::ComputeDimensions()", "NotApplicable",
144 G4Exception("G4VSolid::ComputeDimensions() << 145 FatalException, "Illegal call to case class.");
145 FatalException, message); <<
146 } 146 }
147 147
148 ////////////////////////////////////////////// 148 //////////////////////////////////////////////////////////////////////////
149 // 149 //
150 // Throw exception (warning) for solids not im 150 // Throw exception (warning) for solids not implementing the method
151 151
152 G4ThreeVector G4VSolid::GetPointOnSurface() co 152 G4ThreeVector G4VSolid::GetPointOnSurface() const
153 { 153 {
154 std::ostringstream message; << 154 G4cerr << "WARNING - G4VSolid::GetPointOnSurface()" << G4endl
155 message << "Not implemented for solid: " << 155 << " Not implemented for solid: "
156 << GetEntityType() << " !" << G4en << 156 << this->GetEntityType() << " !" << G4endl;
157 << "Returning origin."; << 157 G4Exception("G4VSolid::GetPointOnSurface()", "NotImplemented",
158 G4Exception("G4VSolid::GetPointOnSurface() << 158 JustWarning, "Not implemented for this solid ! Returning origin.");
159 JustWarning, message); << 159 return G4ThreeVector(0,0,0);
160 return {0,0,0}; <<
161 } <<
162 <<
163 ////////////////////////////////////////////// <<
164 // <<
165 // Returns total number of constituents that w <<
166 // of the solid. For non-Boolean solids the re <<
167 <<
168 G4int G4VSolid::GetNumOfConstituents() const <<
169 { return 1; } <<
170 <<
171 ////////////////////////////////////////////// <<
172 // <<
173 // Returns true if the solid has only planar f <<
174 <<
175 G4bool G4VSolid::IsFaceted() const <<
176 { return false; } <<
177 <<
178 ////////////////////////////////////////////// <<
179 // <<
180 // Dummy implementations ... <<
181 <<
182 const G4VSolid* G4VSolid::GetConstituentSolid( <<
183 { return nullptr; } <<
184 <<
185 G4VSolid* G4VSolid::GetConstituentSolid(G4int) <<
186 { return nullptr; } <<
187 <<
188 const G4DisplacedSolid* G4VSolid::GetDisplaced <<
189 { return nullptr; } <<
190 <<
191 G4DisplacedSolid* G4VSolid::GetDisplacedSolidP <<
192 { return nullptr; } <<
193 <<
194 ////////////////////////////////////////////// <<
195 // <<
196 // Returns an estimation of the solid volume i <<
197 // The number of statistics and error accuracy <<
198 // This method may be overloaded by derived cl <<
199 // exact geometrical quantity for solids where <<
200 // or anyway to cache the computed value. <<
201 // This implementation does NOT cache the comp <<
202 <<
203 G4double G4VSolid::GetCubicVolume() <<
204 { <<
205 G4int cubVolStatistics = 1000000; <<
206 G4double cubVolEpsilon = 0.001; <<
207 return EstimateCubicVolume(cubVolStatistics, <<
208 } <<
209 <<
210 ////////////////////////////////////////////// <<
211 // <<
212 // Calculate cubic volume based on Inside() me <<
213 // Accuracy is limited by the second argument <<
214 // expressed by the first argument. <<
215 // Implementation is courtesy of Vasiliki Desp <<
216 // University of Athens. <<
217 <<
218 G4double G4VSolid::EstimateCubicVolume(G4int n <<
219 { <<
220 G4int iInside=0; <<
221 G4double px,py,pz,minX,maxX,minY,maxY,minZ,m <<
222 G4ThreeVector p; <<
223 EInside in; <<
224 <<
225 // values needed for CalculateExtent signatu <<
226 <<
227 G4VoxelLimits limit; // Unlim <<
228 G4AffineTransform origin; <<
229 <<
230 // min max extents of pSolid along X,Y,Z <<
231 <<
232 CalculateExtent(kXAxis,limit,origin,minX,max <<
233 CalculateExtent(kYAxis,limit,origin,minY,max <<
234 CalculateExtent(kZAxis,limit,origin,minZ,max <<
235 <<
236 // limits <<
237 <<
238 if(nStat < 100) nStat = 100; <<
239 if(epsilon > 0.01) epsilon = 0.01; <<
240 halfepsilon = 0.5*epsilon; <<
241 <<
242 for(auto i = 0; i < nStat; ++i ) <<
243 { <<
244 px = minX-halfepsilon+(maxX-minX+epsilon)* <<
245 py = minY-halfepsilon+(maxY-minY+epsilon)* <<
246 pz = minZ-halfepsilon+(maxZ-minZ+epsilon)* <<
247 p = G4ThreeVector(px,py,pz); <<
248 in = Inside(p); <<
249 if(in != kOutside) ++iInside; <<
250 } <<
251 volume = (maxX-minX+epsilon)*(maxY-minY+epsi <<
252 * (maxZ-minZ+epsilon)*iInside/nStat; <<
253 return volume; <<
254 } <<
255 <<
256 ////////////////////////////////////////////// <<
257 // <<
258 // Returns an estimation of the solid surface <<
259 // The number of statistics and error accuracy <<
260 // This method may be overloaded by derived cl <<
261 // exact geometrical quantity for solids where <<
262 // or anyway to cache the computed value. <<
263 // This implementation does NOT cache the comp <<
264 <<
265 G4double G4VSolid::GetSurfaceArea() <<
266 { <<
267 G4int stat = 1000000; <<
268 G4double ell = -1.; <<
269 return EstimateSurfaceArea(stat,ell); <<
270 } <<
271 <<
272 ////////////////////////////////////////////// <<
273 // <<
274 // Calculate surface area by estimating volume <<
275 // surrounding the surface using Monte-Carlo m <<
276 // Input parameters: <<
277 // nstat - statistics (number of random poi <<
278 // eps - shell thinkness <<
279 <<
280 G4double G4VSolid::EstimateSurfaceArea(G4int n <<
281 { <<
282 static const G4double s2 = 1./std::sqrt(2.); <<
283 static const G4double s3 = 1./std::sqrt(3.); <<
284 static const G4ThreeVector directions[64] = <<
285 { <<
286 G4ThreeVector( 0, 0, 0), G4ThreeVector( <<
287 G4ThreeVector( 1, 0, 0), G4ThreeVector( <<
288 G4ThreeVector( 0, -1, 0), G4ThreeVector( <<
289 G4ThreeVector( s2, -s2, 0), G4ThreeVector( <<
290 <<
291 G4ThreeVector( 0, 1, 0), G4ThreeVector( <<
292 G4ThreeVector( s2, s2, 0), G4ThreeVector( <<
293 G4ThreeVector( 0, -1, 0), G4ThreeVector( <<
294 G4ThreeVector( 1, 0, 0), G4ThreeVector( <<
295 <<
296 G4ThreeVector( 0, 0, -1), G4ThreeVector( <<
297 G4ThreeVector( s2, 0,-s2), G4ThreeVector( <<
298 G4ThreeVector( 0,-s2,-s2), G4ThreeVector( <<
299 G4ThreeVector( s3,-s3,-s3), G4ThreeVector( <<
300 <<
301 G4ThreeVector( 0, s2,-s2), G4ThreeVector( <<
302 G4ThreeVector( s3, s3,-s3), G4ThreeVector( <<
303 G4ThreeVector( 0, 0, -1), G4ThreeVector( <<
304 G4ThreeVector( s2, 0,-s2), G4ThreeVector( <<
305 <<
306 G4ThreeVector( 0, 0, 1), G4ThreeVector( <<
307 G4ThreeVector( s2, 0, s2), G4ThreeVector( <<
308 G4ThreeVector( 0,-s2, s2), G4ThreeVector( <<
309 G4ThreeVector( s3,-s3, s3), G4ThreeVector( <<
310 <<
311 G4ThreeVector( 0, s2, s2), G4ThreeVector( <<
312 G4ThreeVector( s3, s3, s3), G4ThreeVector( <<
313 G4ThreeVector( 0, 0, 1), G4ThreeVector( <<
314 G4ThreeVector( s2, 0, s2), G4ThreeVector( <<
315 <<
316 G4ThreeVector( 0, 0, -1), G4ThreeVector( <<
317 G4ThreeVector( 1, 0, 0), G4ThreeVector( <<
318 G4ThreeVector( 0, -1, 0), G4ThreeVector( <<
319 G4ThreeVector( s2, -s2, 0), G4ThreeVector( <<
320 <<
321 G4ThreeVector( 0, 1, 0), G4ThreeVector( <<
322 G4ThreeVector( s2, s2, 0), G4ThreeVector( <<
323 G4ThreeVector( 0, -1, 0), G4ThreeVector( <<
324 G4ThreeVector( 1, 0, 0), G4ThreeVector( <<
325 }; <<
326 <<
327 G4ThreeVector bmin, bmax; <<
328 BoundingLimits(bmin, bmax); <<
329 <<
330 G4double dX = bmax.x() - bmin.x(); <<
331 G4double dY = bmax.y() - bmin.y(); <<
332 G4double dZ = bmax.z() - bmin.z(); <<
333 <<
334 // Define statistics and shell thickness <<
335 // <<
336 G4int npoints = (nstat < 1000) ? 1000 : nsta <<
337 G4double coeff = 0.5 / std::cbrt(G4double(np <<
338 G4double eps = (ell > 0) ? ell : coeff * std <<
339 G4double del = 1.8 * eps; // shold be more t <<
340 <<
341 G4double minX = bmin.x() - eps; <<
342 G4double minY = bmin.y() - eps; <<
343 G4double minZ = bmin.z() - eps; <<
344 <<
345 G4double dd = 2. * eps; <<
346 dX += dd; <<
347 dY += dd; <<
348 dZ += dd; <<
349 <<
350 // Calculate surface area <<
351 // <<
352 G4int icount = 0; <<
353 for(auto i = 0; i < npoints; ++i) <<
354 { <<
355 G4double px = minX + dX*G4QuickRand(); <<
356 G4double py = minY + dY*G4QuickRand(); <<
357 G4double pz = minZ + dZ*G4QuickRand(); <<
358 G4ThreeVector p = G4ThreeVector(px, py, p <<
359 EInside in = Inside(p); <<
360 G4double dist = 0; <<
361 if (in == kInside) <<
362 { <<
363 if (DistanceToOut(p) >= eps) continue; <<
364 G4int icase = 0; <<
365 if (Inside(G4ThreeVector(px-del, py, pz) <<
366 if (Inside(G4ThreeVector(px+del, py, pz) <<
367 if (Inside(G4ThreeVector(px, py-del, pz) <<
368 if (Inside(G4ThreeVector(px, py+del, pz) <<
369 if (Inside(G4ThreeVector(px, py, pz-del) <<
370 if (Inside(G4ThreeVector(px, py, pz+del) <<
371 if (icase == 0) continue; <<
372 G4ThreeVector v = directions[icase]; <<
373 dist = DistanceToOut(p, v); <<
374 G4ThreeVector n = SurfaceNormal(p + v*di <<
375 dist *= v.dot(n); <<
376 } <<
377 else if (in == kOutside) <<
378 { <<
379 if (DistanceToIn(p) >= eps) continue; <<
380 G4int icase = 0; <<
381 if (Inside(G4ThreeVector(px-del, py, pz) <<
382 if (Inside(G4ThreeVector(px+del, py, pz) <<
383 if (Inside(G4ThreeVector(px, py-del, pz) <<
384 if (Inside(G4ThreeVector(px, py+del, pz) <<
385 if (Inside(G4ThreeVector(px, py, pz-del) <<
386 if (Inside(G4ThreeVector(px, py, pz+del) <<
387 if (icase == 0) continue; <<
388 G4ThreeVector v = directions[icase]; <<
389 dist = DistanceToIn(p, v); <<
390 if (dist == kInfinity) continue; <<
391 G4ThreeVector n = SurfaceNormal(p + v*di <<
392 dist *= -(v.dot(n)); <<
393 } <<
394 if (dist < eps) ++icount; <<
395 } <<
396 return dX*dY*dZ*icount/npoints/dd; <<
397 } <<
398 <<
399 ////////////////////////////////////////////// <<
400 // <<
401 // Returns a pointer of a dynamically allocate <<
402 // Returns NULL pointer with warning in case t <<
403 // implement this method. The caller has respo <<
404 // <<
405 <<
406 G4VSolid* G4VSolid::Clone() const <<
407 { <<
408 std::ostringstream message; <<
409 message << "Clone() method not implemented f <<
410 << GetEntityType() << "!" << G4endl <<
411 << "Returning NULL pointer!"; <<
412 G4Exception("G4VSolid::Clone()", "GeomMgt100 <<
413 return nullptr; <<
414 } 160 }
415 161
416 ////////////////////////////////////////////// 162 ///////////////////////////////////////////////////////////////////////////
417 // << 163 //
418 // Calculate the maximum and minimum extents o 164 // Calculate the maximum and minimum extents of the polygon described
419 // by the vertices: pSectionIndex->pSectionInd 165 // by the vertices: pSectionIndex->pSectionIndex+1->
420 // pSectionIndex+2->pSection 166 // pSectionIndex+2->pSectionIndex+3->pSectionIndex
421 // in the List pVertices 167 // in the List pVertices
422 // 168 //
423 // If the minimum is <pMin pMin is set to the 169 // If the minimum is <pMin pMin is set to the new minimum
424 // If the maximum is >pMax pMax is set to the 170 // If the maximum is >pMax pMax is set to the new maximum
425 // 171 //
426 // No modifications are made to pVertices 172 // No modifications are made to pVertices
427 // 173 //
428 174
429 void G4VSolid::ClipCrossSection( G4Three 175 void G4VSolid::ClipCrossSection( G4ThreeVectorList* pVertices,
430 const G4int p 176 const G4int pSectionIndex,
431 const G4Voxel 177 const G4VoxelLimits& pVoxelLimit,
432 const EAxis p << 178 const EAxis pAxis,
433 G4doubl 179 G4double& pMin, G4double& pMax) const
434 { 180 {
435 181
436 G4ThreeVectorList polygon; 182 G4ThreeVectorList polygon;
437 polygon.reserve(4); 183 polygon.reserve(4);
438 polygon.push_back((*pVertices)[pSectionIndex 184 polygon.push_back((*pVertices)[pSectionIndex]);
439 polygon.push_back((*pVertices)[pSectionIndex 185 polygon.push_back((*pVertices)[pSectionIndex+1]);
440 polygon.push_back((*pVertices)[pSectionIndex 186 polygon.push_back((*pVertices)[pSectionIndex+2]);
441 polygon.push_back((*pVertices)[pSectionIndex 187 polygon.push_back((*pVertices)[pSectionIndex+3]);
>> 188 // G4cout<<"ClipCrossSection: 0-1-2-3"<<G4endl;
442 CalculateClippedPolygonExtent(polygon,pVoxel 189 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
443 return; 190 return;
444 } 191 }
445 192
446 ////////////////////////////////////////////// 193 //////////////////////////////////////////////////////////////////////////////////
447 // 194 //
448 // Calculate the maximum and minimum extents o 195 // Calculate the maximum and minimum extents of the polygons
449 // joining the CrossSections at pSectionIndex- 196 // joining the CrossSections at pSectionIndex->pSectionIndex+3 and
450 // pSectionIndex+ 197 // pSectionIndex+4->pSectionIndex7
451 // 198 //
452 // in the List pVertices, within the boundarie 199 // in the List pVertices, within the boundaries of the voxel limits pVoxelLimit
453 // 200 //
454 // If the minimum is <pMin pMin is set to the 201 // If the minimum is <pMin pMin is set to the new minimum
455 // If the maximum is >pMax pMax is set to the 202 // If the maximum is >pMax pMax is set to the new maximum
456 // 203 //
457 // No modifications are made to pVertices 204 // No modifications are made to pVertices
458 205
459 void G4VSolid::ClipBetweenSections( G4Thr 206 void G4VSolid::ClipBetweenSections( G4ThreeVectorList* pVertices,
460 const G4int 207 const G4int pSectionIndex,
461 const G4Vox 208 const G4VoxelLimits& pVoxelLimit,
462 const EAxis << 209 const EAxis pAxis,
463 G4dou 210 G4double& pMin, G4double& pMax) const
464 { 211 {
465 G4ThreeVectorList polygon; 212 G4ThreeVectorList polygon;
466 polygon.reserve(4); 213 polygon.reserve(4);
467 polygon.push_back((*pVertices)[pSectionIndex 214 polygon.push_back((*pVertices)[pSectionIndex]);
468 polygon.push_back((*pVertices)[pSectionIndex 215 polygon.push_back((*pVertices)[pSectionIndex+4]);
469 polygon.push_back((*pVertices)[pSectionIndex 216 polygon.push_back((*pVertices)[pSectionIndex+5]);
470 polygon.push_back((*pVertices)[pSectionIndex 217 polygon.push_back((*pVertices)[pSectionIndex+1]);
>> 218 // G4cout<<"ClipBetweenSections: 0-4-5-1"<<G4endl;
471 CalculateClippedPolygonExtent(polygon,pVoxel 219 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
472 polygon.clear(); 220 polygon.clear();
473 221
474 polygon.push_back((*pVertices)[pSectionIndex 222 polygon.push_back((*pVertices)[pSectionIndex+1]);
475 polygon.push_back((*pVertices)[pSectionIndex 223 polygon.push_back((*pVertices)[pSectionIndex+5]);
476 polygon.push_back((*pVertices)[pSectionIndex 224 polygon.push_back((*pVertices)[pSectionIndex+6]);
477 polygon.push_back((*pVertices)[pSectionIndex 225 polygon.push_back((*pVertices)[pSectionIndex+2]);
>> 226 // G4cout<<"ClipBetweenSections: 1-5-6-2"<<G4endl;
478 CalculateClippedPolygonExtent(polygon,pVoxel 227 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
479 polygon.clear(); 228 polygon.clear();
480 229
481 polygon.push_back((*pVertices)[pSectionIndex 230 polygon.push_back((*pVertices)[pSectionIndex+2]);
482 polygon.push_back((*pVertices)[pSectionIndex 231 polygon.push_back((*pVertices)[pSectionIndex+6]);
483 polygon.push_back((*pVertices)[pSectionIndex 232 polygon.push_back((*pVertices)[pSectionIndex+7]);
484 polygon.push_back((*pVertices)[pSectionIndex 233 polygon.push_back((*pVertices)[pSectionIndex+3]);
>> 234 // G4cout<<"ClipBetweenSections: 2-6-7-3"<<G4endl;
485 CalculateClippedPolygonExtent(polygon,pVoxel 235 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
486 polygon.clear(); 236 polygon.clear();
487 237
488 polygon.push_back((*pVertices)[pSectionIndex 238 polygon.push_back((*pVertices)[pSectionIndex+3]);
489 polygon.push_back((*pVertices)[pSectionIndex 239 polygon.push_back((*pVertices)[pSectionIndex+7]);
490 polygon.push_back((*pVertices)[pSectionIndex 240 polygon.push_back((*pVertices)[pSectionIndex+4]);
491 polygon.push_back((*pVertices)[pSectionIndex 241 polygon.push_back((*pVertices)[pSectionIndex]);
>> 242 // G4cout<<"ClipBetweenSections: 3-7-4-0"<<G4endl;
492 CalculateClippedPolygonExtent(polygon,pVoxel 243 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
493 return; 244 return;
494 } 245 }
495 246
496 247
497 ////////////////////////////////////////////// 248 ///////////////////////////////////////////////////////////////////////////////
498 // 249 //
499 // Calculate the maximum and minimum extents o 250 // Calculate the maximum and minimum extents of the convex polygon pPolygon
500 // along the axis pAxis, within the limits pVo 251 // along the axis pAxis, within the limits pVoxelLimit
501 // 252 //
502 253
503 void 254 void
504 G4VSolid::CalculateClippedPolygonExtent(G4Thre 255 G4VSolid::CalculateClippedPolygonExtent(G4ThreeVectorList& pPolygon,
505 const G4Voxe 256 const G4VoxelLimits& pVoxelLimit,
506 const EAxis << 257 const EAxis pAxis,
507 G4doub 258 G4double& pMin,
508 G4doub 259 G4double& pMax) const
509 { 260 {
510 G4int noLeft,i; 261 G4int noLeft,i;
511 G4double component; 262 G4double component;
512 << 263 /*
>> 264 G4cout<<G4endl;
>> 265 for(i = 0 ; i < pPolygon.size() ; i++ )
>> 266 {
>> 267 G4cout << i << "\t"
>> 268 << "p.x = " << pPolygon[i].operator()(pAxis) << "\t"
>> 269 // << "p.y = " << pPolygon[i].y() << "\t"
>> 270 // << "p.z = " << pPolygon[i].z() << "\t"
>> 271 << G4endl;
>> 272 }
>> 273 G4cout<<G4endl;
>> 274 */
513 ClipPolygon(pPolygon,pVoxelLimit,pAxis); 275 ClipPolygon(pPolygon,pVoxelLimit,pAxis);
514 noLeft = (G4int)pPolygon.size(); << 276 noLeft = pPolygon.size();
515 277
516 if ( noLeft != 0 ) << 278 if ( noLeft )
517 { 279 {
518 for (i=0; i<noLeft; ++i) << 280 // G4cout<<G4endl;
>> 281 for (i=0;i<noLeft;i++)
519 { 282 {
520 component = pPolygon[i].operator()(pAxis 283 component = pPolygon[i].operator()(pAxis);
521 << 284 // G4cout <<i<<"\t"<<component<<G4endl;
522 if (component < pMin) << 285
523 { << 286 if (component < pMin)
524 pMin = component; << 287 {
>> 288 // G4cout <<i<<"\t"<<"Pmin = "<<component<<G4endl;
>> 289 pMin = component;
525 } 290 }
526 if (component > pMax) 291 if (component > pMax)
527 { << 292 {
528 pMax = component; << 293 // G4cout <<i<<"\t"<<"PMax = "<<component<<G4endl;
529 } << 294 pMax = component;
>> 295 }
530 } 296 }
>> 297 // G4cout<<G4endl;
531 } 298 }
>> 299 // G4cout<<"pMin = "<<pMin<<"\t"<<"pMax = "<<pMax<<G4endl;
532 } 300 }
533 301
534 ////////////////////////////////////////////// 302 /////////////////////////////////////////////////////////////////////////////
535 // 303 //
536 // Clip the convex polygon described by the ve 304 // Clip the convex polygon described by the vertices at
537 // pSectionIndex ->pSectionIndex+3 within pVer 305 // pSectionIndex ->pSectionIndex+3 within pVertices to the limits pVoxelLimit
538 // 306 //
539 // Set pMin to the smallest 307 // Set pMin to the smallest
540 // 308 //
541 // Calculate the extent of the polygon along p 309 // Calculate the extent of the polygon along pAxis, when clipped to the
542 // limits pVoxelLimit. If the polygon exists a 310 // limits pVoxelLimit. If the polygon exists after clippin, set pMin to
543 // the polygon's minimum extent along the axis 311 // the polygon's minimum extent along the axis if <pMin, and set pMax to
544 // the polygon's maximum extent along the axis 312 // the polygon's maximum extent along the axis if >pMax.
545 // 313 //
546 // The polygon is described by a set of vector 314 // The polygon is described by a set of vectors, where each vector represents
547 // a vertex, so that the polygon is described 315 // a vertex, so that the polygon is described by the vertex sequence:
548 // 0th->1st 1st->2nd 2nd->... nth->0th 316 // 0th->1st 1st->2nd 2nd->... nth->0th
549 // 317 //
550 // Modifications to the polygon are made 318 // Modifications to the polygon are made
551 // 319 //
552 // NOTE: Execessive copying during clipping 320 // NOTE: Execessive copying during clipping
553 321
554 void G4VSolid::ClipPolygon( G4ThreeVector 322 void G4VSolid::ClipPolygon( G4ThreeVectorList& pPolygon,
555 const G4VoxelLimits 323 const G4VoxelLimits& pVoxelLimit,
556 const EAxis 324 const EAxis ) const
557 { 325 {
558 G4ThreeVectorList outputPolygon; 326 G4ThreeVectorList outputPolygon;
559 327
560 if ( pVoxelLimit.IsLimited() ) 328 if ( pVoxelLimit.IsLimited() )
561 { 329 {
562 if (pVoxelLimit.IsXLimited() ) // && pAxis 330 if (pVoxelLimit.IsXLimited() ) // && pAxis != kXAxis)
563 { 331 {
564 G4VoxelLimits simpleLimit1; 332 G4VoxelLimits simpleLimit1;
565 simpleLimit1.AddLimit(kXAxis,pVoxelLimit 333 simpleLimit1.AddLimit(kXAxis,pVoxelLimit.GetMinXExtent(),kInfinity);
>> 334 // G4cout<<"MinXExtent()"<<G4endl;
566 ClipPolygonToSimpleLimits(pPolygon,outpu 335 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
567 << 336
568 pPolygon.clear(); 337 pPolygon.clear();
569 338
570 if ( outputPolygon.empty() ) return; << 339 if ( !outputPolygon.size() ) return;
571 340
572 G4VoxelLimits simpleLimit2; 341 G4VoxelLimits simpleLimit2;
>> 342 // G4cout<<"MaxXExtent()"<<G4endl;
573 simpleLimit2.AddLimit(kXAxis,-kInfinity, 343 simpleLimit2.AddLimit(kXAxis,-kInfinity,pVoxelLimit.GetMaxXExtent());
574 ClipPolygonToSimpleLimits(outputPolygon, 344 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
575 345
576 if ( pPolygon.empty() ) return; << 346 if ( !pPolygon.size() ) return;
577 else outputPoly 347 else outputPolygon.clear();
578 } 348 }
579 if ( pVoxelLimit.IsYLimited() ) // && pAxi 349 if ( pVoxelLimit.IsYLimited() ) // && pAxis != kYAxis)
580 { 350 {
581 G4VoxelLimits simpleLimit1; 351 G4VoxelLimits simpleLimit1;
582 simpleLimit1.AddLimit(kYAxis,pVoxelLimit 352 simpleLimit1.AddLimit(kYAxis,pVoxelLimit.GetMinYExtent(),kInfinity);
583 ClipPolygonToSimpleLimits(pPolygon,outpu 353 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
584 354
585 // Must always clear pPolygon - for clip 355 // Must always clear pPolygon - for clip to simpleLimit2 and in case of
586 // early exit 356 // early exit
587 357
588 pPolygon.clear(); 358 pPolygon.clear();
589 359
590 if ( outputPolygon.empty() ) return; << 360 if ( !outputPolygon.size() ) return;
591 361
592 G4VoxelLimits simpleLimit2; 362 G4VoxelLimits simpleLimit2;
593 simpleLimit2.AddLimit(kYAxis,-kInfinity, 363 simpleLimit2.AddLimit(kYAxis,-kInfinity,pVoxelLimit.GetMaxYExtent());
594 ClipPolygonToSimpleLimits(outputPolygon, 364 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
595 365
596 if ( pPolygon.empty() ) return; << 366 if ( !pPolygon.size() ) return;
597 else outputPoly 367 else outputPolygon.clear();
598 } 368 }
599 if ( pVoxelLimit.IsZLimited() ) // && pAxi 369 if ( pVoxelLimit.IsZLimited() ) // && pAxis != kZAxis)
600 { 370 {
601 G4VoxelLimits simpleLimit1; 371 G4VoxelLimits simpleLimit1;
602 simpleLimit1.AddLimit(kZAxis,pVoxelLimit 372 simpleLimit1.AddLimit(kZAxis,pVoxelLimit.GetMinZExtent(),kInfinity);
603 ClipPolygonToSimpleLimits(pPolygon,outpu 373 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
604 374
605 // Must always clear pPolygon - for clip 375 // Must always clear pPolygon - for clip to simpleLimit2 and in case of
606 // early exit 376 // early exit
607 377
608 pPolygon.clear(); 378 pPolygon.clear();
609 379
610 if ( outputPolygon.empty() ) return; << 380 if ( !outputPolygon.size() ) return;
611 381
612 G4VoxelLimits simpleLimit2; 382 G4VoxelLimits simpleLimit2;
613 simpleLimit2.AddLimit(kZAxis,-kInfinity, 383 simpleLimit2.AddLimit(kZAxis,-kInfinity,pVoxelLimit.GetMaxZExtent());
614 ClipPolygonToSimpleLimits(outputPolygon, 384 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
615 385
616 // Return after final clip - no cleanup 386 // Return after final clip - no cleanup
617 } 387 }
618 } 388 }
619 } 389 }
620 390
621 ////////////////////////////////////////////// 391 ////////////////////////////////////////////////////////////////////////////
622 // 392 //
623 // pVoxelLimits must be only limited along one 393 // pVoxelLimits must be only limited along one axis, and either the maximum
624 // along the axis must be +kInfinity, or the m 394 // along the axis must be +kInfinity, or the minimum -kInfinity
625 395
626 void 396 void
627 G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVe 397 G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVectorList& pPolygon,
628 G4ThreeVe 398 G4ThreeVectorList& outputPolygon,
629 const G4VoxelLi 399 const G4VoxelLimits& pVoxelLimit ) const
630 { 400 {
631 G4int i; 401 G4int i;
632 auto noVertices = (G4int)pPolygon.size(); << 402 G4int noVertices=pPolygon.size();
633 G4ThreeVector vEnd,vStart; 403 G4ThreeVector vEnd,vStart;
634 404
635 for (i = 0 ; i < noVertices ; ++i ) << 405 for (i = 0 ; i < noVertices ; i++ )
636 { 406 {
637 vStart = pPolygon[i]; 407 vStart = pPolygon[i];
>> 408 // G4cout << "i = " << i << G4endl;
638 if ( i == noVertices-1 ) vEnd = pPolygo 409 if ( i == noVertices-1 ) vEnd = pPolygon[0];
639 else vEnd = pPolygo 410 else vEnd = pPolygon[i+1];
640 411
641 if ( pVoxelLimit.Inside(vStart) ) 412 if ( pVoxelLimit.Inside(vStart) )
642 { 413 {
643 if (pVoxelLimit.Inside(vEnd)) 414 if (pVoxelLimit.Inside(vEnd))
644 { 415 {
645 // vStart and vEnd inside -> output en 416 // vStart and vEnd inside -> output end point
646 // 417 //
647 outputPolygon.push_back(vEnd); 418 outputPolygon.push_back(vEnd);
648 } 419 }
649 else 420 else
650 { 421 {
651 // vStart inside, vEnd outside -> outp 422 // vStart inside, vEnd outside -> output crossing point
652 // 423 //
>> 424 // G4cout << "vStart inside, vEnd outside" << G4endl;
653 pVoxelLimit.ClipToLimits(vStart,vEnd); 425 pVoxelLimit.ClipToLimits(vStart,vEnd);
654 outputPolygon.push_back(vEnd); 426 outputPolygon.push_back(vEnd);
655 } << 427 }
656 } 428 }
657 else 429 else
658 { 430 {
659 if (pVoxelLimit.Inside(vEnd)) 431 if (pVoxelLimit.Inside(vEnd))
660 { 432 {
661 // vStart outside, vEnd inside -> outp 433 // vStart outside, vEnd inside -> output inside section
662 // 434 //
>> 435 // G4cout << "vStart outside, vEnd inside" << G4endl;
663 pVoxelLimit.ClipToLimits(vStart,vEnd); 436 pVoxelLimit.ClipToLimits(vStart,vEnd);
664 outputPolygon.push_back(vStart); 437 outputPolygon.push_back(vStart);
665 outputPolygon.push_back(vEnd); << 438 outputPolygon.push_back(vEnd);
666 } 439 }
667 else // Both point outside -> no output 440 else // Both point outside -> no output
668 { 441 {
669 // outputPolygon.push_back(vStart); 442 // outputPolygon.push_back(vStart);
670 // outputPolygon.push_back(vEnd); << 443 // outputPolygon.push_back(vEnd);
671 } 444 }
672 } 445 }
673 } 446 }
674 } 447 }
675 448
676 ////////////////////////////////////////////// << 449 const G4VSolid* G4VSolid::GetConstituentSolid(G4int) const
677 // << 450 { return 0; }
678 // Throw exception (warning) for solids not im <<
679 451
680 void G4VSolid::BoundingLimits(G4ThreeVector& p << 452 G4VSolid* G4VSolid::GetConstituentSolid(G4int)
681 { << 453 { return 0; }
682 std::ostringstream message; <<
683 message << "Not implemented for solid: " <<
684 << GetEntityType() << " !" <<
685 << "\nReturning infinite boundinx bo <<
686 G4Exception("G4VSolid::BoundingLimits()", "G <<
687 JustWarning, message); <<
688 454
689 pMin.set(-kInfinity,-kInfinity,-kInfinity); << 455 const G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() const
690 pMax.set( kInfinity, kInfinity, kInfinity); << 456 { return 0; }
691 } <<
692 457
693 ////////////////////////////////////////////// << 458 G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr()
694 // << 459 { return 0; }
695 // Get G4VisExtent - bounding box for graphics <<
696 460
697 G4VisExtent G4VSolid::GetExtent () const << 461 G4VisExtent G4VSolid::GetExtent () const
698 { 462 {
699 G4VisExtent extent; 463 G4VisExtent extent;
700 G4VoxelLimits voxelLimits; // Defaults to " 464 G4VoxelLimits voxelLimits; // Defaults to "infinite" limits.
701 G4AffineTransform affineTransform; 465 G4AffineTransform affineTransform;
702 G4double vmin, vmax; 466 G4double vmin, vmax;
703 CalculateExtent(kXAxis,voxelLimits,affineTra 467 CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
704 extent.SetXmin (vmin); 468 extent.SetXmin (vmin);
705 extent.SetXmax (vmax); 469 extent.SetXmax (vmax);
706 CalculateExtent(kYAxis,voxelLimits,affineTra 470 CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
707 extent.SetYmin (vmin); 471 extent.SetYmin (vmin);
708 extent.SetYmax (vmax); 472 extent.SetYmax (vmax);
709 CalculateExtent(kZAxis,voxelLimits,affineTra 473 CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
710 extent.SetZmin (vmin); 474 extent.SetZmin (vmin);
711 extent.SetZmax (vmax); 475 extent.SetZmax (vmax);
712 return extent; 476 return extent;
713 } 477 }
714 478
715 G4Polyhedron* G4VSolid::CreatePolyhedron () co 479 G4Polyhedron* G4VSolid::CreatePolyhedron () const
716 { 480 {
717 return nullptr; << 481 return 0;
>> 482 }
>> 483
>> 484 G4NURBS* G4VSolid::CreateNURBS () const
>> 485 {
>> 486 return 0;
718 } 487 }
719 488
720 G4Polyhedron* G4VSolid::GetPolyhedron () const 489 G4Polyhedron* G4VSolid::GetPolyhedron () const
721 { 490 {
722 return nullptr; << 491 return 0;
>> 492 }
>> 493
>> 494 ////////////////////////////////////////////////////////////////
>> 495 //
>> 496 // Returns an estimation of the solid volume in internal units.
>> 497 // The number of statistics and error accuracy is fixed.
>> 498 // This method may be overloaded by derived classes to compute the
>> 499 // exact geometrical quantity for solids where this is possible.
>> 500 // or anyway to cache the computed value.
>> 501 // This implementation does NOT cache the computed value.
>> 502
>> 503 G4double G4VSolid::GetCubicVolume()
>> 504 {
>> 505 G4int cubVolStatistics = 1000000;
>> 506 G4double cubVolEpsilon = 0.001;
>> 507 return EstimateCubicVolume(cubVolStatistics, cubVolEpsilon);
>> 508 }
>> 509
>> 510 ////////////////////////////////////////////////////////////////
>> 511 //
>> 512 // Calculate cubic volume based on Inside() method.
>> 513 // Accuracy is limited by the second argument or the statistics
>> 514 // expressed by the first argument.
>> 515 // Implementation is courtesy of Vasiliki Despoina Mitsou,
>> 516 // University of Athens.
>> 517
>> 518 G4double G4VSolid::EstimateCubicVolume(G4int nStat, G4double epsilon) const
>> 519 {
>> 520 G4int iInside=0;
>> 521 G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,volume;
>> 522 G4bool yesno;
>> 523 G4ThreeVector p;
>> 524 EInside in;
>> 525
>> 526 // values needed for CalculateExtent signature
>> 527
>> 528 G4VoxelLimits limit; // Unlimited
>> 529 G4AffineTransform origin;
>> 530
>> 531 // min max extents of pSolid along X,Y,Z
>> 532
>> 533 yesno = this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
>> 534 yesno = this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
>> 535 yesno = this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
>> 536
>> 537 // limits
>> 538
>> 539 if(nStat < 100) nStat = 100;
>> 540 if(epsilon > 0.01) epsilon = 0.01;
>> 541
>> 542 for(G4int i = 0; i < nStat; i++ )
>> 543 {
>> 544 px = minX+(maxX-minX)*G4UniformRand();
>> 545 py = minY+(maxY-minY)*G4UniformRand();
>> 546 pz = minZ+(maxZ-minZ)*G4UniformRand();
>> 547 p = G4ThreeVector(px,py,pz);
>> 548 in = this->Inside(p);
>> 549 if(in != kOutside) iInside++;
>> 550 }
>> 551 volume = (maxX-minX)*(maxY-minY)*(maxZ-minZ)*iInside/nStat;
>> 552 return volume;
>> 553 }
>> 554
>> 555 ////////////////////////////////////////////////////////////////
>> 556 //
>> 557 // Returns an estimation of the solid surface area in internal units.
>> 558 // The number of statistics and error accuracy is fixed.
>> 559 // This method may be overloaded by derived classes to compute the
>> 560 // exact geometrical quantity for solids where this is possible.
>> 561 // or anyway to cache the computed value.
>> 562 // This implementation does NOT cache the computed value.
>> 563
>> 564 G4double G4VSolid::GetSurfaceArea()
>> 565 {
>> 566 G4int stat = 1000000;
>> 567 G4double ell = -1.;
>> 568 return EstimateSurfaceArea(stat,ell);
>> 569 }
>> 570
>> 571 ////////////////////////////////////////////////////////////////
>> 572 //
>> 573 // Estimate surface area based on Inside(), DistanceToIn(), and
>> 574 // DistanceToOut() methods. Accuracy is limited by the statistics
>> 575 // defined by the first argument. Implemented by Mikhail Kosov.
>> 576
>> 577 G4double G4VSolid::EstimateSurfaceArea(G4int nStat, G4double ell) const
>> 578 {
>> 579 G4int inside=0;
>> 580 G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,surf;
>> 581 G4bool yesno;
>> 582 G4ThreeVector p;
>> 583 EInside in;
>> 584
>> 585 // values needed for CalculateExtent signature
>> 586
>> 587 G4VoxelLimits limit; // Unlimited
>> 588 G4AffineTransform origin;
>> 589
>> 590 // min max extents of pSolid along X,Y,Z
>> 591
>> 592 yesno = this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
>> 593 yesno = this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
>> 594 yesno = this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
>> 595
>> 596 // limits
>> 597
>> 598 if(nStat < 100) { nStat = 100; }
>> 599
>> 600 G4double dX=maxX-minX;
>> 601 G4double dY=maxY-minY;
>> 602 G4double dZ=maxZ-minZ;
>> 603 if(ell<=0.) // Automatic definition of skin thickness
>> 604 {
>> 605 G4double minval=dX;
>> 606 if(dY<dX) { minval=dY; }
>> 607 if(dZ<minval) { minval=dZ; }
>> 608 ell=.01*minval;
>> 609 }
>> 610
>> 611 G4double dd=2*ell;
>> 612 minX-=ell; minY-=ell; minZ-=ell; dX+=dd; dY+=dd; dZ+=dd;
>> 613
>> 614 for(G4int i = 0; i < nStat; i++ )
>> 615 {
>> 616 px = minX+dX*G4UniformRand();
>> 617 py = minY+dY*G4UniformRand();
>> 618 pz = minZ+dZ*G4UniformRand();
>> 619 p = G4ThreeVector(px,py,pz);
>> 620 in = this->Inside(p);
>> 621 if(in != kOutside)
>> 622 {
>> 623 if (DistanceToOut(p)<ell) { inside++; }
>> 624 }
>> 625 else if(DistanceToIn(p)<ell) { inside++; }
>> 626 }
>> 627 // @@ The conformal correction can be upgraded
>> 628 surf = dX*dY*dZ*inside/dd/nStat;
>> 629 return surf;
723 } 630 }
724 631