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