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