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