| 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 // G4ExtrudedSolid implementation 26 // G4ExtrudedSolid implementation
27 // 27 //
28 // Author: Ivana Hrivnacova, IPN Orsay 28 // Author: Ivana Hrivnacova, IPN Orsay
29 // 29 //
30 // CHANGE HISTORY 30 // CHANGE HISTORY
31 // -------------- 31 // --------------
32 // 32 //
33 // 31.10.2017 E.Tcherniaev: added implementati 33 // 31.10.2017 E.Tcherniaev: added implementation for a non-convex
34 // right prism 34 // right prism
35 // 08.09.2017 E.Tcherniaev: added implementati 35 // 08.09.2017 E.Tcherniaev: added implementation for a convex
36 // right prism 36 // right prism
37 // 21.10.2016 E.Tcherniaev: reimplemented Calc 37 // 21.10.2016 E.Tcherniaev: reimplemented CalculateExtent(),
38 // used G4GeomTools::PolygonArea() 38 // used G4GeomTools::PolygonArea() to calculate area,
39 // replaced IsConvex() with G4GeomT 39 // replaced IsConvex() with G4GeomTools::IsConvex()
40 // 02.03.2016 E.Tcherniaev: added CheckPolygon 40 // 02.03.2016 E.Tcherniaev: added CheckPolygon() to remove
41 // collinear and coincident points 41 // collinear and coincident points from polygon
42 // ------------------------------------------- 42 // --------------------------------------------------------------------
43 43
44 #include "G4ExtrudedSolid.hh" 44 #include "G4ExtrudedSolid.hh"
45 45
46 #if !defined(G4GEOM_USE_UEXTRUDEDSOLID) 46 #if !defined(G4GEOM_USE_UEXTRUDEDSOLID)
47 47
48 #include <set> 48 #include <set>
49 #include <algorithm> 49 #include <algorithm>
50 #include <cmath> 50 #include <cmath>
51 #include <iomanip> 51 #include <iomanip>
52 52
53 #include "G4GeomTools.hh" 53 #include "G4GeomTools.hh"
54 #include "G4VoxelLimits.hh" 54 #include "G4VoxelLimits.hh"
55 #include "G4AffineTransform.hh" 55 #include "G4AffineTransform.hh"
56 #include "G4BoundingEnvelope.hh" 56 #include "G4BoundingEnvelope.hh"
57 57
58 #include "G4GeometryTolerance.hh" 58 #include "G4GeometryTolerance.hh"
59 #include "G4PhysicalConstants.hh" 59 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh" 60 #include "G4SystemOfUnits.hh"
61 #include "G4TriangularFacet.hh" 61 #include "G4TriangularFacet.hh"
62 #include "G4QuadrangularFacet.hh" 62 #include "G4QuadrangularFacet.hh"
63 63
64 //____________________________________________ 64 //_____________________________________________________________________________
65 65
66 G4ExtrudedSolid::G4ExtrudedSolid( const G4Stri 66 G4ExtrudedSolid::G4ExtrudedSolid( const G4String& pName,
67 const std::v 67 const std::vector<G4TwoVector>& polygon,
68 const std::v 68 const std::vector<ZSection>& zsections)
69 : G4TessellatedSolid(pName), 69 : G4TessellatedSolid(pName),
70 fNv(polygon.size()), 70 fNv(polygon.size()),
71 fNz(zsections.size()), 71 fNz(zsections.size()),
72 fIsConvex(false), 72 fIsConvex(false),
73 fGeometryType("G4ExtrudedSolid"), 73 fGeometryType("G4ExtrudedSolid"),
74 fSolidType(0) 74 fSolidType(0)
75 { 75 {
76 // General constructor 76 // General constructor
77 77
78 // First check input parameters 78 // First check input parameters
79 79
80 if (fNv < 3) 80 if (fNv < 3)
81 { 81 {
82 std::ostringstream message; 82 std::ostringstream message;
83 message << "Number of vertices in polygon 83 message << "Number of vertices in polygon < 3 - " << pName;
84 G4Exception("G4ExtrudedSolid::G4ExtrudedSo 84 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
85 FatalErrorInArgument, message) 85 FatalErrorInArgument, message);
86 } 86 }
87 87
88 if (fNz < 2) 88 if (fNz < 2)
89 { 89 {
90 std::ostringstream message; 90 std::ostringstream message;
91 message << "Number of z-sides < 2 - " << p 91 message << "Number of z-sides < 2 - " << pName;
92 G4Exception("G4ExtrudedSolid::G4ExtrudedSo 92 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
93 FatalErrorInArgument, message) 93 FatalErrorInArgument, message);
94 } 94 }
95 95
96 for ( std::size_t i=0; i<fNz-1; ++i ) 96 for ( std::size_t i=0; i<fNz-1; ++i )
97 { 97 {
98 if ( zsections[i].fZ > zsections[i+1].fZ ) 98 if ( zsections[i].fZ > zsections[i+1].fZ )
99 { 99 {
100 std::ostringstream message; 100 std::ostringstream message;
101 message << "Z-sections have to be ordere 101 message << "Z-sections have to be ordered by z value (z0 < z1 < z2...) - "
102 << pName; 102 << pName;
103 G4Exception("G4ExtrudedSolid::G4Extruded 103 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
104 FatalErrorInArgument, messag 104 FatalErrorInArgument, message);
105 } 105 }
106 if ( std::fabs( zsections[i+1].fZ - zsecti 106 if ( std::fabs( zsections[i+1].fZ - zsections[i].fZ ) < kCarToleranceHalf )
107 { 107 {
108 std::ostringstream message; 108 std::ostringstream message;
109 message << "Z-sections with the same z p 109 message << "Z-sections with the same z position are not supported - "
110 << pName; 110 << pName;
111 G4Exception("G4ExtrudedSolid::G4Extruded 111 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0001",
112 FatalException, message); 112 FatalException, message);
113 } 113 }
114 } 114 }
115 115
116 // Copy polygon 116 // Copy polygon
117 // 117 //
118 fPolygon = polygon; 118 fPolygon = polygon;
119 119
120 // Remove collinear and coincident vertices, 120 // Remove collinear and coincident vertices, if any
121 // 121 //
122 std::vector<G4int> removedVertices; 122 std::vector<G4int> removedVertices;
123 G4GeomTools::RemoveRedundantVertices(fPolygo 123 G4GeomTools::RemoveRedundantVertices(fPolygon,removedVertices,
124 2*kCarT 124 2*kCarTolerance);
125 if (!removedVertices.empty()) << 125 if (removedVertices.size() != 0)
126 { 126 {
127 std::size_t nremoved = removedVertices.siz 127 std::size_t nremoved = removedVertices.size();
128 std::ostringstream message; 128 std::ostringstream message;
129 message << "The following "<< nremoved 129 message << "The following "<< nremoved
130 << " vertices have been removed fr 130 << " vertices have been removed from polygon in " << pName
131 << "\nas collinear or coincident w 131 << "\nas collinear or coincident with other vertices: "
132 << removedVertices[0]; 132 << removedVertices[0];
133 for (std::size_t i=1; i<nremoved; ++i) mes 133 for (std::size_t i=1; i<nremoved; ++i) message << ", " << removedVertices[i];
134 G4Exception("G4ExtrudedSolid::G4ExtrudedSo 134 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001",
135 JustWarning, message); 135 JustWarning, message);
136 } 136 }
137 137
138 fNv = fPolygon.size(); 138 fNv = fPolygon.size();
139 if (fNv < 3) 139 if (fNv < 3)
140 { 140 {
141 std::ostringstream message; 141 std::ostringstream message;
142 message << "Number of vertices in polygon 142 message << "Number of vertices in polygon after removal < 3 - " << pName;
143 G4Exception("G4ExtrudedSolid::G4ExtrudedSo 143 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
144 FatalErrorInArgument, message) 144 FatalErrorInArgument, message);
145 } 145 }
146 146
147 // Check if polygon vertices are defined clo 147 // Check if polygon vertices are defined clockwise
148 // (the area is positive if polygon vertices 148 // (the area is positive if polygon vertices are defined anti-clockwise)
149 // 149 //
150 if (G4GeomTools::PolygonArea(fPolygon) > 0.) 150 if (G4GeomTools::PolygonArea(fPolygon) > 0.)
151 { 151 {
152 // Polygon vertices are defined anti-clock 152 // Polygon vertices are defined anti-clockwise, we revert them
153 // G4Exception("G4ExtrudedSolid::G4Extrude 153 // G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001",
154 // JustWarning, 154 // JustWarning,
155 // "Polygon vertices defined an 155 // "Polygon vertices defined anti-clockwise, reverting polygon");
156 std::reverse(fPolygon.begin(),fPolygon.end 156 std::reverse(fPolygon.begin(),fPolygon.end());
157 } 157 }
158 158
159 // Copy z-sections 159 // Copy z-sections
160 // 160 //
161 fZSections = zsections; 161 fZSections = zsections;
162 162
163 G4bool result = MakeFacets(); 163 G4bool result = MakeFacets();
164 if (!result) 164 if (!result)
165 { 165 {
166 std::ostringstream message; 166 std::ostringstream message;
167 message << "Making facets failed - " << pN 167 message << "Making facets failed - " << pName;
168 G4Exception("G4ExtrudedSolid::G4ExtrudedSo 168 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0003",
169 FatalException, message); 169 FatalException, message);
170 } 170 }
171 fIsConvex = G4GeomTools::IsConvex(fPolygon); 171 fIsConvex = G4GeomTools::IsConvex(fPolygon);
172 172
173 ComputeProjectionParameters(); 173 ComputeProjectionParameters();
174 174
175 // Check if the solid is a right prism, if s 175 // Check if the solid is a right prism, if so then set lateral planes
176 // 176 //
177 if ((fNz == 2) 177 if ((fNz == 2)
178 && (fZSections[0].fScale == 1) && (fZSec 178 && (fZSections[0].fScale == 1) && (fZSections[1].fScale == 1)
179 && (fZSections[0].fOffset == G4TwoVector 179 && (fZSections[0].fOffset == G4TwoVector(0,0))
180 && (fZSections[1].fOffset == G4TwoVector 180 && (fZSections[1].fOffset == G4TwoVector(0,0)))
181 { 181 {
182 fSolidType = (fIsConvex) ? 1 : 2; // 1 - c 182 fSolidType = (fIsConvex) ? 1 : 2; // 1 - convex, 2 - non-convex right prism
183 ComputeLateralPlanes(); 183 ComputeLateralPlanes();
184 } 184 }
185 } 185 }
186 186
187 //____________________________________________ 187 //_____________________________________________________________________________
188 188
189 G4ExtrudedSolid::G4ExtrudedSolid( const G4Stri 189 G4ExtrudedSolid::G4ExtrudedSolid( const G4String& pName,
190 const std::v 190 const std::vector<G4TwoVector>& polygon,
191 G4doub 191 G4double dz,
192 const G4TwoV 192 const G4TwoVector& off1, G4double scale1,
193 const G4TwoV 193 const G4TwoVector& off2, G4double scale2 )
194 : G4TessellatedSolid(pName), 194 : G4TessellatedSolid(pName),
195 fNv(polygon.size()), 195 fNv(polygon.size()),
196 fNz(2), 196 fNz(2),
197 fGeometryType("G4ExtrudedSolid") 197 fGeometryType("G4ExtrudedSolid")
198 { 198 {
199 // Special constructor for solid with 2 z-se 199 // Special constructor for solid with 2 z-sections
200 200
201 // First check input parameters 201 // First check input parameters
202 // 202 //
203 if (fNv < 3) 203 if (fNv < 3)
204 { 204 {
205 std::ostringstream message; 205 std::ostringstream message;
206 message << "Number of vertices in polygon 206 message << "Number of vertices in polygon < 3 - " << pName;
207 G4Exception("G4ExtrudedSolid::G4ExtrudedSo 207 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
208 FatalErrorInArgument, message) 208 FatalErrorInArgument, message);
209 } 209 }
210 210
211 // Copy polygon 211 // Copy polygon
212 // 212 //
213 fPolygon = polygon; 213 fPolygon = polygon;
214 214
215 // Remove collinear and coincident vertices, 215 // Remove collinear and coincident vertices, if any
216 // 216 //
217 std::vector<G4int> removedVertices; 217 std::vector<G4int> removedVertices;
218 G4GeomTools::RemoveRedundantVertices(fPolygo 218 G4GeomTools::RemoveRedundantVertices(fPolygon,removedVertices,
219 2*kCarT 219 2*kCarTolerance);
220 if (!removedVertices.empty()) << 220 if (removedVertices.size() != 0)
221 { 221 {
222 std::size_t nremoved = removedVertices.siz 222 std::size_t nremoved = removedVertices.size();
223 std::ostringstream message; 223 std::ostringstream message;
224 message << "The following "<< nremoved 224 message << "The following "<< nremoved
225 << " vertices have been removed fr 225 << " vertices have been removed from polygon in " << pName
226 << "\nas collinear or coincident w 226 << "\nas collinear or coincident with other vertices: "
227 << removedVertices[0]; 227 << removedVertices[0];
228 for (std::size_t i=1; i<nremoved; ++i) mes 228 for (std::size_t i=1; i<nremoved; ++i) message << ", " << removedVertices[i];
229 G4Exception("G4ExtrudedSolid::G4ExtrudedSo 229 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001",
230 JustWarning, message); 230 JustWarning, message);
231 } 231 }
232 232
233 fNv = fPolygon.size(); 233 fNv = fPolygon.size();
234 if (fNv < 3) 234 if (fNv < 3)
235 { 235 {
236 std::ostringstream message; 236 std::ostringstream message;
237 message << "Number of vertices in polygon 237 message << "Number of vertices in polygon after removal < 3 - " << pName;
238 G4Exception("G4ExtrudedSolid::G4ExtrudedSo 238 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
239 FatalErrorInArgument, message) 239 FatalErrorInArgument, message);
240 } 240 }
241 241
242 // Check if polygon vertices are defined clo 242 // Check if polygon vertices are defined clockwise
243 // (the area is positive if polygon vertices 243 // (the area is positive if polygon vertices are defined anti-clockwise)
244 // 244 //
245 if (G4GeomTools::PolygonArea(fPolygon) > 0.) 245 if (G4GeomTools::PolygonArea(fPolygon) > 0.)
246 { 246 {
247 // Polygon vertices are defined anti-clock 247 // Polygon vertices are defined anti-clockwise, we revert them
248 // G4Exception("G4ExtrudedSolid::G4Extrude 248 // G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001",
249 // JustWarning, 249 // JustWarning,
250 // "Polygon vertices defined an 250 // "Polygon vertices defined anti-clockwise, reverting polygon");
251 std::reverse(fPolygon.begin(),fPolygon.end 251 std::reverse(fPolygon.begin(),fPolygon.end());
252 } 252 }
253 253
254 // Copy z-sections 254 // Copy z-sections
255 // 255 //
256 fZSections.emplace_back(-dz, off1, scale1); << 256 fZSections.push_back(ZSection(-dz, off1, scale1));
257 fZSections.emplace_back( dz, off2, scale2); << 257 fZSections.push_back(ZSection( dz, off2, scale2));
258 258
259 G4bool result = MakeFacets(); 259 G4bool result = MakeFacets();
260 if (!result) 260 if (!result)
261 { 261 {
262 std::ostringstream message; 262 std::ostringstream message;
263 message << "Making facets failed - " << pN 263 message << "Making facets failed - " << pName;
264 G4Exception("G4ExtrudedSolid::G4ExtrudedSo 264 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0003",
265 FatalException, message); 265 FatalException, message);
266 } 266 }
267 fIsConvex = G4GeomTools::IsConvex(fPolygon); 267 fIsConvex = G4GeomTools::IsConvex(fPolygon);
268 268
269 ComputeProjectionParameters(); 269 ComputeProjectionParameters();
270 270
271 // Check if the solid is a right prism, if s 271 // Check if the solid is a right prism, if so then set lateral planes
272 // 272 //
273 if ((scale1 == 1) && (scale2 == 1) 273 if ((scale1 == 1) && (scale2 == 1)
274 && (off1 == G4TwoVector(0,0)) && (off2 = 274 && (off1 == G4TwoVector(0,0)) && (off2 == G4TwoVector(0,0)))
275 { 275 {
276 fSolidType = (fIsConvex) ? 1 : 2; // 1 - c 276 fSolidType = (fIsConvex) ? 1 : 2; // 1 - convex, 2 - non-convex right prism
277 ComputeLateralPlanes(); 277 ComputeLateralPlanes();
278 } 278 }
279 } 279 }
280 280
281 //____________________________________________ 281 //_____________________________________________________________________________
282 282
283 G4ExtrudedSolid::G4ExtrudedSolid( __void__& a 283 G4ExtrudedSolid::G4ExtrudedSolid( __void__& a )
284 : G4TessellatedSolid(a), fGeometryType("G4Ex << 284 : G4TessellatedSolid(a), fNv(0), fNz(0),
>> 285 fGeometryType("G4ExtrudedSolid")
285 { 286 {
286 // Fake default constructor - sets only memb 287 // Fake default constructor - sets only member data and allocates memory
287 // for usage rest 288 // for usage restricted to object persistency.
288 } 289 }
289 290
290 //____________________________________________ 291 //_____________________________________________________________________________
291 292
292 G4ExtrudedSolid::G4ExtrudedSolid(const G4Extru << 293 G4ExtrudedSolid::G4ExtrudedSolid(const G4ExtrudedSolid& rhs)
>> 294 : G4TessellatedSolid(rhs), fNv(rhs.fNv), fNz(rhs.fNz),
>> 295 fPolygon(rhs.fPolygon), fZSections(rhs.fZSections),
>> 296 fTriangles(rhs.fTriangles), fIsConvex(rhs.fIsConvex),
>> 297 fGeometryType(rhs.fGeometryType),
>> 298 fSolidType(rhs.fSolidType), fPlanes(rhs.fPlanes),
>> 299 fLines(rhs.fLines), fLengths(rhs.fLengths),
>> 300 fKScales(rhs.fKScales), fScale0s(rhs.fScale0s),
>> 301 fKOffsets(rhs.fKOffsets), fOffset0s(rhs.fOffset0s)
>> 302 {
>> 303 }
293 304
294 //____________________________________________ 305 //_____________________________________________________________________________
295 306
296 G4ExtrudedSolid& G4ExtrudedSolid::operator = ( 307 G4ExtrudedSolid& G4ExtrudedSolid::operator = (const G4ExtrudedSolid& rhs)
297 { 308 {
298 // Check assignment to self 309 // Check assignment to self
299 // 310 //
300 if (this == &rhs) { return *this; } 311 if (this == &rhs) { return *this; }
301 312
302 // Copy base class data 313 // Copy base class data
303 // 314 //
304 G4TessellatedSolid::operator=(rhs); 315 G4TessellatedSolid::operator=(rhs);
305 316
306 // Copy data 317 // Copy data
307 // 318 //
308 fNv = rhs.fNv; fNz = rhs.fNz; 319 fNv = rhs.fNv; fNz = rhs.fNz;
309 fPolygon = rhs.fPolygon; fZSections = rhs.f 320 fPolygon = rhs.fPolygon; fZSections = rhs.fZSections;
310 fTriangles = rhs.fTriangles; fIsConvex = rh 321 fTriangles = rhs.fTriangles; fIsConvex = rhs.fIsConvex;
311 fGeometryType = rhs.fGeometryType; 322 fGeometryType = rhs.fGeometryType;
312 fSolidType = rhs.fSolidType; fPlanes = rhs. 323 fSolidType = rhs.fSolidType; fPlanes = rhs.fPlanes;
313 fLines = rhs.fLines; fLengths = rhs.fLength 324 fLines = rhs.fLines; fLengths = rhs.fLengths;
314 fKScales = rhs.fKScales; fScale0s = rhs.fSc 325 fKScales = rhs.fKScales; fScale0s = rhs.fScale0s;
315 fKOffsets = rhs.fKOffsets; fOffset0s = rhs. 326 fKOffsets = rhs.fKOffsets; fOffset0s = rhs.fOffset0s;
316 327
317 return *this; 328 return *this;
318 } 329 }
319 330
320 //____________________________________________ 331 //_____________________________________________________________________________
321 332
322 G4ExtrudedSolid::~G4ExtrudedSolid() 333 G4ExtrudedSolid::~G4ExtrudedSolid()
323 { 334 {
324 // Destructor 335 // Destructor
325 } 336 }
326 337
327 //____________________________________________ 338 //_____________________________________________________________________________
328 339
329 void G4ExtrudedSolid::ComputeProjectionParamet 340 void G4ExtrudedSolid::ComputeProjectionParameters()
330 { 341 {
331 // Compute parameters for point projections 342 // Compute parameters for point projections p(z)
332 // to the polygon scale & offset: 343 // to the polygon scale & offset:
333 // scale(z) = k*z + scale0 344 // scale(z) = k*z + scale0
334 // offset(z) = l*z + offset0 345 // offset(z) = l*z + offset0
335 // p(z) = scale(z)*p0 + offset(z) 346 // p(z) = scale(z)*p0 + offset(z)
336 // p0 = (p(z) - offset(z))/scale(z); 347 // p0 = (p(z) - offset(z))/scale(z);
337 // 348 //
338 349
339 for (std::size_t iz=0; iz<fNz-1; ++iz) 350 for (std::size_t iz=0; iz<fNz-1; ++iz)
340 { 351 {
341 G4double z1 = fZSections[iz].fZ; 352 G4double z1 = fZSections[iz].fZ;
342 G4double z2 = fZSections[iz+1].fZ; 353 G4double z2 = fZSections[iz+1].fZ;
343 G4double scale1 = fZSections[iz].fScale; 354 G4double scale1 = fZSections[iz].fScale;
344 G4double scale2 = fZSections[iz+1].fScale 355 G4double scale2 = fZSections[iz+1].fScale;
345 G4TwoVector off1 = fZSections[iz].fOffset; 356 G4TwoVector off1 = fZSections[iz].fOffset;
346 G4TwoVector off2 = fZSections[iz+1].fOffse 357 G4TwoVector off2 = fZSections[iz+1].fOffset;
347 358
348 G4double kscale = (scale2 - scale1)/(z2 - 359 G4double kscale = (scale2 - scale1)/(z2 - z1);
349 G4double scale0 = scale2 - kscale*(z2 - z 360 G4double scale0 = scale2 - kscale*(z2 - z1)/2.0;
350 G4TwoVector koff = (off2 - off1)/(z2 - z1) 361 G4TwoVector koff = (off2 - off1)/(z2 - z1);
351 G4TwoVector off0 = off2 - koff*(z2 - z1)/ 362 G4TwoVector off0 = off2 - koff*(z2 - z1)/2.0;
352 363
353 fKScales.push_back(kscale); 364 fKScales.push_back(kscale);
354 fScale0s.push_back(scale0); 365 fScale0s.push_back(scale0);
355 fKOffsets.push_back(koff); 366 fKOffsets.push_back(koff);
356 fOffset0s.push_back(off0); 367 fOffset0s.push_back(off0);
357 } 368 }
358 } 369 }
359 370
360 //____________________________________________ 371 //_____________________________________________________________________________
361 372
362 void G4ExtrudedSolid::ComputeLateralPlanes() 373 void G4ExtrudedSolid::ComputeLateralPlanes()
363 { 374 {
364 // Compute lateral planes: a*x + b*y + c*z + 375 // Compute lateral planes: a*x + b*y + c*z + d = 0
365 // 376 //
366 std::size_t Nv = fPolygon.size(); 377 std::size_t Nv = fPolygon.size();
367 fPlanes.resize(Nv); 378 fPlanes.resize(Nv);
368 for (std::size_t i=0, k=Nv-1; i<Nv; k=i++) 379 for (std::size_t i=0, k=Nv-1; i<Nv; k=i++)
369 { 380 {
370 G4TwoVector norm = (fPolygon[i] - fPolygon 381 G4TwoVector norm = (fPolygon[i] - fPolygon[k]).unit();
371 fPlanes[i].a = -norm.y(); 382 fPlanes[i].a = -norm.y();
372 fPlanes[i].b = norm.x(); 383 fPlanes[i].b = norm.x();
373 fPlanes[i].c = 0; 384 fPlanes[i].c = 0;
374 fPlanes[i].d = norm.y()*fPolygon[i].x() - 385 fPlanes[i].d = norm.y()*fPolygon[i].x() - norm.x()*fPolygon[i].y();
375 } 386 }
376 387
377 // Compute edge equations: x = k*y + m 388 // Compute edge equations: x = k*y + m
378 // and edge lengths 389 // and edge lengths
379 // 390 //
380 fLines.resize(Nv); 391 fLines.resize(Nv);
381 fLengths.resize(Nv); 392 fLengths.resize(Nv);
382 for (std::size_t i=0, k=Nv-1; i<Nv; k=i++) 393 for (std::size_t i=0, k=Nv-1; i<Nv; k=i++)
383 { 394 {
384 if (fPolygon[k].y() == fPolygon[i].y()) 395 if (fPolygon[k].y() == fPolygon[i].y())
385 { 396 {
386 fLines[i].k = 0; 397 fLines[i].k = 0;
387 fLines[i].m = fPolygon[i].x(); 398 fLines[i].m = fPolygon[i].x();
388 } 399 }
389 else 400 else
390 { 401 {
391 G4double ctg = (fPolygon[k].x()-fPolygon 402 G4double ctg = (fPolygon[k].x()-fPolygon[i].x())/(fPolygon[k].y()-fPolygon[i].y());
392 fLines[i].k = ctg; 403 fLines[i].k = ctg;
393 fLines[i].m = fPolygon[i].x() - ctg*fPol 404 fLines[i].m = fPolygon[i].x() - ctg*fPolygon[i].y();
394 } 405 }
395 fLengths[i] = (fPolygon[i] - fPolygon[k] 406 fLengths[i] = (fPolygon[i] - fPolygon[k]).mag();
396 } 407 }
397 } 408 }
398 409
399 //____________________________________________ 410 //_____________________________________________________________________________
400 411
401 G4ThreeVector G4ExtrudedSolid::GetVertex(G4int 412 G4ThreeVector G4ExtrudedSolid::GetVertex(G4int iz, G4int ind) const
402 { 413 {
403 // Shift and scale vertices 414 // Shift and scale vertices
404 415
405 return { fPolygon[ind].x() * fZSections[iz]. << 416 return G4ThreeVector( fPolygon[ind].x() * fZSections[iz].fScale
406 + fZSections[iz].fOffset.x(), << 417 + fZSections[iz].fOffset.x(),
407 fPolygon[ind].y() * fZSections[iz]. << 418 fPolygon[ind].y() * fZSections[iz].fScale
408 + fZSections[iz].fOffset.y(), << 419 + fZSections[iz].fOffset.y(), fZSections[iz].fZ);
409 fZSections[iz].fZ }; <<
410 } 420 }
411 421
412 //____________________________________________ 422 //_____________________________________________________________________________
413 423
414 G4TwoVector G4ExtrudedSolid::ProjectPoint(cons 424 G4TwoVector G4ExtrudedSolid::ProjectPoint(const G4ThreeVector& point) const
415 { 425 {
416 // Project point in the polygon scale 426 // Project point in the polygon scale
417 // scale(z) = k*z + scale0 427 // scale(z) = k*z + scale0
418 // offset(z) = l*z + offset0 428 // offset(z) = l*z + offset0
419 // p(z) = scale(z)*p0 + offset(z) 429 // p(z) = scale(z)*p0 + offset(z)
420 // p0 = (p(z) - offset(z))/scale(z); 430 // p0 = (p(z) - offset(z))/scale(z);
421 431
422 // Select projection (z-segment of the solid 432 // Select projection (z-segment of the solid) according to p.z()
423 // 433 //
424 std::size_t iz = 0; 434 std::size_t iz = 0;
425 while ( point.z() > fZSections[iz+1].fZ && i 435 while ( point.z() > fZSections[iz+1].fZ && iz < fNz-2 ) { ++iz; }
426 // Loop checking, 13.08.2015, G.Cosmo 436 // Loop checking, 13.08.2015, G.Cosmo
427 437
428 G4double z0 = ( fZSections[iz+1].fZ + fZSect 438 G4double z0 = ( fZSections[iz+1].fZ + fZSections[iz].fZ )/2.0;
429 G4TwoVector p2(point.x(), point.y()); 439 G4TwoVector p2(point.x(), point.y());
430 G4double pscale = fKScales[iz]*(point.z()-z 440 G4double pscale = fKScales[iz]*(point.z()-z0) + fScale0s[iz];
431 G4TwoVector poffset = fKOffsets[iz]*(point.z 441 G4TwoVector poffset = fKOffsets[iz]*(point.z()-z0) + fOffset0s[iz];
432 442
433 // G4cout << point << " projected to " 443 // G4cout << point << " projected to "
434 // << iz << "-th z-segment polygon as 444 // << iz << "-th z-segment polygon as "
435 // << (p2 - poffset)/pscale << G4endl 445 // << (p2 - poffset)/pscale << G4endl;
436 446
437 // pscale is always >0 as it is an interpola 447 // pscale is always >0 as it is an interpolation between two
438 // positive scale values 448 // positive scale values
439 // 449 //
440 return (p2 - poffset)/pscale; 450 return (p2 - poffset)/pscale;
441 } 451 }
442 452
443 //____________________________________________ 453 //_____________________________________________________________________________
444 454
445 G4bool G4ExtrudedSolid::IsSameLine(const G4Two 455 G4bool G4ExtrudedSolid::IsSameLine(const G4TwoVector& p,
446 const G4Two 456 const G4TwoVector& l1,
447 const G4Two 457 const G4TwoVector& l2) const
448 { 458 {
449 // Return true if p is on the line through l 459 // Return true if p is on the line through l1, l2
450 460
451 if ( l1.x() == l2.x() ) 461 if ( l1.x() == l2.x() )
452 { 462 {
453 return std::fabs(p.x() - l1.x()) < kCarTol 463 return std::fabs(p.x() - l1.x()) < kCarToleranceHalf;
454 } 464 }
455 G4double slope= ((l2.y() - l1.y())/(l2.x() 465 G4double slope= ((l2.y() - l1.y())/(l2.x() - l1.x()));
456 G4double predy= l1.y() + slope *(p.x() - l 466 G4double predy= l1.y() + slope *(p.x() - l1.x());
457 G4double dy= p.y() - predy; 467 G4double dy= p.y() - predy;
458 468
459 // Calculate perpendicular distance 469 // Calculate perpendicular distance
460 // 470 //
461 // G4double perpD= std::fabs(dy) / std::sqr 471 // G4double perpD= std::fabs(dy) / std::sqrt( 1 + slope * slope );
462 // G4bool simpleComp= (perpD<kCarToleranc 472 // G4bool simpleComp= (perpD<kCarToleranceHalf);
463 473
464 // Check perpendicular distance vs toleranc 474 // Check perpendicular distance vs tolerance 'directly'
465 // 475 //
466 G4bool squareComp = (dy*dy < (1+slope*slope 476 G4bool squareComp = (dy*dy < (1+slope*slope)
467 * kCarToleranceHalf * kCa 477 * kCarToleranceHalf * kCarToleranceHalf);
468 478
469 // return simpleComp; 479 // return simpleComp;
470 return squareComp; 480 return squareComp;
471 } 481 }
472 482
473 //____________________________________________ 483 //_____________________________________________________________________________
474 484
475 G4bool G4ExtrudedSolid::IsSameLineSegment(cons 485 G4bool G4ExtrudedSolid::IsSameLineSegment(const G4TwoVector& p,
476 cons 486 const G4TwoVector& l1,
477 cons 487 const G4TwoVector& l2) const
478 { 488 {
479 // Return true if p is on the line through l 489 // Return true if p is on the line through l1, l2 and lies between
480 // l1 and l2 490 // l1 and l2
481 491
482 if ( p.x() < std::min(l1.x(), l2.x()) - kCar 492 if ( p.x() < std::min(l1.x(), l2.x()) - kCarToleranceHalf ||
483 p.x() > std::max(l1.x(), l2.x()) + kCar 493 p.x() > std::max(l1.x(), l2.x()) + kCarToleranceHalf ||
484 p.y() < std::min(l1.y(), l2.y()) - kCar 494 p.y() < std::min(l1.y(), l2.y()) - kCarToleranceHalf ||
485 p.y() > std::max(l1.y(), l2.y()) + kCar 495 p.y() > std::max(l1.y(), l2.y()) + kCarToleranceHalf )
486 { 496 {
487 return false; 497 return false;
488 } 498 }
489 499
490 return IsSameLine(p, l1, l2); 500 return IsSameLine(p, l1, l2);
491 } 501 }
492 502
493 //____________________________________________ 503 //_____________________________________________________________________________
494 504
495 G4bool G4ExtrudedSolid::IsSameSide(const G4Two 505 G4bool G4ExtrudedSolid::IsSameSide(const G4TwoVector& p1,
496 const G4Two 506 const G4TwoVector& p2,
497 const G4Two 507 const G4TwoVector& l1,
498 const G4Two 508 const G4TwoVector& l2) const
499 { 509 {
500 // Return true if p1 and p2 are on the same 510 // Return true if p1 and p2 are on the same side of the line through l1, l2
501 511
502 return ( (p1.x() - l1.x()) * (l2.y() - l1. 512 return ( (p1.x() - l1.x()) * (l2.y() - l1.y())
503 - (l2.x() - l1.x()) * (p1.y() - l1.y( 513 - (l2.x() - l1.x()) * (p1.y() - l1.y()) )
504 * ( (p2.x() - l1.x()) * (l2.y() - l1. 514 * ( (p2.x() - l1.x()) * (l2.y() - l1.y())
505 - (l2.x() - l1.x()) * (p2.y() - l1.y( 515 - (l2.x() - l1.x()) * (p2.y() - l1.y()) ) > 0;
506 } 516 }
507 517
508 //____________________________________________ 518 //_____________________________________________________________________________
509 519
510 G4bool G4ExtrudedSolid::IsPointInside(const G4 520 G4bool G4ExtrudedSolid::IsPointInside(const G4TwoVector& a,
511 const G4 521 const G4TwoVector& b,
512 const G4 522 const G4TwoVector& c,
513 const G4 523 const G4TwoVector& p) const
514 { 524 {
515 // Return true if p is inside of triangle ab 525 // Return true if p is inside of triangle abc or on its edges,
516 // else returns false 526 // else returns false
517 527
518 // Check extent first 528 // Check extent first
519 // 529 //
520 if ( ( p.x() < a.x() && p.x() < b.x() && p.x 530 if ( ( p.x() < a.x() && p.x() < b.x() && p.x() < c.x() ) ||
521 ( p.x() > a.x() && p.x() > b.x() && p.x 531 ( p.x() > a.x() && p.x() > b.x() && p.x() > c.x() ) ||
522 ( p.y() < a.y() && p.y() < b.y() && p.y 532 ( p.y() < a.y() && p.y() < b.y() && p.y() < c.y() ) ||
523 ( p.y() > a.y() && p.y() > b.y() && p.y 533 ( p.y() > a.y() && p.y() > b.y() && p.y() > c.y() ) ) return false;
524 534
525 G4bool inside 535 G4bool inside
526 = IsSameSide(p, a, b, c) 536 = IsSameSide(p, a, b, c)
527 && IsSameSide(p, b, a, c) 537 && IsSameSide(p, b, a, c)
528 && IsSameSide(p, c, a, b); 538 && IsSameSide(p, c, a, b);
529 539
530 G4bool onEdge 540 G4bool onEdge
531 = IsSameLineSegment(p, a, b) 541 = IsSameLineSegment(p, a, b)
532 || IsSameLineSegment(p, b, c) 542 || IsSameLineSegment(p, b, c)
533 || IsSameLineSegment(p, c, a); 543 || IsSameLineSegment(p, c, a);
534 544
535 return inside || onEdge; 545 return inside || onEdge;
536 } 546 }
537 547
538 //____________________________________________ 548 //_____________________________________________________________________________
539 549
540 G4double 550 G4double
541 G4ExtrudedSolid::GetAngle(const G4TwoVector& p 551 G4ExtrudedSolid::GetAngle(const G4TwoVector& po,
542 const G4TwoVector& p 552 const G4TwoVector& pa,
543 const G4TwoVector& p 553 const G4TwoVector& pb) const
544 { 554 {
545 // Return the angle of the vertex in po 555 // Return the angle of the vertex in po
546 556
547 G4TwoVector t1 = pa - po; 557 G4TwoVector t1 = pa - po;
548 G4TwoVector t2 = pb - po; 558 G4TwoVector t2 = pb - po;
549 559
550 G4double result = (std::atan2(t1.y(), t1.x() 560 G4double result = (std::atan2(t1.y(), t1.x()) - std::atan2(t2.y(), t2.x()));
551 561
552 if ( result < 0 ) result += 2*pi; 562 if ( result < 0 ) result += 2*pi;
553 563
554 return result; 564 return result;
555 } 565 }
556 566
557 //____________________________________________ 567 //_____________________________________________________________________________
558 568
559 G4VFacet* 569 G4VFacet*
560 G4ExtrudedSolid::MakeDownFacet(G4int ind1, G4i 570 G4ExtrudedSolid::MakeDownFacet(G4int ind1, G4int ind2, G4int ind3) const
561 { 571 {
562 // Create a triangular facet from the polygo 572 // Create a triangular facet from the polygon points given by indices
563 // forming the down side ( the normal goes i 573 // forming the down side ( the normal goes in -z)
564 574
565 std::vector<G4ThreeVector> vertices; 575 std::vector<G4ThreeVector> vertices;
566 vertices.push_back(GetVertex(0, ind1)); 576 vertices.push_back(GetVertex(0, ind1));
567 vertices.push_back(GetVertex(0, ind2)); 577 vertices.push_back(GetVertex(0, ind2));
568 vertices.push_back(GetVertex(0, ind3)); 578 vertices.push_back(GetVertex(0, ind3));
569 579
570 // first vertex most left 580 // first vertex most left
571 // 581 //
572 G4ThreeVector cross 582 G4ThreeVector cross
573 = (vertices[1]-vertices[0]).cross(vertices 583 = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
574 584
575 if ( cross.z() > 0.0 ) 585 if ( cross.z() > 0.0 )
576 { 586 {
577 // vertices ordered clock wise has to be r 587 // vertices ordered clock wise has to be reordered
578 588
579 // G4cout << "G4ExtrudedSolid::MakeDownFac 589 // G4cout << "G4ExtrudedSolid::MakeDownFacet: reordering vertices "
580 // << ind1 << ", " << ind2 << ", " 590 // << ind1 << ", " << ind2 << ", " << ind3 << G4endl;
581 591
582 G4ThreeVector tmp = vertices[1]; 592 G4ThreeVector tmp = vertices[1];
583 vertices[1] = vertices[2]; 593 vertices[1] = vertices[2];
584 vertices[2] = tmp; 594 vertices[2] = tmp;
585 } 595 }
586 596
587 return new G4TriangularFacet(vertices[0], ve 597 return new G4TriangularFacet(vertices[0], vertices[1],
588 vertices[2], AB 598 vertices[2], ABSOLUTE);
589 } 599 }
590 600
591 //____________________________________________ 601 //_____________________________________________________________________________
592 602
593 G4VFacet* 603 G4VFacet*
594 G4ExtrudedSolid::MakeUpFacet(G4int ind1, G4int 604 G4ExtrudedSolid::MakeUpFacet(G4int ind1, G4int ind2, G4int ind3) const
595 { 605 {
596 // Creates a triangular facet from the polyg 606 // Creates a triangular facet from the polygon points given by indices
597 // forming the upper side ( z>0 ) 607 // forming the upper side ( z>0 )
598 608
599 std::vector<G4ThreeVector> vertices; 609 std::vector<G4ThreeVector> vertices;
600 vertices.push_back(GetVertex((G4int)fNz-1, i 610 vertices.push_back(GetVertex((G4int)fNz-1, ind1));
601 vertices.push_back(GetVertex((G4int)fNz-1, i 611 vertices.push_back(GetVertex((G4int)fNz-1, ind2));
602 vertices.push_back(GetVertex((G4int)fNz-1, i 612 vertices.push_back(GetVertex((G4int)fNz-1, ind3));
603 613
604 // first vertex most left 614 // first vertex most left
605 // 615 //
606 G4ThreeVector cross 616 G4ThreeVector cross
607 = (vertices[1]-vertices[0]).cross(vertices 617 = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
608 618
609 if ( cross.z() < 0.0 ) 619 if ( cross.z() < 0.0 )
610 { 620 {
611 // vertices ordered clock wise has to be r 621 // vertices ordered clock wise has to be reordered
612 622
613 // G4cout << "G4ExtrudedSolid::MakeUpFacet 623 // G4cout << "G4ExtrudedSolid::MakeUpFacet: reordering vertices "
614 // << ind1 << ", " << ind2 << ", " 624 // << ind1 << ", " << ind2 << ", " << ind3 << G4endl;
615 625
616 G4ThreeVector tmp = vertices[1]; 626 G4ThreeVector tmp = vertices[1];
617 vertices[1] = vertices[2]; 627 vertices[1] = vertices[2];
618 vertices[2] = tmp; 628 vertices[2] = tmp;
619 } 629 }
620 630
621 return new G4TriangularFacet(vertices[0], ve 631 return new G4TriangularFacet(vertices[0], vertices[1],
622 vertices[2], AB 632 vertices[2], ABSOLUTE);
623 } 633 }
624 634
625 //____________________________________________ 635 //_____________________________________________________________________________
626 636
627 G4bool G4ExtrudedSolid::AddGeneralPolygonFacet 637 G4bool G4ExtrudedSolid::AddGeneralPolygonFacets()
628 { 638 {
629 // Decompose polygonal sides in triangular f 639 // Decompose polygonal sides in triangular facets
630 640
631 typedef std::pair < G4TwoVector, G4int > Ver 641 typedef std::pair < G4TwoVector, G4int > Vertex;
632 642
633 static const G4double kAngTolerance = 643 static const G4double kAngTolerance =
634 G4GeometryTolerance::GetInstance()->GetAng 644 G4GeometryTolerance::GetInstance()->GetAngularTolerance();
635 645
636 // Fill one more vector 646 // Fill one more vector
637 // 647 //
638 std::vector< Vertex > verticesToBeDone; 648 std::vector< Vertex > verticesToBeDone;
639 for ( G4int i=0; i<(G4int)fNv; ++i ) 649 for ( G4int i=0; i<(G4int)fNv; ++i )
640 { 650 {
641 verticesToBeDone.emplace_back(fPolygon[i], << 651 verticesToBeDone.push_back(Vertex(fPolygon[i], i));
642 } 652 }
643 std::vector< Vertex > ears; 653 std::vector< Vertex > ears;
644 654
645 auto c1 = verticesToBeDone.begin(); << 655 std::vector< Vertex >::iterator c1 = verticesToBeDone.begin();
646 auto c2 = c1+1; << 656 std::vector< Vertex >::iterator c2 = c1+1;
647 auto c3 = c1+2; << 657 std::vector< Vertex >::iterator c3 = c1+2;
648 while ( verticesToBeDone.size()>2 ) // Lo 658 while ( verticesToBeDone.size()>2 ) // Loop checking, 13.08.2015, G.Cosmo
649 { 659 {
650 660
651 // G4cout << "Looking at triangle : " 661 // G4cout << "Looking at triangle : "
652 // << c1->second << " " << c2->se 662 // << c1->second << " " << c2->second
653 // << " " << c3->second << G4endl; 663 // << " " << c3->second << G4endl;
654 //G4cout << "Looking at triangle : " 664 //G4cout << "Looking at triangle : "
655 // << c1->first << " " << c2->firs 665 // << c1->first << " " << c2->first
656 // << " " << c3->first << G4endl; 666 // << " " << c3->first << G4endl;
657 667
658 // skip concave vertices 668 // skip concave vertices
659 // 669 //
660 G4double angle = GetAngle(c2->first, c3->f 670 G4double angle = GetAngle(c2->first, c3->first, c1->first);
661 671
662 //G4cout << "angle " << angle << G4endl; 672 //G4cout << "angle " << angle << G4endl;
663 673
664 std::size_t counter = 0; 674 std::size_t counter = 0;
665 while ( angle >= (pi-kAngTolerance) ) // 675 while ( angle >= (pi-kAngTolerance) ) // Loop checking, 13.08.2015, G.Cosmo
666 { 676 {
667 // G4cout << "Skipping concave vertex " 677 // G4cout << "Skipping concave vertex " << c2->second << G4endl;
668 678
669 // try next three consecutive vertices 679 // try next three consecutive vertices
670 // 680 //
671 c1 = c2; 681 c1 = c2;
672 c2 = c3; 682 c2 = c3;
673 ++c3; 683 ++c3;
674 if ( c3 == verticesToBeDone.end() ) { c3 684 if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
675 685
676 //G4cout << "Looking at triangle : " 686 //G4cout << "Looking at triangle : "
677 // << c1->first << " " << c2->firs 687 // << c1->first << " " << c2->first
678 // << " " << c3->first << G4endl 688 // << " " << c3->first << G4endl;
679 689
680 angle = GetAngle(c2->first, c3->first, c 690 angle = GetAngle(c2->first, c3->first, c1->first);
681 //G4cout << "angle " << angle << G4endl 691 //G4cout << "angle " << angle << G4endl;
682 692
683 ++counter; 693 ++counter;
684 694
685 if ( counter > fNv ) 695 if ( counter > fNv )
686 { 696 {
687 G4Exception("G4ExtrudedSolid::AddGener 697 G4Exception("G4ExtrudedSolid::AddGeneralPolygonFacets",
688 "GeomSolids0003", FatalExc 698 "GeomSolids0003", FatalException,
689 "Triangularisation has fai 699 "Triangularisation has failed.");
690 break; 700 break;
691 } 701 }
692 } 702 }
693 703
694 G4bool good = true; 704 G4bool good = true;
695 for ( auto it=verticesToBeDone.cbegin(); i 705 for ( auto it=verticesToBeDone.cbegin(); it!=verticesToBeDone.cend(); ++it )
696 { 706 {
697 // skip vertices of tested triangle 707 // skip vertices of tested triangle
698 // 708 //
699 if ( it == c1 || it == c2 || it == c3 ) 709 if ( it == c1 || it == c2 || it == c3 ) { continue; }
700 710
701 if ( IsPointInside(c1->first, c2->first, 711 if ( IsPointInside(c1->first, c2->first, c3->first, it->first) )
702 { 712 {
703 // G4cout << "Point " << it->second << 713 // G4cout << "Point " << it->second << " is inside" << G4endl;
704 good = false; 714 good = false;
705 715
706 // try next three consecutive vertices 716 // try next three consecutive vertices
707 // 717 //
708 c1 = c2; 718 c1 = c2;
709 c2 = c3; 719 c2 = c3;
710 ++c3; 720 ++c3;
711 if ( c3 == verticesToBeDone.end() ) { 721 if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
712 break; 722 break;
713 } 723 }
714 // else 724 // else
715 // { G4cout << "Point " << it->second 725 // { G4cout << "Point " << it->second << " is outside" << G4endl; }
716 } 726 }
717 if ( good ) 727 if ( good )
718 { 728 {
719 // all points are outside triangle, we c 729 // all points are outside triangle, we can make a facet
720 730
721 // G4cout << "Found triangle : " 731 // G4cout << "Found triangle : "
722 // << c1->second << " " << c2->s 732 // << c1->second << " " << c2->second
723 // << " " << c3->second << G4end 733 // << " " << c3->second << G4endl;
724 734
725 G4bool result; 735 G4bool result;
726 result = AddFacet( MakeDownFacet(c1->sec 736 result = AddFacet( MakeDownFacet(c1->second, c2->second, c3->second) );
727 if ( ! result ) { return false; } 737 if ( ! result ) { return false; }
728 738
729 result = AddFacet( MakeUpFacet(c1->secon 739 result = AddFacet( MakeUpFacet(c1->second, c2->second, c3->second) );
730 if ( ! result ) { return false; } 740 if ( ! result ) { return false; }
731 741
732 std::vector<G4int> triangle(3); 742 std::vector<G4int> triangle(3);
733 triangle[0] = c1->second; 743 triangle[0] = c1->second;
734 triangle[1] = c2->second; 744 triangle[1] = c2->second;
735 triangle[2] = c3->second; 745 triangle[2] = c3->second;
736 fTriangles.push_back(std::move(triangle) << 746 fTriangles.push_back(triangle);
737 747
738 // remove the ear point from verticesToB 748 // remove the ear point from verticesToBeDone
739 // 749 //
740 verticesToBeDone.erase(c2); 750 verticesToBeDone.erase(c2);
741 c1 = verticesToBeDone.begin(); 751 c1 = verticesToBeDone.begin();
742 c2 = c1+1; 752 c2 = c1+1;
743 c3 = c1+2; 753 c3 = c1+2;
744 } 754 }
745 } 755 }
746 return true; 756 return true;
747 } 757 }
748 758
749 //____________________________________________ 759 //_____________________________________________________________________________
750 760
751 G4bool G4ExtrudedSolid::MakeFacets() 761 G4bool G4ExtrudedSolid::MakeFacets()
752 { 762 {
753 // Define facets 763 // Define facets
754 764
755 G4bool good; 765 G4bool good;
756 766
757 // Decomposition of polygonal sides in the f 767 // Decomposition of polygonal sides in the facets
758 // 768 //
759 if ( fNv == 3 ) 769 if ( fNv == 3 )
760 { 770 {
761 good = AddFacet( new G4TriangularFacet( Ge 771 good = AddFacet( new G4TriangularFacet( GetVertex(0, 0), GetVertex(0, 1),
762 Ge 772 GetVertex(0, 2), ABSOLUTE) );
763 if ( ! good ) { return false; } 773 if ( ! good ) { return false; }
764 774
765 good = AddFacet( new G4TriangularFacet( Ge 775 good = AddFacet( new G4TriangularFacet( GetVertex((G4int)fNz-1, 2),
766 Ge 776 GetVertex((G4int)fNz-1, 1),
767 Ge 777 GetVertex((G4int)fNz-1, 0),
768 AB 778 ABSOLUTE) );
769 if ( ! good ) { return false; } 779 if ( ! good ) { return false; }
770 780
771 std::vector<G4int> triangle(3); 781 std::vector<G4int> triangle(3);
772 triangle[0] = 0; 782 triangle[0] = 0;
773 triangle[1] = 1; 783 triangle[1] = 1;
774 triangle[2] = 2; 784 triangle[2] = 2;
775 fTriangles.push_back(std::move(triangle)); << 785 fTriangles.push_back(triangle);
776 } 786 }
777 787
778 else if ( fNv == 4 ) 788 else if ( fNv == 4 )
779 { 789 {
780 good = AddFacet( new G4QuadrangularFacet( 790 good = AddFacet( new G4QuadrangularFacet( GetVertex(0, 0),GetVertex(0, 1),
781 791 GetVertex(0, 2),GetVertex(0, 3),
782 792 ABSOLUTE) );
783 if ( ! good ) { return false; } 793 if ( ! good ) { return false; }
784 794
785 good = AddFacet( new G4QuadrangularFacet( 795 good = AddFacet( new G4QuadrangularFacet( GetVertex((G4int)fNz-1, 3),
786 796 GetVertex((G4int)fNz-1, 2),
787 797 GetVertex((G4int)fNz-1, 1),
788 798 GetVertex((G4int)fNz-1, 0),
789 799 ABSOLUTE) );
790 if ( ! good ) { return false; } 800 if ( ! good ) { return false; }
791 801
792 std::vector<G4int> triangle1(3); 802 std::vector<G4int> triangle1(3);
793 triangle1[0] = 0; 803 triangle1[0] = 0;
794 triangle1[1] = 1; 804 triangle1[1] = 1;
795 triangle1[2] = 2; 805 triangle1[2] = 2;
796 fTriangles.push_back(std::move(triangle1)) << 806 fTriangles.push_back(triangle1);
797 807
798 std::vector<G4int> triangle2(3); 808 std::vector<G4int> triangle2(3);
799 triangle2[0] = 0; 809 triangle2[0] = 0;
800 triangle2[1] = 2; 810 triangle2[1] = 2;
801 triangle2[2] = 3; 811 triangle2[2] = 3;
802 fTriangles.push_back(std::move(triangle2)) << 812 fTriangles.push_back(triangle2);
803 } 813 }
804 else 814 else
805 { 815 {
806 good = AddGeneralPolygonFacets(); 816 good = AddGeneralPolygonFacets();
807 if ( ! good ) { return false; } 817 if ( ! good ) { return false; }
808 } 818 }
809 819
810 // The quadrangular sides 820 // The quadrangular sides
811 // 821 //
812 for ( G4int iz = 0; iz < (G4int)fNz-1; ++iz 822 for ( G4int iz = 0; iz < (G4int)fNz-1; ++iz )
813 { 823 {
814 for ( G4int i = 0; i < (G4int)fNv; ++i ) 824 for ( G4int i = 0; i < (G4int)fNv; ++i )
815 { 825 {
816 G4int j = (i+1) % fNv; 826 G4int j = (i+1) % fNv;
817 good = AddFacet( new G4QuadrangularFacet 827 good = AddFacet( new G4QuadrangularFacet
818 ( GetVertex(iz, j), Ge 828 ( GetVertex(iz, j), GetVertex(iz, i),
819 GetVertex(iz+1, i), 829 GetVertex(iz+1, i), GetVertex(iz+1, j), ABSOLUTE) );
820 if ( ! good ) { return false; } 830 if ( ! good ) { return false; }
821 } 831 }
822 } 832 }
823 833
824 SetSolidClosed(true); 834 SetSolidClosed(true);
825 835
826 return good; 836 return good;
827 } 837 }
828 838
829 //____________________________________________ 839 //_____________________________________________________________________________
830 840
831 G4GeometryType G4ExtrudedSolid::GetEntityType 841 G4GeometryType G4ExtrudedSolid::GetEntityType () const
832 { 842 {
833 // Return entity type 843 // Return entity type
834 844
835 return fGeometryType; 845 return fGeometryType;
836 } 846 }
837 847
838 //____________________________________________ 848 //_____________________________________________________________________________
839 849
840 G4bool G4ExtrudedSolid::IsFaceted () const <<
841 { <<
842 return true; <<
843 } <<
844 <<
845 //____________________________________________ <<
846 <<
847 G4VSolid* G4ExtrudedSolid::Clone() const 850 G4VSolid* G4ExtrudedSolid::Clone() const
848 { 851 {
849 return new G4ExtrudedSolid(*this); 852 return new G4ExtrudedSolid(*this);
850 } 853 }
851 854
852 //____________________________________________ 855 //_____________________________________________________________________________
853 856
854 EInside G4ExtrudedSolid::Inside(const G4ThreeV 857 EInside G4ExtrudedSolid::Inside(const G4ThreeVector &p) const
855 { 858 {
856 switch (fSolidType) 859 switch (fSolidType)
857 { 860 {
858 case 1: // convex right prism 861 case 1: // convex right prism
859 { 862 {
860 G4double dist = std::max(fZSections[0].f 863 G4double dist = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
861 if (dist > kCarToleranceHalf) { return 864 if (dist > kCarToleranceHalf) { return kOutside; }
862 865
863 std::size_t np = fPlanes.size(); 866 std::size_t np = fPlanes.size();
864 for (std::size_t i=0; i<np; ++i) 867 for (std::size_t i=0; i<np; ++i)
865 { 868 {
866 G4double dd = fPlanes[i].a*p.x() + fPl 869 G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
867 if (dd > dist) { dist = dd; } 870 if (dd > dist) { dist = dd; }
868 } 871 }
869 if (dist > kCarToleranceHalf) { return 872 if (dist > kCarToleranceHalf) { return kOutside; }
870 return (dist > -kCarToleranceHalf) ? kSu 873 return (dist > -kCarToleranceHalf) ? kSurface : kInside;
871 } 874 }
872 case 2: // non-convex right prism 875 case 2: // non-convex right prism
873 { 876 {
874 G4double distz = std::max(fZSections[0]. 877 G4double distz = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
875 if (distz > kCarToleranceHalf) { return 878 if (distz > kCarToleranceHalf) { return kOutside; }
876 879
877 G4bool in = PointInPolygon(p); 880 G4bool in = PointInPolygon(p);
878 if (distz > -kCarToleranceHalf && in) { 881 if (distz > -kCarToleranceHalf && in) { return kSurface; }
879 882
880 G4double dd = DistanceToPolygonSqr(p) - 883 G4double dd = DistanceToPolygonSqr(p) - kCarToleranceHalf*kCarToleranceHalf;
881 if (in) 884 if (in)
882 { 885 {
883 return (dd >= 0) ? kInside : kSurface; 886 return (dd >= 0) ? kInside : kSurface;
884 } 887 }
885 else 888 else
886 { 889 {
887 return (dd > 0) ? kOutside : kSurface; 890 return (dd > 0) ? kOutside : kSurface;
888 } 891 }
889 } 892 }
890 } 893 }
891 894
892 // Override the base class function as it f 895 // Override the base class function as it fails in case of concave polygon.
893 // Project the point in the original polygon 896 // Project the point in the original polygon scale and check if it is inside
894 // for each triangle. 897 // for each triangle.
895 898
896 // Check first if outside extent 899 // Check first if outside extent
897 // 900 //
898 if ( p.x() < GetMinXExtent() - kCarTolerance 901 if ( p.x() < GetMinXExtent() - kCarToleranceHalf ||
899 p.x() > GetMaxXExtent() + kCarTolerance 902 p.x() > GetMaxXExtent() + kCarToleranceHalf ||
900 p.y() < GetMinYExtent() - kCarTolerance 903 p.y() < GetMinYExtent() - kCarToleranceHalf ||
901 p.y() > GetMaxYExtent() + kCarTolerance 904 p.y() > GetMaxYExtent() + kCarToleranceHalf ||
902 p.z() < GetMinZExtent() - kCarTolerance 905 p.z() < GetMinZExtent() - kCarToleranceHalf ||
903 p.z() > GetMaxZExtent() + kCarTolerance 906 p.z() > GetMaxZExtent() + kCarToleranceHalf )
904 { 907 {
905 // G4cout << "G4ExtrudedSolid::Outside ext 908 // G4cout << "G4ExtrudedSolid::Outside extent: " << p << G4endl;
906 return kOutside; 909 return kOutside;
907 } 910 }
908 911
909 // Project point p(z) to the polygon scale p 912 // Project point p(z) to the polygon scale p0
910 // 913 //
911 G4TwoVector pscaled = ProjectPoint(p); 914 G4TwoVector pscaled = ProjectPoint(p);
912 915
913 // Check if on surface of polygon 916 // Check if on surface of polygon
914 // 917 //
915 for ( G4int i=0; i<(G4int)fNv; ++i ) 918 for ( G4int i=0; i<(G4int)fNv; ++i )
916 { 919 {
917 G4int j = (i+1) % fNv; 920 G4int j = (i+1) % fNv;
918 if ( IsSameLineSegment(pscaled, fPolygon[i 921 if ( IsSameLineSegment(pscaled, fPolygon[i], fPolygon[j]) )
919 { 922 {
920 // G4cout << "G4ExtrudedSolid::Inside re 923 // G4cout << "G4ExtrudedSolid::Inside return Surface (on polygon) "
921 // << G4endl; 924 // << G4endl;
922 925
923 return kSurface; 926 return kSurface;
924 } 927 }
925 } 928 }
926 929
927 // Now check if inside triangles 930 // Now check if inside triangles
928 // 931 //
929 auto it = fTriangles.cbegin(); 932 auto it = fTriangles.cbegin();
930 G4bool inside = false; 933 G4bool inside = false;
931 do // Loop checking, 13.08.2015, G.Cosmo 934 do // Loop checking, 13.08.2015, G.Cosmo
932 { 935 {
933 if ( IsPointInside(fPolygon[(*it)[0]], fPo 936 if ( IsPointInside(fPolygon[(*it)[0]], fPolygon[(*it)[1]],
934 fPolygon[(*it)[2]], psc 937 fPolygon[(*it)[2]], pscaled) ) { inside = true; }
935 ++it; 938 ++it;
936 } while ( (!inside) && (it != fTriangles.cen << 939 } while ( (inside == false) && (it != fTriangles.cend()) );
937 940
938 if ( inside ) 941 if ( inside )
939 { 942 {
940 // Check if on surface of z sides 943 // Check if on surface of z sides
941 // 944 //
942 if ( std::fabs( p.z() - fZSections[0].fZ ) 945 if ( std::fabs( p.z() - fZSections[0].fZ ) < kCarToleranceHalf ||
943 std::fabs( p.z() - fZSections[fNz-1]. 946 std::fabs( p.z() - fZSections[fNz-1].fZ ) < kCarToleranceHalf )
944 { 947 {
945 // G4cout << "G4ExtrudedSolid::Inside re 948 // G4cout << "G4ExtrudedSolid::Inside return Surface (on z side)"
946 // << G4endl; 949 // << G4endl;
947 950
948 return kSurface; 951 return kSurface;
949 } 952 }
950 953
951 // G4cout << "G4ExtrudedSolid::Inside retu 954 // G4cout << "G4ExtrudedSolid::Inside return Inside" << G4endl;
952 955
953 return kInside; 956 return kInside;
954 } 957 }
955 958
956 // G4cout << "G4ExtrudedSolid::Inside return 959 // G4cout << "G4ExtrudedSolid::Inside return Outside " << G4endl;
957 960
958 return kOutside; 961 return kOutside;
959 } 962 }
960 963
961 //____________________________________________ 964 //_____________________________________________________________________________
962 965
963 G4ThreeVector G4ExtrudedSolid::SurfaceNormal(c 966 G4ThreeVector G4ExtrudedSolid::SurfaceNormal(const G4ThreeVector& p) const
964 { 967 {
965 G4int nsurf = 0; 968 G4int nsurf = 0;
966 G4double nx = 0., ny = 0., nz = 0.; 969 G4double nx = 0., ny = 0., nz = 0.;
967 switch (fSolidType) 970 switch (fSolidType)
968 { 971 {
969 case 1: // convex right prism 972 case 1: // convex right prism
970 { 973 {
971 if (std::abs(p.z() - fZSections[0].fZ) < 974 if (std::abs(p.z() - fZSections[0].fZ) <= kCarToleranceHalf)
972 { 975 {
973 nz = -1; ++nsurf; 976 nz = -1; ++nsurf;
974 } 977 }
975 if (std::abs(p.z() - fZSections[1].fZ) < 978 if (std::abs(p.z() - fZSections[1].fZ) <= kCarToleranceHalf)
976 { 979 {
977 nz = 1; ++nsurf; 980 nz = 1; ++nsurf;
978 } 981 }
979 for (std::size_t i=0; i<fNv; ++i) 982 for (std::size_t i=0; i<fNv; ++i)
980 { 983 {
981 G4double dd = fPlanes[i].a*p.x() + fPl 984 G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
982 if (std::abs(dd) > kCarToleranceHalf) 985 if (std::abs(dd) > kCarToleranceHalf) continue;
983 nx += fPlanes[i].a; 986 nx += fPlanes[i].a;
984 ny += fPlanes[i].b; 987 ny += fPlanes[i].b;
985 ++nsurf; 988 ++nsurf;
986 } 989 }
987 break; 990 break;
988 } 991 }
989 case 2: // non-convex right prism 992 case 2: // non-convex right prism
990 { 993 {
991 if (std::abs(p.z() - fZSections[0].fZ) < 994 if (std::abs(p.z() - fZSections[0].fZ) <= kCarToleranceHalf)
992 { 995 {
993 nz = -1; ++nsurf; 996 nz = -1; ++nsurf;
994 } 997 }
995 if (std::abs(p.z() - fZSections[1].fZ) < 998 if (std::abs(p.z() - fZSections[1].fZ) <= kCarToleranceHalf)
996 { 999 {
997 nz = 1; ++nsurf; 1000 nz = 1; ++nsurf;
998 } 1001 }
999 1002
1000 G4double sqrCarToleranceHalf = kCarTole 1003 G4double sqrCarToleranceHalf = kCarToleranceHalf*kCarToleranceHalf;
1001 for (std::size_t i=0, k=fNv-1; i<fNv; k 1004 for (std::size_t i=0, k=fNv-1; i<fNv; k=i++)
1002 { 1005 {
1003 G4double ix = p.x() - fPolygon[i].x() 1006 G4double ix = p.x() - fPolygon[i].x();
1004 G4double iy = p.y() - fPolygon[i].y() 1007 G4double iy = p.y() - fPolygon[i].y();
1005 G4double u = fPlanes[i].a*iy - fPlan 1008 G4double u = fPlanes[i].a*iy - fPlanes[i].b*ix;
1006 if (u < 0) 1009 if (u < 0)
1007 { 1010 {
1008 if (ix*ix + iy*iy > sqrCarTolerance 1011 if (ix*ix + iy*iy > sqrCarToleranceHalf) continue;
1009 } 1012 }
1010 else if (u > fLengths[i]) 1013 else if (u > fLengths[i])
1011 { 1014 {
1012 G4double kx = p.x() - fPolygon[k].x 1015 G4double kx = p.x() - fPolygon[k].x();
1013 G4double ky = p.y() - fPolygon[k].y 1016 G4double ky = p.y() - fPolygon[k].y();
1014 if (kx*kx + ky*ky > sqrCarTolerance 1017 if (kx*kx + ky*ky > sqrCarToleranceHalf) continue;
1015 } 1018 }
1016 else 1019 else
1017 { 1020 {
1018 G4double dd = fPlanes[i].a*p.x() + 1021 G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
1019 if (dd*dd > sqrCarToleranceHalf) co 1022 if (dd*dd > sqrCarToleranceHalf) continue;
1020 } 1023 }
1021 nx += fPlanes[i].a; 1024 nx += fPlanes[i].a;
1022 ny += fPlanes[i].b; 1025 ny += fPlanes[i].b;
1023 ++nsurf; 1026 ++nsurf;
1024 } 1027 }
1025 break; 1028 break;
1026 } 1029 }
1027 default: 1030 default:
1028 { 1031 {
1029 return G4TessellatedSolid::SurfaceNorma 1032 return G4TessellatedSolid::SurfaceNormal(p);
1030 } 1033 }
1031 } 1034 }
1032 1035
1033 // Return normal (right prism) 1036 // Return normal (right prism)
1034 // 1037 //
1035 if (nsurf == 1) 1038 if (nsurf == 1)
1036 { 1039 {
1037 return { nx,ny,nz }; << 1040 return G4ThreeVector(nx,ny,nz);
1038 } 1041 }
1039 else if (nsurf != 0) // edge or corner 1042 else if (nsurf != 0) // edge or corner
1040 { 1043 {
1041 return G4ThreeVector(nx,ny,nz).unit(); 1044 return G4ThreeVector(nx,ny,nz).unit();
1042 } 1045 }
1043 else 1046 else
1044 { 1047 {
1045 // Point is not on the surface, compute a 1048 // Point is not on the surface, compute approximate normal
1046 // 1049 //
1047 #ifdef G4CSGDEBUG 1050 #ifdef G4CSGDEBUG
1048 std::ostringstream message; 1051 std::ostringstream message;
1049 G4long oldprc = message.precision(16); 1052 G4long oldprc = message.precision(16);
1050 message << "Point p is not on surface (!? 1053 message << "Point p is not on surface (!?) of solid: "
1051 << GetName() << G4endl; 1054 << GetName() << G4endl;
1052 message << "Position:\n"; 1055 message << "Position:\n";
1053 message << " p.x() = " << p.x()/mm << " 1056 message << " p.x() = " << p.x()/mm << " mm\n";
1054 message << " p.y() = " << p.y()/mm << " 1057 message << " p.y() = " << p.y()/mm << " mm\n";
1055 message << " p.z() = " << p.z()/mm << " 1058 message << " p.z() = " << p.z()/mm << " mm";
1056 G4cout.precision(oldprc) ; 1059 G4cout.precision(oldprc) ;
1057 G4Exception("G4TesselatedSolid::SurfaceNo 1060 G4Exception("G4TesselatedSolid::SurfaceNormal(p)", "GeomSolids1002",
1058 JustWarning, message ); 1061 JustWarning, message );
1059 DumpInfo(); 1062 DumpInfo();
1060 #endif 1063 #endif
1061 return ApproxSurfaceNormal(p); 1064 return ApproxSurfaceNormal(p);
1062 } 1065 }
1063 } 1066 }
1064 1067
1065 //___________________________________________ 1068 //_____________________________________________________________________________
1066 1069
1067 G4ThreeVector G4ExtrudedSolid::ApproxSurfaceN 1070 G4ThreeVector G4ExtrudedSolid::ApproxSurfaceNormal(const G4ThreeVector& p) const
1068 { 1071 {
1069 // This method is valid only for right pris 1072 // This method is valid only for right prisms and
1070 // normally should not be called 1073 // normally should not be called
1071 1074
1072 if (fSolidType == 1 || fSolidType == 2) 1075 if (fSolidType == 1 || fSolidType == 2)
1073 { 1076 {
1074 // Find distances to z-planes 1077 // Find distances to z-planes
1075 // 1078 //
1076 G4double dz0 = fZSections[0].fZ - p.z(); 1079 G4double dz0 = fZSections[0].fZ - p.z();
1077 G4double dz1 = p.z() - fZSections[1].fZ; 1080 G4double dz1 = p.z() - fZSections[1].fZ;
1078 G4double ddz0 = dz0*dz0; 1081 G4double ddz0 = dz0*dz0;
1079 G4double ddz1 = dz1*dz1; 1082 G4double ddz1 = dz1*dz1;
1080 1083
1081 // Find nearest lateral side and distance 1084 // Find nearest lateral side and distance to it
1082 // 1085 //
1083 std::size_t iside = 0; 1086 std::size_t iside = 0;
1084 G4double dd = DBL_MAX; 1087 G4double dd = DBL_MAX;
1085 for (std::size_t i=0, k=fNv-1; i<fNv; k=i 1088 for (std::size_t i=0, k=fNv-1; i<fNv; k=i++)
1086 { 1089 {
1087 G4double ix = p.x() - fPolygon[i].x(); 1090 G4double ix = p.x() - fPolygon[i].x();
1088 G4double iy = p.y() - fPolygon[i].y(); 1091 G4double iy = p.y() - fPolygon[i].y();
1089 G4double u = fPlanes[i].a*iy - fPlanes 1092 G4double u = fPlanes[i].a*iy - fPlanes[i].b*ix;
1090 if (u < 0) 1093 if (u < 0)
1091 { 1094 {
1092 G4double tmp = ix*ix + iy*iy; 1095 G4double tmp = ix*ix + iy*iy;
1093 if (tmp < dd) { dd = tmp; iside = i; 1096 if (tmp < dd) { dd = tmp; iside = i; }
1094 } 1097 }
1095 else if (u > fLengths[i]) 1098 else if (u > fLengths[i])
1096 { 1099 {
1097 G4double kx = p.x() - fPolygon[k].x() 1100 G4double kx = p.x() - fPolygon[k].x();
1098 G4double ky = p.y() - fPolygon[k].y() 1101 G4double ky = p.y() - fPolygon[k].y();
1099 G4double tmp = kx*kx + ky*ky; 1102 G4double tmp = kx*kx + ky*ky;
1100 if (tmp < dd) { dd = tmp; iside = i; 1103 if (tmp < dd) { dd = tmp; iside = i; }
1101 } 1104 }
1102 else 1105 else
1103 { 1106 {
1104 G4double tmp = fPlanes[i].a*p.x() + f 1107 G4double tmp = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
1105 tmp *= tmp; 1108 tmp *= tmp;
1106 if (tmp < dd) { dd = tmp; iside = i; 1109 if (tmp < dd) { dd = tmp; iside = i; }
1107 } 1110 }
1108 } 1111 }
1109 1112
1110 // Find region 1113 // Find region
1111 // 1114 //
1112 // 3 | 1 | 3 1115 // 3 | 1 | 3
1113 // ----+-------+---- 1116 // ----+-------+----
1114 // 2 | 0 | 2 1117 // 2 | 0 | 2
1115 // ----+-------+---- 1118 // ----+-------+----
1116 // 3 | 1 | 3 1119 // 3 | 1 | 3
1117 // 1120 //
1118 G4int iregion = 0; 1121 G4int iregion = 0;
1119 if (std::max(dz0,dz1) > 0) iregion = 1; 1122 if (std::max(dz0,dz1) > 0) iregion = 1;
1120 1123
1121 G4bool in = PointInPolygon(p); 1124 G4bool in = PointInPolygon(p);
1122 if (!in) iregion += 2; 1125 if (!in) iregion += 2;
1123 1126
1124 // Return normal 1127 // Return normal
1125 // 1128 //
1126 switch (iregion) 1129 switch (iregion)
1127 { 1130 {
1128 case 0: 1131 case 0:
1129 { 1132 {
1130 if (ddz0 <= ddz1 && ddz0 <= dd) retur << 1133 if (ddz0 <= ddz1 && ddz0 <= dd) return G4ThreeVector(0, 0,-1);
1131 if (ddz1 <= ddz0 && ddz1 <= dd) retur << 1134 if (ddz1 <= ddz0 && ddz1 <= dd) return G4ThreeVector(0, 0, 1);
1132 return { fPlanes[iside].a, fPlanes[is << 1135 return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, 0);
1133 } 1136 }
1134 case 1: 1137 case 1:
1135 { 1138 {
1136 return { 0, 0, (G4double)((dz0 > dz1) << 1139 return G4ThreeVector(0, 0, (dz0 > dz1) ? -1 : 1);
1137 } 1140 }
1138 case 2: 1141 case 2:
1139 { 1142 {
1140 return { fPlanes[iside].a, fPlanes[is << 1143 return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, 0);
1141 } 1144 }
1142 case 3: 1145 case 3:
1143 { 1146 {
1144 G4double dzmax = std::max(dz0,dz1); 1147 G4double dzmax = std::max(dz0,dz1);
1145 if (dzmax*dzmax > dd) return { 0, 0, << 1148 if (dzmax*dzmax > dd) return G4ThreeVector(0,0,(dz0 > dz1) ? -1 : 1);
1146 return { fPlanes[iside].a, fPlanes[is << 1149 return G4ThreeVector(fPlanes[iside].a,fPlanes[iside].b, 0);
1147 } 1150 }
1148 } 1151 }
1149 } 1152 }
1150 return {0,0,0}; << 1153 return G4ThreeVector(0,0,0);
1151 } 1154 }
1152 1155
1153 //___________________________________________ 1156 //_____________________________________________________________________________
1154 1157
1155 G4double G4ExtrudedSolid::DistanceToIn(const 1158 G4double G4ExtrudedSolid::DistanceToIn(const G4ThreeVector& p,
1156 const 1159 const G4ThreeVector& v) const
1157 { 1160 {
1158 G4double z0 = fZSections[0].fZ; 1161 G4double z0 = fZSections[0].fZ;
1159 G4double z1 = fZSections[fNz-1].fZ; 1162 G4double z1 = fZSections[fNz-1].fZ;
1160 if ((p.z() <= z0 + kCarToleranceHalf) && v. 1163 if ((p.z() <= z0 + kCarToleranceHalf) && v.z() <= 0) return kInfinity;
1161 if ((p.z() >= z1 - kCarToleranceHalf) && v. 1164 if ((p.z() >= z1 - kCarToleranceHalf) && v.z() >= 0) return kInfinity;
1162 1165
1163 switch (fSolidType) 1166 switch (fSolidType)
1164 { 1167 {
1165 case 1: // convex right prism 1168 case 1: // convex right prism
1166 { 1169 {
1167 // Intersection with Z planes 1170 // Intersection with Z planes
1168 // 1171 //
1169 G4double dz = (z1 - z0)*0.5; 1172 G4double dz = (z1 - z0)*0.5;
1170 G4double pz = p.z() - dz - z0; 1173 G4double pz = p.z() - dz - z0;
1171 1174
1172 G4double invz = (v.z() == 0) ? DBL_MAX 1175 G4double invz = (v.z() == 0) ? DBL_MAX : -1./v.z();
1173 G4double ddz = (invz < 0) ? dz : -dz; 1176 G4double ddz = (invz < 0) ? dz : -dz;
1174 G4double tzmin = (pz + ddz)*invz; 1177 G4double tzmin = (pz + ddz)*invz;
1175 G4double tzmax = (pz - ddz)*invz; 1178 G4double tzmax = (pz - ddz)*invz;
1176 1179
1177 // Intersection with lateral planes 1180 // Intersection with lateral planes
1178 // 1181 //
1179 std::size_t np = fPlanes.size(); 1182 std::size_t np = fPlanes.size();
1180 G4double txmin = tzmin, txmax = tzmax; 1183 G4double txmin = tzmin, txmax = tzmax;
1181 for (std::size_t i=0; i<np; ++i) 1184 for (std::size_t i=0; i<np; ++i)
1182 { 1185 {
1183 G4double cosa = fPlanes[i].a*v.x()+fP 1186 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y();
1184 G4double dist = fPlanes[i].a*p.x()+fP 1187 G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].d;
1185 if (dist >= -kCarToleranceHalf) 1188 if (dist >= -kCarToleranceHalf)
1186 { 1189 {
1187 if (cosa >= 0) { return kInfinity; 1190 if (cosa >= 0) { return kInfinity; }
1188 G4double tmp = -dist/cosa; 1191 G4double tmp = -dist/cosa;
1189 if (txmin < tmp) { txmin = tmp; } 1192 if (txmin < tmp) { txmin = tmp; }
1190 } 1193 }
1191 else if (cosa > 0) 1194 else if (cosa > 0)
1192 { 1195 {
1193 G4double tmp = -dist/cosa; 1196 G4double tmp = -dist/cosa;
1194 if (txmax > tmp) { txmax = tmp; } 1197 if (txmax > tmp) { txmax = tmp; }
1195 } 1198 }
1196 } 1199 }
1197 1200
1198 // Find distance 1201 // Find distance
1199 // 1202 //
1200 G4double tmin = txmin, tmax = txmax; 1203 G4double tmin = txmin, tmax = txmax;
1201 if (tmax <= tmin + kCarToleranceHalf) 1204 if (tmax <= tmin + kCarToleranceHalf) // touch or no hit
1202 { 1205 {
1203 return kInfinity; 1206 return kInfinity;
1204 } 1207 }
1205 return (tmin < kCarToleranceHalf) ? 0. 1208 return (tmin < kCarToleranceHalf) ? 0. : tmin;
1206 } 1209 }
1207 case 2: // non-convex right prism 1210 case 2: // non-convex right prism
1208 { 1211 {
1209 } 1212 }
1210 } 1213 }
1211 return G4TessellatedSolid::DistanceToIn(p,v 1214 return G4TessellatedSolid::DistanceToIn(p,v);
1212 } 1215 }
1213 1216
1214 //___________________________________________ 1217 //_____________________________________________________________________________
1215 1218
1216 G4double G4ExtrudedSolid::DistanceToIn (const 1219 G4double G4ExtrudedSolid::DistanceToIn (const G4ThreeVector& p) const
1217 { 1220 {
1218 switch (fSolidType) 1221 switch (fSolidType)
1219 { 1222 {
1220 case 1: // convex right prism 1223 case 1: // convex right prism
1221 { 1224 {
1222 G4double dist = std::max(fZSections[0]. 1225 G4double dist = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
1223 std::size_t np = fPlanes.size(); 1226 std::size_t np = fPlanes.size();
1224 for (std::size_t i=0; i<np; ++i) 1227 for (std::size_t i=0; i<np; ++i)
1225 { 1228 {
1226 G4double dd = fPlanes[i].a*p.x() + fP 1229 G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
1227 if (dd > dist) dist = dd; 1230 if (dd > dist) dist = dd;
1228 } 1231 }
1229 return (dist > 0) ? dist : 0.; 1232 return (dist > 0) ? dist : 0.;
1230 } 1233 }
1231 case 2: // non-convex right prism 1234 case 2: // non-convex right prism
1232 { 1235 {
1233 G4bool in = PointInPolygon(p); 1236 G4bool in = PointInPolygon(p);
1234 if (in) 1237 if (in)
1235 { 1238 {
1236 G4double distz= std::max(fZSections[0 1239 G4double distz= std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
1237 return (distz > 0) ? distz : 0; 1240 return (distz > 0) ? distz : 0;
1238 } 1241 }
1239 else 1242 else
1240 { 1243 {
1241 G4double distz= std::max(fZSections[0 1244 G4double distz= std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
1242 G4double dd = DistanceToPolygonSqr(p) 1245 G4double dd = DistanceToPolygonSqr(p);
1243 if (distz > 0) dd += distz*distz; 1246 if (distz > 0) dd += distz*distz;
1244 return std::sqrt(dd); 1247 return std::sqrt(dd);
1245 } 1248 }
1246 } 1249 }
1247 } 1250 }
1248 1251
1249 // General case: use tessellated solid 1252 // General case: use tessellated solid
1250 return G4TessellatedSolid::DistanceToIn(p); 1253 return G4TessellatedSolid::DistanceToIn(p);
1251 } 1254 }
1252 1255
1253 //___________________________________________ 1256 //_____________________________________________________________________________
1254 1257
1255 G4double G4ExtrudedSolid::DistanceToOut (cons 1258 G4double G4ExtrudedSolid::DistanceToOut (const G4ThreeVector &p,
1256 cons 1259 const G4ThreeVector &v,
1257 cons 1260 const G4bool calcNorm,
1258 1261 G4bool* validNorm,
1259 1262 G4ThreeVector* n) const
1260 { 1263 {
1261 G4bool getnorm = calcNorm; 1264 G4bool getnorm = calcNorm;
1262 if (getnorm) *validNorm = true; 1265 if (getnorm) *validNorm = true;
1263 1266
1264 G4double z0 = fZSections[0].fZ; 1267 G4double z0 = fZSections[0].fZ;
1265 G4double z1 = fZSections[fNz-1].fZ; 1268 G4double z1 = fZSections[fNz-1].fZ;
1266 if ((p.z() <= z0 + kCarToleranceHalf) && v. 1269 if ((p.z() <= z0 + kCarToleranceHalf) && v.z() < 0)
1267 { 1270 {
1268 if (getnorm) n->set(0,0,-1); 1271 if (getnorm) n->set(0,0,-1);
1269 return 0; 1272 return 0;
1270 } 1273 }
1271 if ((p.z() >= z1 - kCarToleranceHalf) && v. 1274 if ((p.z() >= z1 - kCarToleranceHalf) && v.z() > 0)
1272 { 1275 {
1273 if (getnorm) n->set(0,0,1); 1276 if (getnorm) n->set(0,0,1);
1274 return 0; 1277 return 0;
1275 } 1278 }
1276 1279
1277 switch (fSolidType) 1280 switch (fSolidType)
1278 { 1281 {
1279 case 1: // convex right prism 1282 case 1: // convex right prism
1280 { 1283 {
1281 // Intersection with Z planes 1284 // Intersection with Z planes
1282 // 1285 //
1283 G4double dz = (z1 - z0)*0.5; 1286 G4double dz = (z1 - z0)*0.5;
1284 G4double pz = p.z() - 0.5 * (z0 + z1); 1287 G4double pz = p.z() - 0.5 * (z0 + z1);
1285 1288
1286 G4double vz = v.z(); 1289 G4double vz = v.z();
1287 G4double tmax = (vz == 0) ? DBL_MAX : ( 1290 G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(dz,vz) - pz)/vz;
1288 G4int iside = (vz < 0) ? -4 : -2; // li 1291 G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
1289 1292
1290 // Intersection with lateral planes 1293 // Intersection with lateral planes
1291 // 1294 //
1292 std::size_t np = fPlanes.size(); 1295 std::size_t np = fPlanes.size();
1293 for (std::size_t i=0; i<np; ++i) 1296 for (std::size_t i=0; i<np; ++i)
1294 { 1297 {
1295 G4double cosa = fPlanes[i].a*v.x()+fP 1298 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y();
1296 if (cosa > 0) 1299 if (cosa > 0)
1297 { 1300 {
1298 G4double dist = fPlanes[i].a*p.x()+ 1301 G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].d;
1299 if (dist >= -kCarToleranceHalf) 1302 if (dist >= -kCarToleranceHalf)
1300 { 1303 {
1301 if (getnorm) n->set(fPlanes[i].a, 1304 if (getnorm) n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c);
1302 return 0; 1305 return 0;
1303 } 1306 }
1304 G4double tmp = -dist/cosa; 1307 G4double tmp = -dist/cosa;
1305 if (tmax > tmp) { tmax = tmp; iside 1308 if (tmax > tmp) { tmax = tmp; iside = (G4int)i; }
1306 } 1309 }
1307 } 1310 }
1308 1311
1309 // Set normal, if required, and return 1312 // Set normal, if required, and return distance
1310 // 1313 //
1311 if (getnorm) 1314 if (getnorm)
1312 { 1315 {
1313 if (iside < 0) 1316 if (iside < 0)
1314 { n->set(0, 0, iside + 3); } // (-4 1317 { n->set(0, 0, iside + 3); } // (-4+3)=-1, (-2+3)=+1
1315 else 1318 else
1316 { n->set(fPlanes[iside].a, fPlanes[ 1319 { n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c); }
1317 } 1320 }
1318 return tmax; 1321 return tmax;
1319 } 1322 }
1320 case 2: // non-convex right prism 1323 case 2: // non-convex right prism
1321 { 1324 {
1322 } 1325 }
1323 } 1326 }
1324 1327
1325 // Override the base class function to rede 1328 // Override the base class function to redefine validNorm
1326 // (the solid can be concave) 1329 // (the solid can be concave)
1327 1330
1328 G4double distOut = 1331 G4double distOut =
1329 G4TessellatedSolid::DistanceToOut(p, v, c 1332 G4TessellatedSolid::DistanceToOut(p, v, calcNorm, validNorm, n);
1330 if (validNorm != nullptr) { *validNorm = fI << 1333 if (validNorm) { *validNorm = fIsConvex; }
1331 1334
1332 return distOut; 1335 return distOut;
1333 } 1336 }
1334 1337
1335 //___________________________________________ 1338 //_____________________________________________________________________________
1336 1339
1337 G4double G4ExtrudedSolid::DistanceToOut(const 1340 G4double G4ExtrudedSolid::DistanceToOut(const G4ThreeVector& p) const
1338 { 1341 {
1339 switch (fSolidType) 1342 switch (fSolidType)
1340 { 1343 {
1341 case 1: // convex right prism 1344 case 1: // convex right prism
1342 { 1345 {
1343 G4double dist = std::max(fZSections[0]. 1346 G4double dist = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
1344 std::size_t np = fPlanes.size(); 1347 std::size_t np = fPlanes.size();
1345 for (std::size_t i=0; i<np; ++i) 1348 for (std::size_t i=0; i<np; ++i)
1346 { 1349 {
1347 G4double dd = fPlanes[i].a*p.x() + fP 1350 G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
1348 if (dd > dist) dist = dd; 1351 if (dd > dist) dist = dd;
1349 } 1352 }
1350 return (dist < 0) ? -dist : 0.; 1353 return (dist < 0) ? -dist : 0.;
1351 } 1354 }
1352 case 2: // non-convex right prism 1355 case 2: // non-convex right prism
1353 { 1356 {
1354 G4double distz = std::max(fZSections[0] 1357 G4double distz = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
1355 G4bool in = PointInPolygon(p); 1358 G4bool in = PointInPolygon(p);
1356 if (distz >= 0 || (!in)) return 0; // p 1359 if (distz >= 0 || (!in)) return 0; // point is outside
1357 return std::min(-distz,std::sqrt(Distan 1360 return std::min(-distz,std::sqrt(DistanceToPolygonSqr(p)));
1358 } 1361 }
1359 } 1362 }
1360 1363
1361 // General case: use tessellated solid 1364 // General case: use tessellated solid
1362 return G4TessellatedSolid::DistanceToOut(p) 1365 return G4TessellatedSolid::DistanceToOut(p);
1363 } 1366 }
1364 1367
1365 //___________________________________________ 1368 //_____________________________________________________________________________
1366 // Get bounding box 1369 // Get bounding box
1367 1370
1368 void G4ExtrudedSolid::BoundingLimits(G4ThreeV 1371 void G4ExtrudedSolid::BoundingLimits(G4ThreeVector& pMin,
1369 G4ThreeV 1372 G4ThreeVector& pMax) const
1370 { 1373 {
1371 G4double xmin0 = kInfinity, xmax0 = -kInfin 1374 G4double xmin0 = kInfinity, xmax0 = -kInfinity;
1372 G4double ymin0 = kInfinity, ymax0 = -kInfin 1375 G4double ymin0 = kInfinity, ymax0 = -kInfinity;
1373 1376
1374 for (G4int i=0; i<GetNofVertices(); ++i) 1377 for (G4int i=0; i<GetNofVertices(); ++i)
1375 { 1378 {
1376 G4double x = fPolygon[i].x(); 1379 G4double x = fPolygon[i].x();
1377 if (x < xmin0) xmin0 = x; 1380 if (x < xmin0) xmin0 = x;
1378 if (x > xmax0) xmax0 = x; 1381 if (x > xmax0) xmax0 = x;
1379 G4double y = fPolygon[i].y(); 1382 G4double y = fPolygon[i].y();
1380 if (y < ymin0) ymin0 = y; 1383 if (y < ymin0) ymin0 = y;
1381 if (y > ymax0) ymax0 = y; 1384 if (y > ymax0) ymax0 = y;
1382 } 1385 }
1383 1386
1384 G4double xmin = kInfinity, xmax = -kInfinit 1387 G4double xmin = kInfinity, xmax = -kInfinity;
1385 G4double ymin = kInfinity, ymax = -kInfinit 1388 G4double ymin = kInfinity, ymax = -kInfinity;
1386 1389
1387 G4int nsect = GetNofZSections(); 1390 G4int nsect = GetNofZSections();
1388 for (G4int i=0; i<nsect; ++i) 1391 for (G4int i=0; i<nsect; ++i)
1389 { 1392 {
1390 ZSection zsect = GetZSection(i); 1393 ZSection zsect = GetZSection(i);
1391 G4double dx = zsect.fOffset.x(); 1394 G4double dx = zsect.fOffset.x();
1392 G4double dy = zsect.fOffset.y(); 1395 G4double dy = zsect.fOffset.y();
1393 G4double scale = zsect.fScale; 1396 G4double scale = zsect.fScale;
1394 xmin = std::min(xmin,xmin0*scale+dx); 1397 xmin = std::min(xmin,xmin0*scale+dx);
1395 xmax = std::max(xmax,xmax0*scale+dx); 1398 xmax = std::max(xmax,xmax0*scale+dx);
1396 ymin = std::min(ymin,ymin0*scale+dy); 1399 ymin = std::min(ymin,ymin0*scale+dy);
1397 ymax = std::max(ymax,ymax0*scale+dy); 1400 ymax = std::max(ymax,ymax0*scale+dy);
1398 } 1401 }
1399 1402
1400 G4double zmin = GetZSection(0).fZ; 1403 G4double zmin = GetZSection(0).fZ;
1401 G4double zmax = GetZSection(nsect-1).fZ; 1404 G4double zmax = GetZSection(nsect-1).fZ;
1402 1405
1403 pMin.set(xmin,ymin,zmin); 1406 pMin.set(xmin,ymin,zmin);
1404 pMax.set(xmax,ymax,zmax); 1407 pMax.set(xmax,ymax,zmax);
1405 1408
1406 // Check correctness of the bounding box 1409 // Check correctness of the bounding box
1407 // 1410 //
1408 if (pMin.x() >= pMax.x() || pMin.y() >= pMa 1411 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
1409 { 1412 {
1410 std::ostringstream message; 1413 std::ostringstream message;
1411 message << "Bad bounding box (min >= max) 1414 message << "Bad bounding box (min >= max) for solid: "
1412 << GetName() << " !" 1415 << GetName() << " !"
1413 << "\npMin = " << pMin 1416 << "\npMin = " << pMin
1414 << "\npMax = " << pMax; 1417 << "\npMax = " << pMax;
1415 G4Exception("G4ExtrudedSolid::BoundingLim 1418 G4Exception("G4ExtrudedSolid::BoundingLimits()",
1416 "GeomMgt0001", JustWarning, m 1419 "GeomMgt0001", JustWarning, message);
1417 DumpInfo(); 1420 DumpInfo();
1418 } 1421 }
1419 } 1422 }
1420 1423
1421 //___________________________________________ 1424 //_____________________________________________________________________________
1422 // Calculate extent under transform and speci 1425 // Calculate extent under transform and specified limit
1423 1426
1424 G4bool 1427 G4bool
1425 G4ExtrudedSolid::CalculateExtent(const EAxis 1428 G4ExtrudedSolid::CalculateExtent(const EAxis pAxis,
1426 const G4Voxe 1429 const G4VoxelLimits& pVoxelLimit,
1427 const G4Affi 1430 const G4AffineTransform& pTransform,
1428 G4doub 1431 G4double& pMin, G4double& pMax) const
1429 { 1432 {
1430 G4ThreeVector bmin, bmax; 1433 G4ThreeVector bmin, bmax;
1431 G4bool exist; 1434 G4bool exist;
1432 1435
1433 // Check bounding box (bbox) 1436 // Check bounding box (bbox)
1434 // 1437 //
1435 BoundingLimits(bmin,bmax); 1438 BoundingLimits(bmin,bmax);
1436 G4BoundingEnvelope bbox(bmin,bmax); 1439 G4BoundingEnvelope bbox(bmin,bmax);
1437 #ifdef G4BBOX_EXTENT 1440 #ifdef G4BBOX_EXTENT
1438 return bbox.CalculateExtent(pAxis,pVoxelLim 1441 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
1439 #endif 1442 #endif
1440 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVo 1443 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
1441 { 1444 {
1442 return exist = pMin < pMax; << 1445 return exist = (pMin < pMax) ? true : false;
1443 } 1446 }
1444 1447
1445 // To find the extent, the base polygon is 1448 // To find the extent, the base polygon is subdivided in triangles.
1446 // The extent is calculated as cumulative e 1449 // The extent is calculated as cumulative extent of the parts
1447 // formed by extrusion of the triangles 1450 // formed by extrusion of the triangles
1448 // 1451 //
1449 G4TwoVectorList triangles; 1452 G4TwoVectorList triangles;
1450 G4double eminlim = pVoxelLimit.GetMinExtent 1453 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
1451 G4double emaxlim = pVoxelLimit.GetMaxExtent 1454 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
1452 1455
1453 // triangulate the base polygon 1456 // triangulate the base polygon
1454 if (!G4GeomTools::TriangulatePolygon(fPolyg 1457 if (!G4GeomTools::TriangulatePolygon(fPolygon,triangles))
1455 { 1458 {
1456 std::ostringstream message; 1459 std::ostringstream message;
1457 message << "Triangulation of the base pol 1460 message << "Triangulation of the base polygon has failed for solid: "
1458 << GetName() << " !" 1461 << GetName() << " !"
1459 << "\nExtent has been calculated 1462 << "\nExtent has been calculated using boundary box";
1460 G4Exception("G4ExtrudedSolid::CalculateEx 1463 G4Exception("G4ExtrudedSolid::CalculateExtent()",
1461 "GeomMgt1002",JustWarning,mes 1464 "GeomMgt1002",JustWarning,message);
1462 return bbox.CalculateExtent(pAxis,pVoxelL 1465 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
1463 } 1466 }
1464 1467
1465 // allocate vector lists 1468 // allocate vector lists
1466 G4int nsect = GetNofZSections(); 1469 G4int nsect = GetNofZSections();
1467 std::vector<const G4ThreeVectorList *> poly 1470 std::vector<const G4ThreeVectorList *> polygons;
1468 polygons.resize(nsect); 1471 polygons.resize(nsect);
1469 for (G4int k=0; k<nsect; ++k) { polygons[k] 1472 for (G4int k=0; k<nsect; ++k) { polygons[k] = new G4ThreeVectorList(3); }
1470 1473
1471 // main loop along triangles 1474 // main loop along triangles
1472 pMin = kInfinity; 1475 pMin = kInfinity;
1473 pMax = -kInfinity; 1476 pMax = -kInfinity;
1474 G4int ntria = (G4int)triangles.size()/3; 1477 G4int ntria = (G4int)triangles.size()/3;
1475 for (G4int i=0; i<ntria; ++i) 1478 for (G4int i=0; i<ntria; ++i)
1476 { 1479 {
1477 G4int i3 = i*3; 1480 G4int i3 = i*3;
1478 for (G4int k=0; k<nsect; ++k) // extrude 1481 for (G4int k=0; k<nsect; ++k) // extrude triangle
1479 { 1482 {
1480 ZSection zsect = GetZSection(k); 1483 ZSection zsect = GetZSection(k);
1481 G4double z = zsect.fZ; 1484 G4double z = zsect.fZ;
1482 G4double dx = zsect.fOffset.x(); 1485 G4double dx = zsect.fOffset.x();
1483 G4double dy = zsect.fOffset.y(); 1486 G4double dy = zsect.fOffset.y();
1484 G4double scale = zsect.fScale; 1487 G4double scale = zsect.fScale;
1485 1488
1486 auto ptr = const_cast<G4ThreeVectorList << 1489 G4ThreeVectorList* ptr = const_cast<G4ThreeVectorList*>(polygons[k]);
1487 auto iter = ptr->begin(); << 1490 G4ThreeVectorList::iterator iter = ptr->begin();
1488 G4double x0 = triangles[i3+0].x()*scale 1491 G4double x0 = triangles[i3+0].x()*scale+dx;
1489 G4double y0 = triangles[i3+0].y()*scale 1492 G4double y0 = triangles[i3+0].y()*scale+dy;
1490 iter->set(x0,y0,z); 1493 iter->set(x0,y0,z);
1491 iter++; 1494 iter++;
1492 G4double x1 = triangles[i3+1].x()*scale 1495 G4double x1 = triangles[i3+1].x()*scale+dx;
1493 G4double y1 = triangles[i3+1].y()*scale 1496 G4double y1 = triangles[i3+1].y()*scale+dy;
1494 iter->set(x1,y1,z); 1497 iter->set(x1,y1,z);
1495 iter++; 1498 iter++;
1496 G4double x2 = triangles[i3+2].x()*scale 1499 G4double x2 = triangles[i3+2].x()*scale+dx;
1497 G4double y2 = triangles[i3+2].y()*scale 1500 G4double y2 = triangles[i3+2].y()*scale+dy;
1498 iter->set(x2,y2,z); 1501 iter->set(x2,y2,z);
1499 } 1502 }
1500 1503
1501 // set sub-envelope and adjust extent 1504 // set sub-envelope and adjust extent
1502 G4double emin,emax; 1505 G4double emin,emax;
1503 G4BoundingEnvelope benv(polygons); 1506 G4BoundingEnvelope benv(polygons);
1504 if (!benv.CalculateExtent(pAxis,pVoxelLim 1507 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
1505 if (emin < pMin) pMin = emin; 1508 if (emin < pMin) pMin = emin;
1506 if (emax > pMax) pMax = emax; 1509 if (emax > pMax) pMax = emax;
1507 if (eminlim > pMin && emaxlim < pMax) bre 1510 if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
1508 } 1511 }
1509 // free memory 1512 // free memory
1510 for (G4int k=0; k<nsect; ++k) { delete poly << 1513 for (G4int k=0; k<nsect; ++k) { delete polygons[k]; polygons[k]=0;}
1511 return (pMin < pMax); 1514 return (pMin < pMax);
1512 } 1515 }
1513 1516
1514 //___________________________________________ 1517 //_____________________________________________________________________________
1515 1518
1516 std::ostream& G4ExtrudedSolid::StreamInfo(std 1519 std::ostream& G4ExtrudedSolid::StreamInfo(std::ostream &os) const
1517 { 1520 {
1518 G4long oldprc = os.precision(16); 1521 G4long oldprc = os.precision(16);
1519 os << "------------------------------------ 1522 os << "-----------------------------------------------------------\n"
1520 << " *** Dump for solid - " << GetNam 1523 << " *** Dump for solid - " << GetName() << " ***\n"
1521 << " ================================ 1524 << " ===================================================\n"
1522 << " Solid geometry type: " << fGeometry 1525 << " Solid geometry type: " << fGeometryType << G4endl;
1523 1526
1524 if ( fIsConvex) 1527 if ( fIsConvex)
1525 { os << " Convex polygon; list of vertice 1528 { os << " Convex polygon; list of vertices:" << G4endl; }
1526 else 1529 else
1527 { os << " Concave polygon; list of vertic 1530 { os << " Concave polygon; list of vertices:" << G4endl; }
1528 1531
1529 for ( std::size_t i=0; i<fNv; ++i ) 1532 for ( std::size_t i=0; i<fNv; ++i )
1530 { 1533 {
1531 os << std::setw(5) << "#" << i 1534 os << std::setw(5) << "#" << i
1532 << " vx = " << fPolygon[i].x()/mm << 1535 << " vx = " << fPolygon[i].x()/mm << " mm"
1533 << " vy = " << fPolygon[i].y()/mm << 1536 << " vy = " << fPolygon[i].y()/mm << " mm" << G4endl;
1534 } 1537 }
1535 1538
1536 os << " Sections:" << G4endl; 1539 os << " Sections:" << G4endl;
1537 for ( std::size_t iz=0; iz<fNz; ++iz ) 1540 for ( std::size_t iz=0; iz<fNz; ++iz )
1538 { 1541 {
1539 os << " z = " << fZSections[iz].fZ/mm 1542 os << " z = " << fZSections[iz].fZ/mm << " mm "
1540 << " x0= " << fZSections[iz].fOffs 1543 << " x0= " << fZSections[iz].fOffset.x()/mm << " mm "
1541 << " y0= " << fZSections[iz].fOffs 1544 << " y0= " << fZSections[iz].fOffset.y()/mm << " mm "
1542 << " scale= " << fZSections[iz].fScal 1545 << " scale= " << fZSections[iz].fScale << G4endl;
1543 } 1546 }
1544 1547
1545 /* 1548 /*
1546 // Triangles (for debugging) 1549 // Triangles (for debugging)
1547 os << G4endl; 1550 os << G4endl;
1548 os << " Triangles:" << G4endl; 1551 os << " Triangles:" << G4endl;
1549 os << " Triangle # vertex1 vertex2 ve 1552 os << " Triangle # vertex1 vertex2 vertex3" << G4endl;
1550 1553
1551 G4int counter = 0; 1554 G4int counter = 0;
1552 std::vector< std::vector<G4int> >::const_it 1555 std::vector< std::vector<G4int> >::const_iterator it;
1553 for ( it = fTriangles.begin(); it != fTrian 1556 for ( it = fTriangles.begin(); it != fTriangles.end(); it++ ) {
1554 std::vector<G4int> triangle = *it; 1557 std::vector<G4int> triangle = *it;
1555 os << std::setw(10) << counter++ 1558 os << std::setw(10) << counter++
1556 << std::setw(10) << triangle[0] << st 1559 << std::setw(10) << triangle[0] << std::setw(10) << triangle[1]
1557 << std::setw(10) << triangle[2] 1560 << std::setw(10) << triangle[2]
1558 << G4endl; 1561 << G4endl;
1559 } 1562 }
1560 */ 1563 */
1561 os.precision(oldprc); 1564 os.precision(oldprc);
1562 1565
1563 return os; 1566 return os;
1564 } 1567 }
1565 1568
1566 #endif 1569 #endif
1567 1570