| 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 // G4 Polyhedron library 26 // G4 Polyhedron library
27 // 27 //
28 // History: 28 // History:
29 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@ce 29 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
30 // 30 //
31 // 30.09.96 E.Chernyaev 31 // 30.09.96 E.Chernyaev
32 // - added GetNextVertexIndex, GetVertex by Ya 32 // - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
33 // - added GetNextUnitNormal, GetNextEdgeIndic 33 // - added GetNextUnitNormal, GetNextEdgeIndices, GetNextEdge
34 // 34 //
35 // 15.12.96 E.Chernyaev 35 // 15.12.96 E.Chernyaev
36 // - added GetNumberOfRotationSteps, RotateEdg 36 // - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
37 // - rewritten G4PolyhedronCons; 37 // - rewritten G4PolyhedronCons;
38 // - added G4PolyhedronPara, ...Trap, ...Pgon, 38 // - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
39 // 39 //
40 // 01.06.97 E.Chernyaev 40 // 01.06.97 E.Chernyaev
41 // - modified RotateAroundZ, added SetSideFace 41 // - modified RotateAroundZ, added SetSideFacets
42 // 42 //
43 // 19.03.00 E.Chernyaev 43 // 19.03.00 E.Chernyaev
44 // - implemented boolean operations (add, subt 44 // - implemented boolean operations (add, subtract, intersect) on polyhedra;
45 // 45 //
46 // 25.05.01 E.Chernyaev 46 // 25.05.01 E.Chernyaev
47 // - added GetSurfaceArea() and GetVolume() 47 // - added GetSurfaceArea() and GetVolume()
48 // 48 //
49 // 05.11.02 E.Chernyaev 49 // 05.11.02 E.Chernyaev
50 // - added createTwistedTrap() and createPolyh 50 // - added createTwistedTrap() and createPolyhedron()
51 // 51 //
52 // 20.06.05 G.Cosmo 52 // 20.06.05 G.Cosmo
53 // - added HepPolyhedronEllipsoid 53 // - added HepPolyhedronEllipsoid
54 // 54 //
55 // 18.07.07 T.Nikitina 55 // 18.07.07 T.Nikitina
56 // - added HepPolyhedronParaboloid 56 // - added HepPolyhedronParaboloid
57 // 57 //
58 // 22.02.20 E.Chernyaev 58 // 22.02.20 E.Chernyaev
59 // - added HepPolyhedronTet, HepPolyhedronHybe 59 // - added HepPolyhedronTet, HepPolyhedronHyberbolicMirror
60 // 60 //
61 // 12.05.21 E.Chernyaev 61 // 12.05.21 E.Chernyaev
62 // - added TriangulatePolygon(), RotateContour 62 // - added TriangulatePolygon(), RotateContourAroundZ()
63 // - added HepPolyhedronPgon, HepPolyhedronPco 63 // - added HepPolyhedronPgon, HepPolyhedronPcon given by rz-contour
64 // 64 //
65 // 26.03.22 E.Chernyaev 65 // 26.03.22 E.Chernyaev
66 // - added SetVertex(), SetFacet() 66 // - added SetVertex(), SetFacet()
67 // - added HepPolyhedronTetMesh 67 // - added HepPolyhedronTetMesh
68 // 68 //
69 // 04.04.22 E.Chernyaev 69 // 04.04.22 E.Chernyaev
70 // - added JoinCoplanarFacets() 70 // - added JoinCoplanarFacets()
71 // 71 //
72 // 07.04.22 E.Chernyaev 72 // 07.04.22 E.Chernyaev
73 // - added HepPolyhedronBoxMesh 73 // - added HepPolyhedronBoxMesh
74 74
75 #include "HepPolyhedron.h" 75 #include "HepPolyhedron.h"
76 #include "G4PhysicalConstants.hh" 76 #include "G4PhysicalConstants.hh"
77 #include "G4Vector3D.hh" 77 #include "G4Vector3D.hh"
78 78
79 #include <cstdlib> // Required on some compil 79 #include <cstdlib> // Required on some compilers for std::abs(int) ...
80 #include <cmath> 80 #include <cmath>
81 #include <algorithm> 81 #include <algorithm>
82 82
83 using CLHEP::perMillion; 83 using CLHEP::perMillion;
84 using CLHEP::deg; 84 using CLHEP::deg;
85 using CLHEP::pi; 85 using CLHEP::pi;
86 using CLHEP::twopi; 86 using CLHEP::twopi;
87 using CLHEP::nm; 87 using CLHEP::nm;
88 const G4double spatialTolerance = 0.01*nm; 88 const G4double spatialTolerance = 0.01*nm;
89 89
90 /********************************************* 90 /***********************************************************************
91 * 91 * *
92 * Name: HepPolyhedron operator << 92 * Name: HepPolyhedron operator << Date: 09.05.96 *
93 * Author: E.Chernyaev (IHEP/Protvino) 93 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
94 * 94 * *
95 * Function: Print contents of G4 polyhedron 95 * Function: Print contents of G4 polyhedron *
96 * 96 * *
97 ********************************************* 97 ***********************************************************************/
98 std::ostream & operator<<(std::ostream & ostr, 98 std::ostream & operator<<(std::ostream & ostr, const G4Facet & facet) {
99 for (const auto& edge : facet.edge) { 99 for (const auto& edge : facet.edge) {
100 ostr << " " << edge.v << "/" << edge.f; 100 ostr << " " << edge.v << "/" << edge.f;
101 } 101 }
102 return ostr; 102 return ostr;
103 } 103 }
104 104
105 std::ostream & operator<<(std::ostream & ostr, 105 std::ostream & operator<<(std::ostream & ostr, const HepPolyhedron & ph) {
106 ostr << std::endl; 106 ostr << std::endl;
107 ostr << "Nvertices=" << ph.nvert << ", Nface 107 ostr << "Nvertices=" << ph.nvert << ", Nfacets=" << ph.nface << std::endl;
108 G4int i; 108 G4int i;
109 for (i=1; i<=ph.nvert; i++) { 109 for (i=1; i<=ph.nvert; i++) {
110 ostr << "xyz(" << i << ")=" 110 ostr << "xyz(" << i << ")="
111 << ph.pV[i].x() << ' ' << ph.pV[i].y 111 << ph.pV[i].x() << ' ' << ph.pV[i].y() << ' ' << ph.pV[i].z()
112 << std::endl; 112 << std::endl;
113 } 113 }
114 for (i=1; i<=ph.nface; i++) { 114 for (i=1; i<=ph.nface; i++) {
115 ostr << "face(" << i << ")=" << ph.pF[i] < 115 ostr << "face(" << i << ")=" << ph.pF[i] << std::endl;
116 } 116 }
117 return ostr; 117 return ostr;
118 } 118 }
119 119
120 HepPolyhedron::HepPolyhedron(G4int Nvert, G4in 120 HepPolyhedron::HepPolyhedron(G4int Nvert, G4int Nface)
121 /********************************************* 121 /***********************************************************************
122 * 122 * *
123 * Name: HepPolyhedron constructor with 123 * Name: HepPolyhedron constructor with Date: 26.03.2022 *
124 * allocation of memory 124 * allocation of memory Revised: *
125 * Author: E.Tcherniaev (E.Chernyaev) 125 * Author: E.Tcherniaev (E.Chernyaev) *
126 * 126 * *
127 ********************************************* 127 ***********************************************************************/
128 : nvert(0), nface(0), pV(nullptr), pF(nullptr) 128 : nvert(0), nface(0), pV(nullptr), pF(nullptr)
129 { 129 {
130 AllocateMemory(Nvert, Nface); 130 AllocateMemory(Nvert, Nface);
131 } 131 }
132 132
133 HepPolyhedron::HepPolyhedron(const HepPolyhedr 133 HepPolyhedron::HepPolyhedron(const HepPolyhedron &from)
134 /********************************************* 134 /***********************************************************************
135 * 135 * *
136 * Name: HepPolyhedron copy constructor 136 * Name: HepPolyhedron copy constructor Date: 23.07.96 *
137 * Author: E.Chernyaev (IHEP/Protvino) 137 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
138 * 138 * *
139 ********************************************* 139 ***********************************************************************/
140 : nvert(0), nface(0), pV(nullptr), pF(nullptr) 140 : nvert(0), nface(0), pV(nullptr), pF(nullptr)
141 { 141 {
142 AllocateMemory(from.nvert, from.nface); 142 AllocateMemory(from.nvert, from.nface);
143 for (G4int i=1; i<=nvert; i++) pV[i] = from. 143 for (G4int i=1; i<=nvert; i++) pV[i] = from.pV[i];
144 for (G4int k=1; k<=nface; k++) pF[k] = from. 144 for (G4int k=1; k<=nface; k++) pF[k] = from.pF[k];
145 } 145 }
146 146
147 HepPolyhedron::HepPolyhedron(HepPolyhedron&& f 147 HepPolyhedron::HepPolyhedron(HepPolyhedron&& from)
148 /********************************************* 148 /***********************************************************************
149 * 149 * *
150 * Name: HepPolyhedron move constructor 150 * Name: HepPolyhedron move constructor Date: 04.11.2019 *
151 * Author: E.Tcherniaev (E.Chernyaev) 151 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
152 * 152 * *
153 ********************************************* 153 ***********************************************************************/
154 : nvert(0), nface(0), pV(nullptr), pF(nullptr) 154 : nvert(0), nface(0), pV(nullptr), pF(nullptr)
155 { 155 {
156 nvert = from.nvert; 156 nvert = from.nvert;
157 nface = from.nface; 157 nface = from.nface;
158 pV = from.pV; 158 pV = from.pV;
159 pF = from.pF; 159 pF = from.pF;
160 160
161 // Release the data from the source object 161 // Release the data from the source object
162 from.nvert = 0; 162 from.nvert = 0;
163 from.nface = 0; 163 from.nface = 0;
164 from.pV = nullptr; 164 from.pV = nullptr;
165 from.pF = nullptr; 165 from.pF = nullptr;
166 } 166 }
167 167
168 HepPolyhedron & HepPolyhedron::operator=(const 168 HepPolyhedron & HepPolyhedron::operator=(const HepPolyhedron &from)
169 /********************************************* 169 /***********************************************************************
170 * 170 * *
171 * Name: HepPolyhedron operator = 171 * Name: HepPolyhedron operator = Date: 23.07.96 *
172 * Author: E.Chernyaev (IHEP/Protvino) 172 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
173 * 173 * *
174 * Function: Copy contents of one polyhedron t 174 * Function: Copy contents of one polyhedron to another *
175 * 175 * *
176 ********************************************* 176 ***********************************************************************/
177 { 177 {
178 if (this != &from) { 178 if (this != &from) {
179 AllocateMemory(from.nvert, from.nface); 179 AllocateMemory(from.nvert, from.nface);
180 for (G4int i=1; i<=nvert; i++) pV[i] = fro 180 for (G4int i=1; i<=nvert; i++) pV[i] = from.pV[i];
181 for (G4int k=1; k<=nface; k++) pF[k] = fro 181 for (G4int k=1; k<=nface; k++) pF[k] = from.pF[k];
182 } 182 }
183 return *this; 183 return *this;
184 } 184 }
185 185
186 HepPolyhedron & HepPolyhedron::operator=(HepPo 186 HepPolyhedron & HepPolyhedron::operator=(HepPolyhedron&& from)
187 /********************************************* 187 /***********************************************************************
188 * 188 * *
189 * Name: HepPolyhedron move operator = 189 * Name: HepPolyhedron move operator = Date: 04.11.2019 *
190 * Author: E.Tcherniaev (E.Chernyaev) 190 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
191 * 191 * *
192 * Function: Move contents of one polyhedron t 192 * Function: Move contents of one polyhedron to another *
193 * 193 * *
194 ********************************************* 194 ***********************************************************************/
195 { 195 {
196 if (this != &from) { 196 if (this != &from) {
197 delete [] pV; 197 delete [] pV;
198 delete [] pF; 198 delete [] pF;
199 nvert = from.nvert; 199 nvert = from.nvert;
200 nface = from.nface; 200 nface = from.nface;
201 pV = from.pV; 201 pV = from.pV;
202 pF = from.pF; 202 pF = from.pF;
203 203
204 // Release the data from the source object 204 // Release the data from the source object
205 from.nvert = 0; 205 from.nvert = 0;
206 from.nface = 0; 206 from.nface = 0;
207 from.pV = nullptr; 207 from.pV = nullptr;
208 from.pF = nullptr; 208 from.pF = nullptr;
209 } 209 }
210 return *this; 210 return *this;
211 } 211 }
212 212
213 G4int 213 G4int
214 HepPolyhedron::FindNeighbour(G4int iFace, G4in 214 HepPolyhedron::FindNeighbour(G4int iFace, G4int iNode, G4int iOrder) const
215 /********************************************* 215 /***********************************************************************
216 * 216 * *
217 * Name: HepPolyhedron::FindNeighbour 217 * Name: HepPolyhedron::FindNeighbour Date: 22.11.99 *
218 * Author: E.Chernyaev 218 * Author: E.Chernyaev Revised: *
219 * 219 * *
220 * Function: Find neighbouring face 220 * Function: Find neighbouring face *
221 * 221 * *
222 ********************************************* 222 ***********************************************************************/
223 { 223 {
224 G4int i; 224 G4int i;
225 for (i=0; i<4; i++) { 225 for (i=0; i<4; i++) {
226 if (iNode == std::abs(pF[iFace].edge[i].v) 226 if (iNode == std::abs(pF[iFace].edge[i].v)) break;
227 } 227 }
228 if (i == 4) { 228 if (i == 4) {
229 std::cerr 229 std::cerr
230 << "HepPolyhedron::FindNeighbour: face " 230 << "HepPolyhedron::FindNeighbour: face " << iFace
231 << " has no node " << iNode 231 << " has no node " << iNode
232 << std::endl; 232 << std::endl;
233 return 0; 233 return 0;
234 } 234 }
235 if (iOrder < 0) { 235 if (iOrder < 0) {
236 if ( --i < 0) i = 3; 236 if ( --i < 0) i = 3;
237 if (pF[iFace].edge[i].v == 0) i = 2; 237 if (pF[iFace].edge[i].v == 0) i = 2;
238 } 238 }
239 return (pF[iFace].edge[i].v > 0) ? 0 : pF[iF 239 return (pF[iFace].edge[i].v > 0) ? 0 : pF[iFace].edge[i].f;
240 } 240 }
241 241
242 G4Normal3D HepPolyhedron::FindNodeNormal(G4int 242 G4Normal3D HepPolyhedron::FindNodeNormal(G4int iFace, G4int iNode) const
243 /********************************************* 243 /***********************************************************************
244 * 244 * *
245 * Name: HepPolyhedron::FindNodeNormal 245 * Name: HepPolyhedron::FindNodeNormal Date: 22.11.99 *
246 * Author: E.Chernyaev 246 * Author: E.Chernyaev Revised: *
247 * 247 * *
248 * Function: Find normal at given node 248 * Function: Find normal at given node *
249 * 249 * *
250 ********************************************* 250 ***********************************************************************/
251 { 251 {
252 G4Normal3D normal = GetUnitNormal(iFace); 252 G4Normal3D normal = GetUnitNormal(iFace);
253 G4int k = iFace, iOrder = 1; 253 G4int k = iFace, iOrder = 1;
254 254
255 for(;;) { 255 for(;;) {
256 k = FindNeighbour(k, iNode, iOrder); 256 k = FindNeighbour(k, iNode, iOrder);
257 if (k == iFace) break; 257 if (k == iFace) break;
258 if (k > 0) { 258 if (k > 0) {
259 normal += GetUnitNormal(k); 259 normal += GetUnitNormal(k);
260 }else{ 260 }else{
261 if (iOrder < 0) break; 261 if (iOrder < 0) break;
262 k = iFace; 262 k = iFace;
263 iOrder = -iOrder; 263 iOrder = -iOrder;
264 } 264 }
265 } 265 }
266 return normal.unit(); 266 return normal.unit();
267 } 267 }
268 268
269 G4int HepPolyhedron::GetNumberOfRotationSteps( 269 G4int HepPolyhedron::GetNumberOfRotationSteps()
270 /********************************************* 270 /***********************************************************************
271 * 271 * *
272 * Name: HepPolyhedron::GetNumberOfRotationSte 272 * Name: HepPolyhedron::GetNumberOfRotationSteps Date: 24.06.97 *
273 * Author: J.Allison (Manchester University) 273 * Author: J.Allison (Manchester University) Revised: *
274 * 274 * *
275 * Function: Get number of steps for whole cir 275 * Function: Get number of steps for whole circle *
276 * 276 * *
277 ********************************************* 277 ***********************************************************************/
278 { 278 {
279 return fNumberOfRotationSteps; 279 return fNumberOfRotationSteps;
280 } 280 }
281 281
282 void HepPolyhedron::SetVertex(G4int index, con 282 void HepPolyhedron::SetVertex(G4int index, const G4Point3D& v)
283 /********************************************* 283 /***********************************************************************
284 * 284 * *
285 * Name: HepPolyhedron::SetVertex 285 * Name: HepPolyhedron::SetVertex Date: 26.03.22 *
286 * Author: E.Tcherniaev (E.Chernyaev) 286 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
287 * 287 * *
288 * Function: Set vertex 288 * Function: Set vertex *
289 * 289 * *
290 ********************************************* 290 ***********************************************************************/
291 { 291 {
292 if (index < 1 || index > nvert) 292 if (index < 1 || index > nvert)
293 { 293 {
294 std::cerr 294 std::cerr
295 << "HepPolyhedron::SetVertex: vertex ind 295 << "HepPolyhedron::SetVertex: vertex index = " << index
296 << " is out of range\n" 296 << " is out of range\n"
297 << " N. of vertices = " << nvert << "\ 297 << " N. of vertices = " << nvert << "\n"
298 << " N. of facets = " << nface << std: 298 << " N. of facets = " << nface << std::endl;
299 return; 299 return;
300 } 300 }
301 pV[index] = v; 301 pV[index] = v;
302 } 302 }
303 303
304 void 304 void
305 HepPolyhedron::SetFacet(G4int index, G4int iv1 305 HepPolyhedron::SetFacet(G4int index, G4int iv1, G4int iv2, G4int iv3, G4int iv4)
306 /********************************************* 306 /***********************************************************************
307 * 307 * *
308 * Name: HepPolyhedron::SetFacet 308 * Name: HepPolyhedron::SetFacet Date: 26.03.22 *
309 * Author: E.Tcherniaev (E.Chernyaev) 309 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
310 * 310 * *
311 * Function: Set facet 311 * Function: Set facet *
312 * 312 * *
313 ********************************************* 313 ***********************************************************************/
314 { 314 {
315 if (index < 1 || index > nface) 315 if (index < 1 || index > nface)
316 { 316 {
317 std::cerr 317 std::cerr
318 << "HepPolyhedron::SetFacet: facet index 318 << "HepPolyhedron::SetFacet: facet index = " << index
319 << " is out of range\n" 319 << " is out of range\n"
320 << " N. of vertices = " << nvert << "\ 320 << " N. of vertices = " << nvert << "\n"
321 << " N. of facets = " << nface << std: 321 << " N. of facets = " << nface << std::endl;
322 return; 322 return;
323 } 323 }
324 if (iv1 < 1 || iv1 > nvert || 324 if (iv1 < 1 || iv1 > nvert ||
325 iv2 < 1 || iv2 > nvert || 325 iv2 < 1 || iv2 > nvert ||
326 iv3 < 1 || iv3 > nvert || 326 iv3 < 1 || iv3 > nvert ||
327 iv4 < 0 || iv4 > nvert) 327 iv4 < 0 || iv4 > nvert)
328 { 328 {
329 std::cerr 329 std::cerr
330 << "HepPolyhedron::SetFacet: incorrectly 330 << "HepPolyhedron::SetFacet: incorrectly specified facet"
331 << " (" << iv1 << ", " << iv2 << ", " << 331 << " (" << iv1 << ", " << iv2 << ", " << iv3 << ", " << iv4 << ")\n"
332 << " N. of vertices = " << nvert << "\ 332 << " N. of vertices = " << nvert << "\n"
333 << " N. of facets = " << nface << std: 333 << " N. of facets = " << nface << std::endl;
334 return; 334 return;
335 } 335 }
336 pF[index] = G4Facet(iv1, 0, iv2, 0, iv3, 0, 336 pF[index] = G4Facet(iv1, 0, iv2, 0, iv3, 0, iv4, 0);
337 } 337 }
338 338
339 void HepPolyhedron::SetNumberOfRotationSteps(G 339 void HepPolyhedron::SetNumberOfRotationSteps(G4int n)
340 /********************************************* 340 /***********************************************************************
341 * 341 * *
342 * Name: HepPolyhedron::SetNumberOfRotationSte 342 * Name: HepPolyhedron::SetNumberOfRotationSteps Date: 24.06.97 *
343 * Author: J.Allison (Manchester University) 343 * Author: J.Allison (Manchester University) Revised: *
344 * 344 * *
345 * Function: Set number of steps for whole cir 345 * Function: Set number of steps for whole circle *
346 * 346 * *
347 ********************************************* 347 ***********************************************************************/
348 { 348 {
349 const G4int nMin = 3; 349 const G4int nMin = 3;
350 if (n < nMin) { 350 if (n < nMin) {
351 std::cerr 351 std::cerr
352 << "HepPolyhedron::SetNumberOfRotationSt 352 << "HepPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
353 << "number of steps per circle < " << nM 353 << "number of steps per circle < " << nMin << "; forced to " << nMin
354 << std::endl; 354 << std::endl;
355 fNumberOfRotationSteps = nMin; 355 fNumberOfRotationSteps = nMin;
356 }else{ 356 }else{
357 fNumberOfRotationSteps = n; 357 fNumberOfRotationSteps = n;
358 } 358 }
359 } 359 }
360 360
361 void HepPolyhedron::ResetNumberOfRotationSteps 361 void HepPolyhedron::ResetNumberOfRotationSteps()
362 /********************************************* 362 /***********************************************************************
363 * 363 * *
364 * Name: HepPolyhedron::GetNumberOfRotationSte 364 * Name: HepPolyhedron::GetNumberOfRotationSteps Date: 24.06.97 *
365 * Author: J.Allison (Manchester University) 365 * Author: J.Allison (Manchester University) Revised: *
366 * 366 * *
367 * Function: Reset number of steps for whole c 367 * Function: Reset number of steps for whole circle to default value *
368 * 368 * *
369 ********************************************* 369 ***********************************************************************/
370 { 370 {
371 fNumberOfRotationSteps = DEFAULT_NUMBER_OF_S 371 fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
372 } 372 }
373 373
374 void HepPolyhedron::AllocateMemory(G4int Nvert 374 void HepPolyhedron::AllocateMemory(G4int Nvert, G4int Nface)
375 /********************************************* 375 /***********************************************************************
376 * 376 * *
377 * Name: HepPolyhedron::AllocateMemory 377 * Name: HepPolyhedron::AllocateMemory Date: 19.06.96 *
378 * Author: E.Chernyaev (IHEP/Protvino) 378 * Author: E.Chernyaev (IHEP/Protvino) Revised: 05.11.02 *
379 * 379 * *
380 * Function: Allocate memory for GEANT4 polyhe 380 * Function: Allocate memory for GEANT4 polyhedron *
381 * 381 * *
382 * Input: Nvert - number of nodes 382 * Input: Nvert - number of nodes *
383 * Nface - number of faces 383 * Nface - number of faces *
384 * 384 * *
385 ********************************************* 385 ***********************************************************************/
386 { 386 {
387 if (nvert == Nvert && nface == Nface) return 387 if (nvert == Nvert && nface == Nface) return;
388 delete [] pV; 388 delete [] pV;
389 delete [] pF; 389 delete [] pF;
390 if (Nvert > 0 && Nface > 0) { 390 if (Nvert > 0 && Nface > 0) {
391 nvert = Nvert; 391 nvert = Nvert;
392 nface = Nface; 392 nface = Nface;
393 pV = new G4Point3D[nvert+1]; 393 pV = new G4Point3D[nvert+1];
394 pF = new G4Facet[nface+1]; 394 pF = new G4Facet[nface+1];
395 }else{ 395 }else{
396 nvert = 0; nface = 0; pV = nullptr; pF = n 396 nvert = 0; nface = 0; pV = nullptr; pF = nullptr;
397 } 397 }
398 } 398 }
399 399
400 void HepPolyhedron::CreatePrism() 400 void HepPolyhedron::CreatePrism()
401 /********************************************* 401 /***********************************************************************
402 * 402 * *
403 * Name: HepPolyhedron::CreatePrism 403 * Name: HepPolyhedron::CreatePrism Date: 15.07.96 *
404 * Author: E.Chernyaev (IHEP/Protvino) 404 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
405 * 405 * *
406 * Function: Set facets for a prism 406 * Function: Set facets for a prism *
407 * 407 * *
408 ********************************************* 408 ***********************************************************************/
409 { 409 {
410 enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRON 410 enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
411 411
412 pF[1] = G4Facet(1,LEFT, 4,BACK, 3,RIGHT, 412 pF[1] = G4Facet(1,LEFT, 4,BACK, 3,RIGHT, 2,FRONT);
413 pF[2] = G4Facet(5,TOP, 8,BACK, 4,BOTTOM, 413 pF[2] = G4Facet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
414 pF[3] = G4Facet(8,TOP, 7,RIGHT, 3,BOTTOM, 414 pF[3] = G4Facet(8,TOP, 7,RIGHT, 3,BOTTOM, 4,LEFT);
415 pF[4] = G4Facet(7,TOP, 6,FRONT, 2,BOTTOM, 415 pF[4] = G4Facet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
416 pF[5] = G4Facet(6,TOP, 5,LEFT, 1,BOTTOM, 416 pF[5] = G4Facet(6,TOP, 5,LEFT, 1,BOTTOM, 2,RIGHT);
417 pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK, 417 pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK, 8,LEFT);
418 } 418 }
419 419
420 void HepPolyhedron::RotateEdge(G4int k1, G4int 420 void HepPolyhedron::RotateEdge(G4int k1, G4int k2, G4double r1, G4double r2,
421 G4int v1, G4int 421 G4int v1, G4int v2, G4int vEdge,
422 G4bool ifWholeCi 422 G4bool ifWholeCircle, G4int nds, G4int &kface)
423 /********************************************* 423 /***********************************************************************
424 * 424 * *
425 * Name: HepPolyhedron::RotateEdge 425 * Name: HepPolyhedron::RotateEdge Date: 05.12.96 *
426 * Author: E.Chernyaev (IHEP/Protvino) 426 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
427 * 427 * *
428 * Function: Create set of facets by rotation 428 * Function: Create set of facets by rotation of an edge around Z-axis *
429 * 429 * *
430 * Input: k1, k2 - end vertices of the edge 430 * Input: k1, k2 - end vertices of the edge *
431 * r1, r2 - radiuses of the end vertice 431 * r1, r2 - radiuses of the end vertices *
432 * v1, v2 - visibility of edges produce 432 * v1, v2 - visibility of edges produced by rotation of the end *
433 * vertices 433 * vertices *
434 * vEdge - visibility of the edge 434 * vEdge - visibility of the edge *
435 * ifWholeCircle - is true in case of w 435 * ifWholeCircle - is true in case of whole circle rotation *
436 * nds - number of discrete steps 436 * nds - number of discrete steps *
437 * r[] - r-coordinates 437 * r[] - r-coordinates *
438 * kface - current free cell in the pF 438 * kface - current free cell in the pF array *
439 * 439 * *
440 ********************************************* 440 ***********************************************************************/
441 { 441 {
442 if (r1 == 0. && r2 == 0.) return; 442 if (r1 == 0. && r2 == 0.) return;
443 443
444 G4int i; 444 G4int i;
445 G4int i1 = k1; 445 G4int i1 = k1;
446 G4int i2 = k2; 446 G4int i2 = k2;
447 G4int ii1 = ifWholeCircle ? i1 : i1+nds; 447 G4int ii1 = ifWholeCircle ? i1 : i1+nds;
448 G4int ii2 = ifWholeCircle ? i2 : i2+nds; 448 G4int ii2 = ifWholeCircle ? i2 : i2+nds;
449 G4int vv = ifWholeCircle ? vEdge : 1; 449 G4int vv = ifWholeCircle ? vEdge : 1;
450 450
451 if (nds == 1) { 451 if (nds == 1) {
452 if (r1 == 0.) { 452 if (r1 == 0.) {
453 pF[kface++] = G4Facet(i1,0, v2*i2,0 453 pF[kface++] = G4Facet(i1,0, v2*i2,0, (i2+1),0);
454 }else if (r2 == 0.) { 454 }else if (r2 == 0.) {
455 pF[kface++] = G4Facet(i1,0, i2,0, 455 pF[kface++] = G4Facet(i1,0, i2,0, v1*(i1+1),0);
456 }else{ 456 }else{
457 pF[kface++] = G4Facet(i1,0, v2*i2,0 457 pF[kface++] = G4Facet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),0);
458 } 458 }
459 }else{ 459 }else{
460 if (r1 == 0.) { 460 if (r1 == 0.) {
461 pF[kface++] = G4Facet(vv*i1,0, v2*i 461 pF[kface++] = G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0);
462 for (i2++,i=1; i<nds-1; i2++,i++) { 462 for (i2++,i=1; i<nds-1; i2++,i++) {
463 pF[kface++] = G4Facet(vEdge*i1,0, v2*i 463 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
464 } 464 }
465 pF[kface++] = G4Facet(vEdge*i1,0, v2*i 465 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
466 }else if (r2 == 0.) { 466 }else if (r2 == 0.) {
467 pF[kface++] = G4Facet(vv*i1,0, vEdg 467 pF[kface++] = G4Facet(vv*i1,0, vEdge*i2,0, v1*(i1+1),0);
468 for (i1++,i=1; i<nds-1; i1++,i++) { 468 for (i1++,i=1; i<nds-1; i1++,i++) {
469 pF[kface++] = G4Facet(vEdge*i1,0, vEdg 469 pF[kface++] = G4Facet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
470 } 470 }
471 pF[kface++] = G4Facet(vEdge*i1,0, vv*i 471 pF[kface++] = G4Facet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
472 }else{ 472 }else{
473 pF[kface++] = G4Facet(vv*i1,0, v2*i 473 pF[kface++] = G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
474 for (i1++,i2++,i=1; i<nds-1; i1++,i2++,i 474 for (i1++,i2++,i=1; i<nds-1; i1++,i2++,i++) {
475 pF[kface++] = G4Facet(vEdge*i1,0, v2*i 475 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
476 } 476 }
477 pF[kface++] = G4Facet(vEdge*i1,0, v2*i 477 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
478 } 478 }
479 } 479 }
480 } 480 }
481 481
482 void HepPolyhedron::SetSideFacets(G4int ii[4], 482 void HepPolyhedron::SetSideFacets(G4int ii[4], G4int vv[4],
483 G4int *kk, G4 483 G4int *kk, G4double *r,
484 G4double dphi 484 G4double dphi, G4int nds, G4int &kface)
485 /********************************************* 485 /***********************************************************************
486 * 486 * *
487 * Name: HepPolyhedron::SetSideFacets 487 * Name: HepPolyhedron::SetSideFacets Date: 20.05.97 *
488 * Author: E.Chernyaev (IHEP/Protvino) 488 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
489 * 489 * *
490 * Function: Set side facets for the case of i 490 * Function: Set side facets for the case of incomplete rotation *
491 * 491 * *
492 * Input: ii[4] - indices of original vertices 492 * Input: ii[4] - indices of original vertices *
493 * vv[4] - visibility of edges 493 * vv[4] - visibility of edges *
494 * kk[] - indices of nodes 494 * kk[] - indices of nodes *
495 * r[] - radiuses 495 * r[] - radiuses *
496 * dphi - delta phi 496 * dphi - delta phi *
497 * nds - number of discrete steps 497 * nds - number of discrete steps *
498 * kface - current free cell in the pF 498 * kface - current free cell in the pF array *
499 * 499 * *
500 ********************************************* 500 ***********************************************************************/
501 { 501 {
502 G4int k1, k2, k3, k4; 502 G4int k1, k2, k3, k4;
503 503
504 if (std::abs(dphi-pi) < perMillion) { // hal 504 if (std::abs(dphi-pi) < perMillion) { // half a circle
505 for (G4int i=0; i<4; i++) { 505 for (G4int i=0; i<4; i++) {
506 k1 = ii[i]; 506 k1 = ii[i];
507 k2 = ii[(i+1)%4]; 507 k2 = ii[(i+1)%4];
508 if (r[k1] == 0. && r[k2] == 0.) vv[i] = 508 if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
509 } 509 }
510 } 510 }
511 511
512 if (ii[1] == ii[2]) { 512 if (ii[1] == ii[2]) {
513 k1 = kk[ii[0]]; 513 k1 = kk[ii[0]];
514 k2 = kk[ii[2]]; 514 k2 = kk[ii[2]];
515 k3 = kk[ii[3]]; 515 k3 = kk[ii[3]];
516 pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2 516 pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
517 if (r[ii[0]] != 0.) k1 += nds; 517 if (r[ii[0]] != 0.) k1 += nds;
518 if (r[ii[2]] != 0.) k2 += nds; 518 if (r[ii[2]] != 0.) k2 += nds;
519 if (r[ii[3]] != 0.) k3 += nds; 519 if (r[ii[3]] != 0.) k3 += nds;
520 pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2 520 pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
521 }else if (kk[ii[0]] == kk[ii[1]]) { 521 }else if (kk[ii[0]] == kk[ii[1]]) {
522 k1 = kk[ii[0]]; 522 k1 = kk[ii[0]];
523 k2 = kk[ii[2]]; 523 k2 = kk[ii[2]];
524 k3 = kk[ii[3]]; 524 k3 = kk[ii[3]];
525 pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2 525 pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
526 if (r[ii[0]] != 0.) k1 += nds; 526 if (r[ii[0]] != 0.) k1 += nds;
527 if (r[ii[2]] != 0.) k2 += nds; 527 if (r[ii[2]] != 0.) k2 += nds;
528 if (r[ii[3]] != 0.) k3 += nds; 528 if (r[ii[3]] != 0.) k3 += nds;
529 pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2 529 pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
530 }else if (kk[ii[2]] == kk[ii[3]]) { 530 }else if (kk[ii[2]] == kk[ii[3]]) {
531 k1 = kk[ii[0]]; 531 k1 = kk[ii[0]];
532 k2 = kk[ii[1]]; 532 k2 = kk[ii[1]];
533 k3 = kk[ii[2]]; 533 k3 = kk[ii[2]];
534 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2 534 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
535 if (r[ii[0]] != 0.) k1 += nds; 535 if (r[ii[0]] != 0.) k1 += nds;
536 if (r[ii[1]] != 0.) k2 += nds; 536 if (r[ii[1]] != 0.) k2 += nds;
537 if (r[ii[2]] != 0.) k3 += nds; 537 if (r[ii[2]] != 0.) k3 += nds;
538 pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2 538 pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
539 }else{ 539 }else{
540 k1 = kk[ii[0]]; 540 k1 = kk[ii[0]];
541 k2 = kk[ii[1]]; 541 k2 = kk[ii[1]];
542 k3 = kk[ii[2]]; 542 k3 = kk[ii[2]];
543 k4 = kk[ii[3]]; 543 k4 = kk[ii[3]];
544 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2 544 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
545 if (r[ii[0]] != 0.) k1 += nds; 545 if (r[ii[0]] != 0.) k1 += nds;
546 if (r[ii[1]] != 0.) k2 += nds; 546 if (r[ii[1]] != 0.) k2 += nds;
547 if (r[ii[2]] != 0.) k3 += nds; 547 if (r[ii[2]] != 0.) k3 += nds;
548 if (r[ii[3]] != 0.) k4 += nds; 548 if (r[ii[3]] != 0.) k4 += nds;
549 pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3 549 pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
550 } 550 }
551 } 551 }
552 552
553 void HepPolyhedron::RotateAroundZ(G4int nstep, 553 void HepPolyhedron::RotateAroundZ(G4int nstep, G4double phi, G4double dphi,
554 G4int np1, G4 554 G4int np1, G4int np2,
555 const G4doubl 555 const G4double *z, G4double *r,
556 G4int nodeVis 556 G4int nodeVis, G4int edgeVis)
557 /********************************************* 557 /***********************************************************************
558 * 558 * *
559 * Name: HepPolyhedron::RotateAroundZ 559 * Name: HepPolyhedron::RotateAroundZ Date: 27.11.96 *
560 * Author: E.Chernyaev (IHEP/Protvino) 560 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
561 * 561 * *
562 * Function: Create HepPolyhedron for a solid 562 * Function: Create HepPolyhedron for a solid produced by rotation of *
563 * two polylines around Z-axis 563 * two polylines around Z-axis *
564 * 564 * *
565 * Input: nstep - number of discrete steps, if 565 * Input: nstep - number of discrete steps, if 0 then default *
566 * phi - starting phi angle 566 * phi - starting phi angle *
567 * dphi - delta phi 567 * dphi - delta phi *
568 * np1 - number of points in external 568 * np1 - number of points in external polyline *
569 * (must be negative in case of 569 * (must be negative in case of closed polyline) *
570 * np2 - number of points in internal 570 * np2 - number of points in internal polyline (may be 1) *
571 * z[] - z-coordinates (+z >>> -z for 571 * z[] - z-coordinates (+z >>> -z for both polylines) *
572 * r[] - r-coordinates 572 * r[] - r-coordinates *
573 * nodeVis - how to Draw edges joing co 573 * nodeVis - how to Draw edges joing consecutive positions of *
574 * node during rotation 574 * node during rotation *
575 * edgeVis - how to Draw edges 575 * edgeVis - how to Draw edges *
576 * 576 * *
577 ********************************************* 577 ***********************************************************************/
578 { 578 {
579 static const G4double wholeCircle = twopi; 579 static const G4double wholeCircle = twopi;
580 580
581 // S E T R O T A T I O N P A R A M E T 581 // S E T R O T A T I O N P A R A M E T E R S
582 582
583 G4bool ifWholeCircle = std::abs(dphi-wholeCi 583 G4bool ifWholeCircle = std::abs(dphi-wholeCircle) < perMillion;
584 G4double delPhi = ifWholeCircle ? wholeCircl 584 G4double delPhi = ifWholeCircle ? wholeCircle : dphi;
585 G4int nSphi = nstep; 585 G4int nSphi = nstep;
586 if (nSphi <= 0) nSphi = GetNumberOfRotationS 586 if (nSphi <= 0) nSphi = GetNumberOfRotationSteps()*delPhi/wholeCircle + 0.5;
587 if (nSphi == 0) nSphi = 1; 587 if (nSphi == 0) nSphi = 1;
588 G4int nVphi = ifWholeCircle ? nSphi : nSphi 588 G4int nVphi = ifWholeCircle ? nSphi : nSphi + 1;
589 G4bool ifClosed = np1 <= 0; // true if exter 589 G4bool ifClosed = np1 <= 0; // true if external polyline is closed
590 590
591 // C O U N T V E R T I C E S 591 // C O U N T V E R T I C E S
592 592
593 G4int absNp1 = std::abs(np1); 593 G4int absNp1 = std::abs(np1);
594 G4int absNp2 = std::abs(np2); 594 G4int absNp2 = std::abs(np2);
595 G4int i1beg = 0; 595 G4int i1beg = 0;
596 G4int i1end = absNp1-1; 596 G4int i1end = absNp1-1;
597 G4int i2beg = absNp1; 597 G4int i2beg = absNp1;
598 G4int i2end = absNp1+absNp2-1; 598 G4int i2end = absNp1+absNp2-1;
599 G4int i, j, k; 599 G4int i, j, k;
600 600
601 for(i=i1beg; i<=i2end; i++) { 601 for(i=i1beg; i<=i2end; i++) {
602 if (std::abs(r[i]) < spatialTolerance) r[i 602 if (std::abs(r[i]) < spatialTolerance) r[i] = 0.;
603 } 603 }
604 604
605 // external polyline - check position of nod 605 // external polyline - check position of nodes relative to Z
606 // 606 //
607 G4int Nverts = 0; 607 G4int Nverts = 0;
608 for (i=i1beg; i<=i1end; i++) { 608 for (i=i1beg; i<=i1end; i++) {
609 Nverts += (r[i] == 0.) ? 1 : nVphi; 609 Nverts += (r[i] == 0.) ? 1 : nVphi;
610 } 610 }
611 611
612 // internal polyline 612 // internal polyline
613 // 613 //
614 G4bool ifSide1 = false; // whether to create 614 G4bool ifSide1 = false; // whether to create bottom faces
615 G4bool ifSide2 = false; // whether to create 615 G4bool ifSide2 = false; // whether to create top faces
616 616
617 if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1 617 if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) { // first node
618 Nverts += (r[i2beg] == 0.) ? 1 : nVphi; 618 Nverts += (r[i2beg] == 0.) ? 1 : nVphi;
619 ifSide1 = true; 619 ifSide1 = true;
620 } 620 }
621 621
622 for(i=i2beg+1; i<i2end; i++) { // intermedia 622 for(i=i2beg+1; i<i2end; i++) { // intermediate nodes
623 Nverts += (r[i] == 0.) ? 1 : nVphi; 623 Nverts += (r[i] == 0.) ? 1 : nVphi;
624 } 624 }
625 625
626 if (r[i2end] != r[i1end] || z[i2end] != z[i1 626 if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) { // last node
627 if (absNp2 > 1) Nverts += (r[i2end] == 0.) 627 if (absNp2 > 1) Nverts += (r[i2end] == 0.) ? 1 : nVphi;
628 ifSide2 = true; 628 ifSide2 = true;
629 } 629 }
630 630
631 // C O U N T F A C E S 631 // C O U N T F A C E S
632 632
633 // external lateral faces 633 // external lateral faces
634 // 634 //
635 G4int Nfaces = ifClosed ? absNp1*nSphi : (ab 635 G4int Nfaces = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi;
636 636
637 // internal lateral faces 637 // internal lateral faces
638 // 638 //
639 if (absNp2 > 1) { 639 if (absNp2 > 1) {
640 for(i=i2beg; i<i2end; i++) { 640 for(i=i2beg; i<i2end; i++) {
641 if (r[i] > 0. || r[i+1] > 0.) Nfaces += 641 if (r[i] > 0. || r[i+1] > 0.) Nfaces += nSphi;
642 } 642 }
643 643
644 if (ifClosed) { 644 if (ifClosed) {
645 if (r[i2end] > 0. || r[i2beg] > 0.) Nfac 645 if (r[i2end] > 0. || r[i2beg] > 0.) Nfaces += nSphi;
646 } 646 }
647 } 647 }
648 648
649 // bottom and top faces 649 // bottom and top faces
650 // 650 //
651 if (!ifClosed) { 651 if (!ifClosed) {
652 if (ifSide1 && (r[i1beg] > 0. || r[i2beg] 652 if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) Nfaces += nSphi;
653 if (ifSide2 && (r[i1end] > 0. || r[i2end] 653 if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) Nfaces += nSphi;
654 } 654 }
655 655
656 // phi_wedge faces 656 // phi_wedge faces
657 // 657 //
658 if (!ifWholeCircle) { 658 if (!ifWholeCircle) {
659 Nfaces += ifClosed ? 2*absNp1 : 2*(absNp1- 659 Nfaces += ifClosed ? 2*absNp1 : 2*(absNp1-1);
660 } 660 }
661 661
662 // A L L O C A T E M E M O R Y 662 // A L L O C A T E M E M O R Y
663 663
664 AllocateMemory(Nverts, Nfaces); 664 AllocateMemory(Nverts, Nfaces);
665 if (pV == nullptr || pF == nullptr) return; 665 if (pV == nullptr || pF == nullptr) return;
666 666
667 // G E N E R A T E V E R T I C E S 667 // G E N E R A T E V E R T I C E S
668 668
669 G4int *kk; // array of start indices along p 669 G4int *kk; // array of start indices along polylines
670 kk = new G4int[absNp1+absNp2]; 670 kk = new G4int[absNp1+absNp2];
671 671
672 // external polyline 672 // external polyline
673 // 673 //
674 k = 1; // free position in array of vertices 674 k = 1; // free position in array of vertices pV
675 for(i=i1beg; i<=i1end; i++) { 675 for(i=i1beg; i<=i1end; i++) {
676 kk[i] = k; 676 kk[i] = k;
677 if (r[i] == 0.) 677 if (r[i] == 0.)
678 { pV[k++] = G4Point3D(0, 0, z[i]); } else 678 { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
679 } 679 }
680 680
681 // first point of internal polyline 681 // first point of internal polyline
682 // 682 //
683 i = i2beg; 683 i = i2beg;
684 if (ifSide1) { 684 if (ifSide1) {
685 kk[i] = k; 685 kk[i] = k;
686 if (r[i] == 0.) 686 if (r[i] == 0.)
687 { pV[k++] = G4Point3D(0, 0, z[i]); } else 687 { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
688 }else{ 688 }else{
689 kk[i] = kk[i1beg]; 689 kk[i] = kk[i1beg];
690 } 690 }
691 691
692 // intermediate points of internal polyline 692 // intermediate points of internal polyline
693 // 693 //
694 for(i=i2beg+1; i<i2end; i++) { 694 for(i=i2beg+1; i<i2end; i++) {
695 kk[i] = k; 695 kk[i] = k;
696 if (r[i] == 0.) 696 if (r[i] == 0.)
697 { pV[k++] = G4Point3D(0, 0, z[i]); } else 697 { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
698 } 698 }
699 699
700 // last point of internal polyline 700 // last point of internal polyline
701 // 701 //
702 if (absNp2 > 1) { 702 if (absNp2 > 1) {
703 i = i2end; 703 i = i2end;
704 if (ifSide2) { 704 if (ifSide2) {
705 kk[i] = k; 705 kk[i] = k;
706 if (r[i] == 0.) pV[k] = G4Point3D(0, 0, 706 if (r[i] == 0.) pV[k] = G4Point3D(0, 0, z[i]);
707 }else{ 707 }else{
708 kk[i] = kk[i1end]; 708 kk[i] = kk[i1end];
709 } 709 }
710 } 710 }
711 711
712 // set vertices 712 // set vertices
713 // 713 //
714 G4double cosPhi, sinPhi; 714 G4double cosPhi, sinPhi;
715 715
716 for(j=0; j<nVphi; j++) { 716 for(j=0; j<nVphi; j++) {
717 cosPhi = std::cos(phi+j*delPhi/nSphi); 717 cosPhi = std::cos(phi+j*delPhi/nSphi);
718 sinPhi = std::sin(phi+j*delPhi/nSphi); 718 sinPhi = std::sin(phi+j*delPhi/nSphi);
719 for(i=i1beg; i<=i2end; i++) { 719 for(i=i1beg; i<=i2end; i++) {
720 if (r[i] != 0.) 720 if (r[i] != 0.)
721 pV[kk[i]+j] = G4Point3D(r[i]*cosPhi,r[ 721 pV[kk[i]+j] = G4Point3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
722 } 722 }
723 } 723 }
724 724
725 // G E N E R A T E F A C E S 725 // G E N E R A T E F A C E S
726 726
727 // external faces 727 // external faces
728 // 728 //
729 G4int v1,v2; 729 G4int v1,v2;
730 730
731 k = 1; // free position in array of faces pF 731 k = 1; // free position in array of faces pF
732 v2 = ifClosed ? nodeVis : 1; 732 v2 = ifClosed ? nodeVis : 1;
733 for(i=i1beg; i<i1end; i++) { 733 for(i=i1beg; i<i1end; i++) {
734 v1 = v2; 734 v1 = v2;
735 if (!ifClosed && i == i1end-1) { 735 if (!ifClosed && i == i1end-1) {
736 v2 = 1; 736 v2 = 1;
737 }else{ 737 }else{
738 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2] 738 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
739 } 739 }
740 RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v 740 RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
741 edgeVis, ifWholeCircle, nSphi, 741 edgeVis, ifWholeCircle, nSphi, k);
742 } 742 }
743 if (ifClosed) { 743 if (ifClosed) {
744 RotateEdge(kk[i1end], kk[i1beg], r[i1end], 744 RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
745 edgeVis, ifWholeCircle, nSphi, 745 edgeVis, ifWholeCircle, nSphi, k);
746 } 746 }
747 747
748 // internal faces 748 // internal faces
749 // 749 //
750 if (absNp2 > 1) { 750 if (absNp2 > 1) {
751 v2 = ifClosed ? nodeVis : 1; 751 v2 = ifClosed ? nodeVis : 1;
752 for(i=i2beg; i<i2end; i++) { 752 for(i=i2beg; i<i2end; i++) {
753 v1 = v2; 753 v1 = v2;
754 if (!ifClosed && i==i2end-1) { 754 if (!ifClosed && i==i2end-1) {
755 v2 = 1; 755 v2 = 1;
756 }else{ 756 }else{
757 v2 = (r[i] == r[i+1] && r[i+1] == r[i+ 757 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
758 } 758 }
759 RotateEdge(kk[i+1], kk[i], r[i+1], r[i], 759 RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
760 edgeVis, ifWholeCircle, nSphi 760 edgeVis, ifWholeCircle, nSphi, k);
761 } 761 }
762 if (ifClosed) { 762 if (ifClosed) {
763 RotateEdge(kk[i2beg], kk[i2end], r[i2beg 763 RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
764 edgeVis, ifWholeCircle, nSphi 764 edgeVis, ifWholeCircle, nSphi, k);
765 } 765 }
766 } 766 }
767 767
768 // bottom and top faces 768 // bottom and top faces
769 // 769 //
770 if (!ifClosed) { 770 if (!ifClosed) {
771 if (ifSide1) { 771 if (ifSide1) {
772 RotateEdge(kk[i2beg], kk[i1beg], r[i2beg 772 RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
773 -1, ifWholeCircle, nSphi, k); 773 -1, ifWholeCircle, nSphi, k);
774 } 774 }
775 if (ifSide2) { 775 if (ifSide2) {
776 RotateEdge(kk[i1end], kk[i2end], r[i1end 776 RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
777 -1, ifWholeCircle, nSphi, k); 777 -1, ifWholeCircle, nSphi, k);
778 } 778 }
779 } 779 }
780 780
781 // phi_wedge faces in case of incomplete cir 781 // phi_wedge faces in case of incomplete circle
782 // 782 //
783 if (!ifWholeCircle) { 783 if (!ifWholeCircle) {
784 784
785 G4int ii[4], vv[4]; 785 G4int ii[4], vv[4];
786 786
787 if (ifClosed) { 787 if (ifClosed) {
788 for (i=i1beg; i<=i1end; i++) { 788 for (i=i1beg; i<=i1end; i++) {
789 ii[0] = i; 789 ii[0] = i;
790 ii[3] = (i == i1end) ? i1beg : i+1; 790 ii[3] = (i == i1end) ? i1beg : i+1;
791 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+ 791 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
792 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+ 792 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
793 vv[0] = -1; 793 vv[0] = -1;
794 vv[1] = 1; 794 vv[1] = 1;
795 vv[2] = -1; 795 vv[2] = -1;
796 vv[3] = 1; 796 vv[3] = 1;
797 SetSideFacets(ii, vv, kk, r, delPhi, n 797 SetSideFacets(ii, vv, kk, r, delPhi, nSphi, k);
798 } 798 }
799 }else{ 799 }else{
800 for (i=i1beg; i<i1end; i++) { 800 for (i=i1beg; i<i1end; i++) {
801 ii[0] = i; 801 ii[0] = i;
802 ii[3] = i+1; 802 ii[3] = i+1;
803 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+ 803 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
804 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+ 804 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
805 vv[0] = (i == i1beg) ? 1 : -1; 805 vv[0] = (i == i1beg) ? 1 : -1;
806 vv[1] = 1; 806 vv[1] = 1;
807 vv[2] = (i == i1end-1) ? 1 : -1; 807 vv[2] = (i == i1end-1) ? 1 : -1;
808 vv[3] = 1; 808 vv[3] = 1;
809 SetSideFacets(ii, vv, kk, r, delPhi, n 809 SetSideFacets(ii, vv, kk, r, delPhi, nSphi, k);
810 } 810 }
811 } 811 }
812 } 812 }
813 813
814 delete [] kk; // free memory 814 delete [] kk; // free memory
815 815
816 // final check 816 // final check
817 // 817 //
818 if (k-1 != nface) { 818 if (k-1 != nface) {
819 std::cerr 819 std::cerr
820 << "HepPolyhedron::RotateAroundZ: number 820 << "HepPolyhedron::RotateAroundZ: number of generated faces ("
821 << k-1 << ") is not equal to the number 821 << k-1 << ") is not equal to the number of allocated faces ("
822 << nface << ")" 822 << nface << ")"
823 << std::endl; 823 << std::endl;
824 } 824 }
825 } 825 }
826 826
827 void 827 void
828 HepPolyhedron::RotateContourAroundZ(G4int nste 828 HepPolyhedron::RotateContourAroundZ(G4int nstep,
829 G4double p 829 G4double phi,
830 G4double d 830 G4double dphi,
831 const std: 831 const std::vector<G4TwoVector> &rz,
832 G4int node 832 G4int nodeVis,
833 G4int edge 833 G4int edgeVis)
834 /********************************************* 834 /***********************************************************************
835 * 835 * *
836 * Name: HepPolyhedron::RotateContourAroundZ 836 * Name: HepPolyhedron::RotateContourAroundZ Date: 12.05.21 *
837 * Author: E.Tcherniaev (E.Chernyaev) 837 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
838 * 838 * *
839 * Function: Create HepPolyhedron for a solid 839 * Function: Create HepPolyhedron for a solid produced by rotation of *
840 * a closed polyline (rz-contour) ar 840 * a closed polyline (rz-contour) around Z-axis *
841 * 841 * *
842 * Input: nstep - number of discrete steps, if 842 * Input: nstep - number of discrete steps, if 0 then default *
843 * phi - starting phi angle 843 * phi - starting phi angle *
844 * dphi - delta phi 844 * dphi - delta phi *
845 * rz - rz-contour 845 * rz - rz-contour *
846 * nodeVis - how to Draw edges joing co 846 * nodeVis - how to Draw edges joing consecutive positions of *
847 * node during rotation 847 * node during rotation *
848 * edgeVis - how to Draw edges 848 * edgeVis - how to Draw edges *
849 * 849 * *
850 ********************************************* 850 ***********************************************************************/
851 { 851 {
852 // S E T R O T A T I O N P A R A M E T 852 // S E T R O T A T I O N P A R A M E T E R S
853 853
854 G4bool ifWholeCircle = std::abs(dphi - twopi 854 G4bool ifWholeCircle = std::abs(dphi - twopi) < perMillion;
855 G4double delPhi = (ifWholeCircle) ? twopi : 855 G4double delPhi = (ifWholeCircle) ? twopi : dphi;
856 G4int nSphi = nstep; 856 G4int nSphi = nstep;
857 if (nSphi <= 0) nSphi = GetNumberOfRotationS 857 if (nSphi <= 0) nSphi = GetNumberOfRotationSteps()*delPhi/twopi + 0.5;
858 if (nSphi == 0) nSphi = 1; 858 if (nSphi == 0) nSphi = 1;
859 G4int nVphi = (ifWholeCircle) ? nSphi : nSph 859 G4int nVphi = (ifWholeCircle) ? nSphi : nSphi + 1;
860 860
861 // C A L C U L A T E A R E A 861 // C A L C U L A T E A R E A
862 862
863 G4int Nrz = (G4int)rz.size(); 863 G4int Nrz = (G4int)rz.size();
864 G4double area = 0; 864 G4double area = 0;
865 for (G4int i = 0; i < Nrz; ++i) 865 for (G4int i = 0; i < Nrz; ++i)
866 { 866 {
867 G4int k = (i == 0) ? Nrz - 1 : i - 1; 867 G4int k = (i == 0) ? Nrz - 1 : i - 1;
868 area += rz[k].x()*rz[i].y() - rz[i].x()*rz 868 area += rz[k].x()*rz[i].y() - rz[i].x()*rz[k].y();
869 } 869 }
870 870
871 // P R E P A R E P O L Y L I N E 871 // P R E P A R E P O L Y L I N E
872 872
873 auto r = new G4double[Nrz]; 873 auto r = new G4double[Nrz];
874 auto z = new G4double[Nrz]; 874 auto z = new G4double[Nrz];
875 for (G4int i = 0; i < Nrz; ++i) 875 for (G4int i = 0; i < Nrz; ++i)
876 { 876 {
877 r[i] = rz[i].x(); 877 r[i] = rz[i].x();
878 z[i] = rz[i].y(); 878 z[i] = rz[i].y();
879 if (std::abs(r[i]) < spatialTolerance) r[i 879 if (std::abs(r[i]) < spatialTolerance) r[i] = 0.;
880 } 880 }
881 881
882 // C O U N T V E R T I C E S A N D F 882 // C O U N T V E R T I C E S A N D F A C E S
883 883
884 G4int Nverts = 0; 884 G4int Nverts = 0;
885 for(G4int i = 0; i < Nrz; ++i) Nverts += (r[ 885 for(G4int i = 0; i < Nrz; ++i) Nverts += (r[i] == 0.) ? 1 : nVphi;
886 886
887 G4int Nedges = Nrz; 887 G4int Nedges = Nrz;
888 for (G4int i = 0; i < Nrz; ++i) 888 for (G4int i = 0; i < Nrz; ++i)
889 { 889 {
890 G4int k = (i == 0) ? Nrz - 1 : i - 1; 890 G4int k = (i == 0) ? Nrz - 1 : i - 1;
891 Nedges -= static_cast<int>(r[k] == 0 && r[ 891 Nedges -= static_cast<int>(r[k] == 0 && r[i] == 0);
892 } 892 }
893 893
894 G4int Nfaces = Nedges*nSphi; / 894 G4int Nfaces = Nedges*nSphi; // lateral faces
895 if (!ifWholeCircle) Nfaces += 2*(Nrz - 2); / 895 if (!ifWholeCircle) Nfaces += 2*(Nrz - 2); // phi_wedge faces
896 896
897 // A L L O C A T E M E M O R Y 897 // A L L O C A T E M E M O R Y
898 898
899 AllocateMemory(Nverts, Nfaces); 899 AllocateMemory(Nverts, Nfaces);
900 if (pV == nullptr || pF == nullptr) 900 if (pV == nullptr || pF == nullptr)
901 { 901 {
902 delete [] r; 902 delete [] r;
903 delete [] z; 903 delete [] z;
904 return; 904 return;
905 } 905 }
906 906
907 // S E T V E R T I C E S 907 // S E T V E R T I C E S
908 908
909 auto kk = new G4int[Nrz]; // start indices a 909 auto kk = new G4int[Nrz]; // start indices along contour
910 G4int kfree = 1; // current free position in 910 G4int kfree = 1; // current free position in array of vertices pV
911 911
912 // set start indices, set vertices for nodes 912 // set start indices, set vertices for nodes with r == 0
913 for(G4int i = 0; i < Nrz; ++i) 913 for(G4int i = 0; i < Nrz; ++i)
914 { 914 {
915 kk[i] = kfree; 915 kk[i] = kfree;
916 if (r[i] == 0.) pV[kfree++] = G4Point3D(0, 916 if (r[i] == 0.) pV[kfree++] = G4Point3D(0, 0, z[i]);
917 if (r[i] != 0.) kfree += nVphi; 917 if (r[i] != 0.) kfree += nVphi;
918 } 918 }
919 919
920 // set vertices by rotating r 920 // set vertices by rotating r
921 for(G4int j = 0; j < nVphi; ++j) 921 for(G4int j = 0; j < nVphi; ++j)
922 { 922 {
923 G4double cosPhi = std::cos(phi + j*delPhi/ 923 G4double cosPhi = std::cos(phi + j*delPhi/nSphi);
924 G4double sinPhi = std::sin(phi + j*delPhi/ 924 G4double sinPhi = std::sin(phi + j*delPhi/nSphi);
925 for(G4int i = 0; i < Nrz; ++i) 925 for(G4int i = 0; i < Nrz; ++i)
926 { 926 {
927 if (r[i] != 0.) 927 if (r[i] != 0.)
928 pV[kk[i] + j] = G4Point3D(r[i]*cosPhi, 928 pV[kk[i] + j] = G4Point3D(r[i]*cosPhi, r[i]*sinPhi, z[i]);
929 } 929 }
930 } 930 }
931 931
932 // S E T F A C E S 932 // S E T F A C E S
933 933
934 kfree = 1; // current free position in array 934 kfree = 1; // current free position in array of faces pF
935 for(G4int i = 0; i < Nrz; ++i) 935 for(G4int i = 0; i < Nrz; ++i)
936 { 936 {
937 G4int i1 = (i < Nrz - 1) ? i + 1 : 0; // i 937 G4int i1 = (i < Nrz - 1) ? i + 1 : 0; // inverse order if area > 0
938 G4int i2 = i; 938 G4int i2 = i;
939 if (area < 0.) std::swap(i1, i2); 939 if (area < 0.) std::swap(i1, i2);
940 RotateEdge(kk[i1], kk[i2], r[i1], r[i2], n 940 RotateEdge(kk[i1], kk[i2], r[i1], r[i2], nodeVis, nodeVis,
941 edgeVis, ifWholeCircle, nSphi, 941 edgeVis, ifWholeCircle, nSphi, kfree);
942 } 942 }
943 943
944 // S E T P H I _ W E D G E F A C E S 944 // S E T P H I _ W E D G E F A C E S
945 945
946 if (!ifWholeCircle) 946 if (!ifWholeCircle)
947 { 947 {
948 std::vector<G4int> triangles; 948 std::vector<G4int> triangles;
949 TriangulatePolygon(rz, triangles); 949 TriangulatePolygon(rz, triangles);
950 950
951 G4int ii[4], vv[4]; 951 G4int ii[4], vv[4];
952 G4int ntria = G4int(triangles.size()/3); 952 G4int ntria = G4int(triangles.size()/3);
953 for (G4int i = 0; i < ntria; ++i) 953 for (G4int i = 0; i < ntria; ++i)
954 { 954 {
955 G4int i1 = triangles[0 + i*3]; 955 G4int i1 = triangles[0 + i*3];
956 G4int i2 = triangles[1 + i*3]; 956 G4int i2 = triangles[1 + i*3];
957 G4int i3 = triangles[2 + i*3]; 957 G4int i3 = triangles[2 + i*3];
958 if (area < 0.) std::swap(i1, i3); 958 if (area < 0.) std::swap(i1, i3);
959 G4int v1 = (std::abs(i2-i1) == 1 || std: 959 G4int v1 = (std::abs(i2-i1) == 1 || std::abs(i2-i1) == Nrz-1) ? 1 : -1;
960 G4int v2 = (std::abs(i3-i2) == 1 || std: 960 G4int v2 = (std::abs(i3-i2) == 1 || std::abs(i3-i2) == Nrz-1) ? 1 : -1;
961 G4int v3 = (std::abs(i1-i3) == 1 || std: 961 G4int v3 = (std::abs(i1-i3) == 1 || std::abs(i1-i3) == Nrz-1) ? 1 : -1;
962 ii[0] = i1; ii[1] = i2; ii[2] = i2; ii[3 962 ii[0] = i1; ii[1] = i2; ii[2] = i2; ii[3] = i3;
963 vv[0] = v1; vv[1] = -1; vv[2] = v2; vv[3 963 vv[0] = v1; vv[1] = -1; vv[2] = v2; vv[3] = v3;
964 SetSideFacets(ii, vv, kk, r, delPhi, nSp 964 SetSideFacets(ii, vv, kk, r, delPhi, nSphi, kfree);
965 } 965 }
966 } 966 }
967 967
968 // free memory 968 // free memory
969 delete [] r; 969 delete [] r;
970 delete [] z; 970 delete [] z;
971 delete [] kk; 971 delete [] kk;
972 972
973 // final check 973 // final check
974 if (kfree - 1 != nface) 974 if (kfree - 1 != nface)
975 { 975 {
976 std::cerr 976 std::cerr
977 << "HepPolyhedron::RotateContourAroundZ: 977 << "HepPolyhedron::RotateContourAroundZ: number of generated faces ("
978 << kfree-1 << ") is not equal to the num 978 << kfree-1 << ") is not equal to the number of allocated faces ("
979 << nface << ")" 979 << nface << ")"
980 << std::endl; 980 << std::endl;
981 } 981 }
982 } 982 }
983 983
984 G4bool 984 G4bool
985 HepPolyhedron::TriangulatePolygon(const std::v 985 HepPolyhedron::TriangulatePolygon(const std::vector<G4TwoVector> &polygon,
986 std::vector< 986 std::vector<G4int> &result)
987 /********************************************* 987 /***********************************************************************
988 * 988 * *
989 * Name: HepPolyhedron::TriangulatePolygon 989 * Name: HepPolyhedron::TriangulatePolygon Date: 12.05.21 *
990 * Author: E.Tcherniaev (E.Chernyaev) 990 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
991 * 991 * *
992 * Function: Simple implementation of "ear cli 992 * Function: Simple implementation of "ear clipping" algorithm for *
993 * triangulation of a simple contour 993 * triangulation of a simple contour/polygon, it places *
994 * the result in a std::vector as tr 994 * the result in a std::vector as triplets of vertex indices *
995 * 995 * *
996 * If triangulation is sucsessfull t 996 * If triangulation is sucsessfull then the function *
997 * returns true, otherwise false 997 * returns true, otherwise false *
998 * 998 * *
999 * Remark: It's a copy of G4GeomTools::Trian 999 * Remark: It's a copy of G4GeomTools::TriangulatePolygon() *
1000 * 1000 * *
1001 ******************************************** 1001 ***********************************************************************/
1002 { 1002 {
1003 result.resize(0); 1003 result.resize(0);
1004 G4int n = (G4int)polygon.size(); 1004 G4int n = (G4int)polygon.size();
1005 if (n < 3) return false; 1005 if (n < 3) return false;
1006 1006
1007 // calculate area 1007 // calculate area
1008 // 1008 //
1009 G4double area = 0.; 1009 G4double area = 0.;
1010 for(G4int i = 0; i < n; ++i) 1010 for(G4int i = 0; i < n; ++i)
1011 { 1011 {
1012 G4int k = (i == 0) ? n - 1 : i - 1; 1012 G4int k = (i == 0) ? n - 1 : i - 1;
1013 area += polygon[k].x()*polygon[i].y() - p 1013 area += polygon[k].x()*polygon[i].y() - polygon[i].x()*polygon[k].y();
1014 } 1014 }
1015 1015
1016 // allocate and initialize list of Vertices 1016 // allocate and initialize list of Vertices
1017 // we want a counter-clockwise polygon in V 1017 // we want a counter-clockwise polygon in V
1018 // 1018 //
1019 auto V = new G4int[n]; 1019 auto V = new G4int[n];
1020 if (area > 0.) 1020 if (area > 0.)
1021 for (G4int i = 0; i < n; ++i) V[i] = i; 1021 for (G4int i = 0; i < n; ++i) V[i] = i;
1022 else 1022 else
1023 for (G4int i = 0; i < n; ++i) V[i] = (n - 1023 for (G4int i = 0; i < n; ++i) V[i] = (n - 1) - i;
1024 1024
1025 // Triangulation: remove nv-2 Vertices, cr 1025 // Triangulation: remove nv-2 Vertices, creating 1 triangle every time
1026 // 1026 //
1027 G4int nv = n; 1027 G4int nv = n;
1028 G4int count = 2*nv; // error detection coun 1028 G4int count = 2*nv; // error detection counter
1029 for(G4int b = nv - 1; nv > 2; ) 1029 for(G4int b = nv - 1; nv > 2; )
1030 { 1030 {
1031 // ERROR: if we loop, it is probably a no 1031 // ERROR: if we loop, it is probably a non-simple polygon
1032 if ((count--) <= 0) 1032 if ((count--) <= 0)
1033 { 1033 {
1034 delete [] V; 1034 delete [] V;
1035 if (area < 0.) std::reverse(result.begi 1035 if (area < 0.) std::reverse(result.begin(),result.end());
1036 return false; 1036 return false;
1037 } 1037 }
1038 1038
1039 // three consecutive vertices in current 1039 // three consecutive vertices in current polygon, <a,b,c>
1040 G4int a = (b < nv) ? b : 0; // previo 1040 G4int a = (b < nv) ? b : 0; // previous
1041 b = (a+1 < nv) ? a+1 : 0; // curren 1041 b = (a+1 < nv) ? a+1 : 0; // current
1042 G4int c = (b+1 < nv) ? b+1 : 0; // next 1042 G4int c = (b+1 < nv) ? b+1 : 0; // next
1043 1043
1044 if (CheckSnip(polygon, a,b,c, nv,V)) 1044 if (CheckSnip(polygon, a,b,c, nv,V))
1045 { 1045 {
1046 // output Triangle 1046 // output Triangle
1047 result.push_back(V[a]); 1047 result.push_back(V[a]);
1048 result.push_back(V[b]); 1048 result.push_back(V[b]);
1049 result.push_back(V[c]); 1049 result.push_back(V[c]);
1050 1050
1051 // remove vertex b from remaining polyg 1051 // remove vertex b from remaining polygon
1052 nv--; 1052 nv--;
1053 for(G4int i = b; i < nv; ++i) V[i] = V[ 1053 for(G4int i = b; i < nv; ++i) V[i] = V[i+1];
1054 1054
1055 count = 2*nv; // resest error detection 1055 count = 2*nv; // resest error detection counter
1056 } 1056 }
1057 } 1057 }
1058 delete [] V; 1058 delete [] V;
1059 if (area < 0.) std::reverse(result.begin(), 1059 if (area < 0.) std::reverse(result.begin(),result.end());
1060 return true; 1060 return true;
1061 } 1061 }
1062 1062
1063 G4bool HepPolyhedron::CheckSnip(const std::ve 1063 G4bool HepPolyhedron::CheckSnip(const std::vector<G4TwoVector> &contour,
1064 G4int a, G4in 1064 G4int a, G4int b, G4int c,
1065 G4int n, cons 1065 G4int n, const G4int* V)
1066 /******************************************** 1066 /***********************************************************************
1067 * 1067 * *
1068 * Name: HepPolyhedron::CheckSnip 1068 * Name: HepPolyhedron::CheckSnip Date: 12.05.21 *
1069 * Author: E.Tcherniaev (E.Chernyaev) 1069 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
1070 * 1070 * *
1071 * Function: Check for a valid snip, 1071 * Function: Check for a valid snip, *
1072 * it is a helper functionfor Trian 1072 * it is a helper functionfor TriangulatePolygon() *
1073 * 1073 * *
1074 ******************************************** 1074 ***********************************************************************/
1075 { 1075 {
1076 static const G4double kCarTolerance = 1.e-9 1076 static const G4double kCarTolerance = 1.e-9;
1077 1077
1078 // check orientation of Triangle 1078 // check orientation of Triangle
1079 G4double Ax = contour[V[a]].x(), Ay = conto 1079 G4double Ax = contour[V[a]].x(), Ay = contour[V[a]].y();
1080 G4double Bx = contour[V[b]].x(), By = conto 1080 G4double Bx = contour[V[b]].x(), By = contour[V[b]].y();
1081 G4double Cx = contour[V[c]].x(), Cy = conto 1081 G4double Cx = contour[V[c]].x(), Cy = contour[V[c]].y();
1082 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) < kCa 1082 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) < kCarTolerance) return false;
1083 1083
1084 // check that there is no point inside Tria 1084 // check that there is no point inside Triangle
1085 G4double xmin = std::min(std::min(Ax,Bx),Cx 1085 G4double xmin = std::min(std::min(Ax,Bx),Cx);
1086 G4double xmax = std::max(std::max(Ax,Bx),Cx 1086 G4double xmax = std::max(std::max(Ax,Bx),Cx);
1087 G4double ymin = std::min(std::min(Ay,By),Cy 1087 G4double ymin = std::min(std::min(Ay,By),Cy);
1088 G4double ymax = std::max(std::max(Ay,By),Cy 1088 G4double ymax = std::max(std::max(Ay,By),Cy);
1089 1089
1090 for (G4int i=0; i<n; ++i) 1090 for (G4int i=0; i<n; ++i)
1091 { 1091 {
1092 if((i == a) || (i == b) || (i == c)) cont 1092 if((i == a) || (i == b) || (i == c)) continue;
1093 G4double Px = contour[V[i]].x(); 1093 G4double Px = contour[V[i]].x();
1094 if (Px < xmin || Px > xmax) continue; 1094 if (Px < xmin || Px > xmax) continue;
1095 G4double Py = contour[V[i]].y(); 1095 G4double Py = contour[V[i]].y();
1096 if (Py < ymin || Py > ymax) continue; 1096 if (Py < ymin || Py > ymax) continue;
1097 // if (PointInTriangle(Ax,Ay,Bx,By,Cx,Cy, 1097 // if (PointInTriangle(Ax,Ay,Bx,By,Cx,Cy,Px,Py)) return false;
1098 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0 1098 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0.)
1099 { 1099 {
1100 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < 1100 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < 0.) continue;
1101 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < 1101 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < 0.) continue;
1102 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < 1102 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < 0.) continue;
1103 } 1103 }
1104 else 1104 else
1105 { 1105 {
1106 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > 1106 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > 0.) continue;
1107 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > 1107 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > 0.) continue;
1108 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > 1108 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > 0.) continue;
1109 } 1109 }
1110 return false; 1110 return false;
1111 } 1111 }
1112 return true; 1112 return true;
1113 } 1113 }
1114 1114
1115 void HepPolyhedron::SetReferences() 1115 void HepPolyhedron::SetReferences()
1116 /******************************************** 1116 /***********************************************************************
1117 * 1117 * *
1118 * Name: HepPolyhedron::SetReferences 1118 * Name: HepPolyhedron::SetReferences Date: 04.12.96 *
1119 * Author: E.Chernyaev (IHEP/Protvino) 1119 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1120 * 1120 * *
1121 * Function: For each edge set reference to n 1121 * Function: For each edge set reference to neighbouring facet *
1122 * 1122 * *
1123 ******************************************** 1123 ***********************************************************************/
1124 { 1124 {
1125 if (nface <= 0) return; 1125 if (nface <= 0) return;
1126 1126
1127 struct edgeListMember { 1127 struct edgeListMember {
1128 edgeListMember *next; 1128 edgeListMember *next;
1129 G4int v2; 1129 G4int v2;
1130 G4int iface; 1130 G4int iface;
1131 G4int iedge; 1131 G4int iedge;
1132 } *edgeList, *freeList, **headList; 1132 } *edgeList, *freeList, **headList;
1133 1133
1134 1134
1135 // A L L O C A T E A N D I N I T I A 1135 // A L L O C A T E A N D I N I T I A T E L I S T S
1136 1136
1137 edgeList = new edgeListMember[2*nface]; 1137 edgeList = new edgeListMember[2*nface];
1138 headList = new edgeListMember*[nvert]; 1138 headList = new edgeListMember*[nvert];
1139 1139
1140 G4int i; 1140 G4int i;
1141 for (i=0; i<nvert; i++) { 1141 for (i=0; i<nvert; i++) {
1142 headList[i] = nullptr; 1142 headList[i] = nullptr;
1143 } 1143 }
1144 freeList = edgeList; 1144 freeList = edgeList;
1145 for (i=0; i<2*nface-1; i++) { 1145 for (i=0; i<2*nface-1; i++) {
1146 edgeList[i].next = &edgeList[i+1]; 1146 edgeList[i].next = &edgeList[i+1];
1147 } 1147 }
1148 edgeList[2*nface-1].next = nullptr; 1148 edgeList[2*nface-1].next = nullptr;
1149 1149
1150 // L O O P A L O N G E D G E S 1150 // L O O P A L O N G E D G E S
1151 1151
1152 G4int iface, iedge, nedge, i1, i2, k1, k2; 1152 G4int iface, iedge, nedge, i1, i2, k1, k2;
1153 edgeListMember *prev, *cur; 1153 edgeListMember *prev, *cur;
1154 1154
1155 for(iface=1; iface<=nface; iface++) { 1155 for(iface=1; iface<=nface; iface++) {
1156 nedge = (pF[iface].edge[3].v == 0) ? 3 : 1156 nedge = (pF[iface].edge[3].v == 0) ? 3 : 4;
1157 for (iedge=0; iedge<nedge; iedge++) { 1157 for (iedge=0; iedge<nedge; iedge++) {
1158 i1 = iedge; 1158 i1 = iedge;
1159 i2 = (iedge < nedge-1) ? iedge+1 : 0; 1159 i2 = (iedge < nedge-1) ? iedge+1 : 0;
1160 i1 = std::abs(pF[iface].edge[i1].v); 1160 i1 = std::abs(pF[iface].edge[i1].v);
1161 i2 = std::abs(pF[iface].edge[i2].v); 1161 i2 = std::abs(pF[iface].edge[i2].v);
1162 k1 = (i1 < i2) ? i1 : i2; // k 1162 k1 = (i1 < i2) ? i1 : i2; // k1 = ::min(i1,i2);
1163 k2 = (i1 > i2) ? i1 : i2; // k 1163 k2 = (i1 > i2) ? i1 : i2; // k2 = ::max(i1,i2);
1164 1164
1165 // check head of the List corresponding 1165 // check head of the List corresponding to k1
1166 cur = headList[k1]; 1166 cur = headList[k1];
1167 if (cur == nullptr) { 1167 if (cur == nullptr) {
1168 headList[k1] = freeList; 1168 headList[k1] = freeList;
1169 if (freeList == nullptr) { 1169 if (freeList == nullptr) {
1170 std::cerr 1170 std::cerr
1171 << "Polyhedron::SetReferences: bad 1171 << "Polyhedron::SetReferences: bad link "
1172 << std::endl; 1172 << std::endl;
1173 break; 1173 break;
1174 } 1174 }
1175 freeList = freeList->next; 1175 freeList = freeList->next;
1176 cur = headList[k1]; 1176 cur = headList[k1];
1177 cur->next = nullptr; 1177 cur->next = nullptr;
1178 cur->v2 = k2; 1178 cur->v2 = k2;
1179 cur->iface = iface; 1179 cur->iface = iface;
1180 cur->iedge = iedge; 1180 cur->iedge = iedge;
1181 continue; 1181 continue;
1182 } 1182 }
1183 1183
1184 if (cur->v2 == k2) { 1184 if (cur->v2 == k2) {
1185 headList[k1] = cur->next; 1185 headList[k1] = cur->next;
1186 cur->next = freeList; 1186 cur->next = freeList;
1187 freeList = cur; 1187 freeList = cur;
1188 pF[iface].edge[iedge].f = cur->iface; 1188 pF[iface].edge[iedge].f = cur->iface;
1189 pF[cur->iface].edge[cur->iedge].f = i 1189 pF[cur->iface].edge[cur->iedge].f = iface;
1190 i1 = (pF[iface].edge[iedge].v < 0) ? 1190 i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
1191 i2 = (pF[cur->iface].edge[cur->iedge] 1191 i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
1192 if (i1 != i2) { 1192 if (i1 != i2) {
1193 std::cerr 1193 std::cerr
1194 << "Polyhedron::SetReferences: di 1194 << "Polyhedron::SetReferences: different edge visibility "
1195 << iface << "/" << iedge << "/" 1195 << iface << "/" << iedge << "/"
1196 << pF[iface].edge[iedge].v << " a 1196 << pF[iface].edge[iedge].v << " and "
1197 << cur->iface << "/" << cur->iedg 1197 << cur->iface << "/" << cur->iedge << "/"
1198 << pF[cur->iface].edge[cur->iedge 1198 << pF[cur->iface].edge[cur->iedge].v
1199 << std::endl; 1199 << std::endl;
1200 } 1200 }
1201 continue; 1201 continue;
1202 } 1202 }
1203 1203
1204 // check List itself 1204 // check List itself
1205 for (;;) { 1205 for (;;) {
1206 prev = cur; 1206 prev = cur;
1207 cur = prev->next; 1207 cur = prev->next;
1208 if (cur == nullptr) { 1208 if (cur == nullptr) {
1209 prev->next = freeList; 1209 prev->next = freeList;
1210 if (freeList == nullptr) { 1210 if (freeList == nullptr) {
1211 std::cerr 1211 std::cerr
1212 << "Polyhedron::SetReferences: ba 1212 << "Polyhedron::SetReferences: bad link "
1213 << std::endl; 1213 << std::endl;
1214 break; 1214 break;
1215 } 1215 }
1216 freeList = freeList->next; 1216 freeList = freeList->next;
1217 cur = prev->next; 1217 cur = prev->next;
1218 cur->next = nullptr; 1218 cur->next = nullptr;
1219 cur->v2 = k2; 1219 cur->v2 = k2;
1220 cur->iface = iface; 1220 cur->iface = iface;
1221 cur->iedge = iedge; 1221 cur->iedge = iedge;
1222 break; 1222 break;
1223 } 1223 }
1224 1224
1225 if (cur->v2 == k2) { 1225 if (cur->v2 == k2) {
1226 prev->next = cur->next; 1226 prev->next = cur->next;
1227 cur->next = freeList; 1227 cur->next = freeList;
1228 freeList = cur; 1228 freeList = cur;
1229 pF[iface].edge[iedge].f = cur->ifac 1229 pF[iface].edge[iedge].f = cur->iface;
1230 pF[cur->iface].edge[cur->iedge].f = 1230 pF[cur->iface].edge[cur->iedge].f = iface;
1231 i1 = (pF[iface].edge[iedge].v < 0) 1231 i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
1232 i2 = (pF[cur->iface].edge[cur->iedg 1232 i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
1233 if (i1 != i2) { 1233 if (i1 != i2) {
1234 std::cerr 1234 std::cerr
1235 << "Polyhedron::SetReferences 1235 << "Polyhedron::SetReferences: different edge visibility "
1236 << iface << "/" << iedge << " 1236 << iface << "/" << iedge << "/"
1237 << pF[iface].edge[iedge].v << 1237 << pF[iface].edge[iedge].v << " and "
1238 << cur->iface << "/" << cur-> 1238 << cur->iface << "/" << cur->iedge << "/"
1239 << pF[cur->iface].edge[cur->i 1239 << pF[cur->iface].edge[cur->iedge].v
1240 << std::endl; 1240 << std::endl;
1241 } 1241 }
1242 break; 1242 break;
1243 } 1243 }
1244 } 1244 }
1245 } 1245 }
1246 } 1246 }
1247 1247
1248 // C H E C K T H A T A L L L I S T S 1248 // C H E C K T H A T A L L L I S T S A R E E M P T Y
1249 1249
1250 for (i=0; i<nvert; i++) { 1250 for (i=0; i<nvert; i++) {
1251 if (headList[i] != nullptr) { 1251 if (headList[i] != nullptr) {
1252 std::cerr 1252 std::cerr
1253 << "Polyhedron::SetReferences: List " 1253 << "Polyhedron::SetReferences: List " << i << " is not empty"
1254 << std::endl; 1254 << std::endl;
1255 } 1255 }
1256 } 1256 }
1257 1257
1258 // F R E E M E M O R Y 1258 // F R E E M E M O R Y
1259 1259
1260 delete [] edgeList; 1260 delete [] edgeList;
1261 delete [] headList; 1261 delete [] headList;
1262 } 1262 }
1263 1263
1264 void HepPolyhedron::JoinCoplanarFacets(G4doub 1264 void HepPolyhedron::JoinCoplanarFacets(G4double tolerance)
1265 /******************************************** 1265 /***********************************************************************
1266 * 1266 * *
1267 * Name: HepPolyhedron::JoinCoplanarFacets 1267 * Name: HepPolyhedron::JoinCoplanarFacets Date: 04.04.22 *
1268 * Author: E.Tcherniaev (E.Chernyaev) 1268 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
1269 * 1269 * *
1270 * Function: Join couples of triangular facet 1270 * Function: Join couples of triangular facets to quadrangular facets *
1271 * where it is possible 1271 * where it is possible *
1272 * 1272 * *
1273 ******************************************** 1273 ***********************************************************************/
1274 { 1274 {
1275 G4int njoin = 0; 1275 G4int njoin = 0;
1276 for (G4int icur = 1; icur <= nface; ++icur) 1276 for (G4int icur = 1; icur <= nface; ++icur)
1277 { 1277 {
1278 // skip if already joined or quadrangle 1278 // skip if already joined or quadrangle
1279 if (pF[icur].edge[0].v == 0) continue; 1279 if (pF[icur].edge[0].v == 0) continue;
1280 if (pF[icur].edge[3].v != 0) continue; 1280 if (pF[icur].edge[3].v != 0) continue;
1281 // skip if all references point to alread 1281 // skip if all references point to already checked facets
1282 if (pF[icur].edge[0].f < icur && 1282 if (pF[icur].edge[0].f < icur &&
1283 pF[icur].edge[1].f < icur && 1283 pF[icur].edge[1].f < icur &&
1284 pF[icur].edge[2].f < icur) continue; 1284 pF[icur].edge[2].f < icur) continue;
1285 // compute plane equation 1285 // compute plane equation
1286 G4Normal3D norm = GetUnitNormal(icur); 1286 G4Normal3D norm = GetUnitNormal(icur);
1287 G4double dd = norm.dot(pV[pF[icur].edge[0 1287 G4double dd = norm.dot(pV[pF[icur].edge[0].v]);
1288 G4int vcur0 = std::abs(pF[icur].edge[0].v 1288 G4int vcur0 = std::abs(pF[icur].edge[0].v);
1289 G4int vcur1 = std::abs(pF[icur].edge[1].v 1289 G4int vcur1 = std::abs(pF[icur].edge[1].v);
1290 G4int vcur2 = std::abs(pF[icur].edge[2].v 1290 G4int vcur2 = std::abs(pF[icur].edge[2].v);
1291 // select neighbouring facet 1291 // select neighbouring facet
1292 G4int kcheck = 0, icheck = 0, vcheck = 0; 1292 G4int kcheck = 0, icheck = 0, vcheck = 0;
1293 G4double dist = DBL_MAX; 1293 G4double dist = DBL_MAX;
1294 for (G4int k = 0; k < 3; ++k) 1294 for (G4int k = 0; k < 3; ++k)
1295 { 1295 {
1296 G4int itmp = pF[icur].edge[k].f; 1296 G4int itmp = pF[icur].edge[k].f;
1297 // skip if already checked, joined or q 1297 // skip if already checked, joined or quadrangle
1298 if (itmp < icur) continue; 1298 if (itmp < icur) continue;
1299 if (pF[itmp].edge[0].v == 0 || 1299 if (pF[itmp].edge[0].v == 0 ||
1300 pF[itmp].edge[3].v != 0) continue; 1300 pF[itmp].edge[3].v != 0) continue;
1301 // get candidate vertex 1301 // get candidate vertex
1302 G4int vtmp = 0; 1302 G4int vtmp = 0;
1303 for (G4int j = 0; j < 3; ++j) 1303 for (G4int j = 0; j < 3; ++j)
1304 { 1304 {
1305 vtmp = std::abs(pF[itmp].edge[j].v); 1305 vtmp = std::abs(pF[itmp].edge[j].v);
1306 if (vtmp != vcur0 && vtmp != vcur1 && vtmp 1306 if (vtmp != vcur0 && vtmp != vcur1 && vtmp != vcur2) break;
1307 } 1307 }
1308 // check distance to the plane 1308 // check distance to the plane
1309 G4double dtmp = std::abs(norm.dot(pV[vt 1309 G4double dtmp = std::abs(norm.dot(pV[vtmp]) - dd);
1310 if (dtmp > tolerance || dtmp >= dist) c 1310 if (dtmp > tolerance || dtmp >= dist) continue;
1311 dist = dtmp; 1311 dist = dtmp;
1312 kcheck = k; 1312 kcheck = k;
1313 icheck = itmp; 1313 icheck = itmp;
1314 vcheck = vtmp; 1314 vcheck = vtmp;
1315 } 1315 }
1316 if (icheck == 0) continue; // no facet se 1316 if (icheck == 0) continue; // no facet selected
1317 // join facets 1317 // join facets
1318 njoin++; 1318 njoin++;
1319 pF[icheck].edge[0].v = 0; // mark facet a 1319 pF[icheck].edge[0].v = 0; // mark facet as joined
1320 if (kcheck == 0) 1320 if (kcheck == 0)
1321 { 1321 {
1322 pF[icur].edge[3].v = pF[icur].edge[2].v 1322 pF[icur].edge[3].v = pF[icur].edge[2].v;
1323 pF[icur].edge[2].v = pF[icur].edge[1].v 1323 pF[icur].edge[2].v = pF[icur].edge[1].v;
1324 pF[icur].edge[1].v = vcheck; 1324 pF[icur].edge[1].v = vcheck;
1325 } 1325 }
1326 else if (kcheck == 1) 1326 else if (kcheck == 1)
1327 { 1327 {
1328 pF[icur].edge[3].v = pF[icur].edge[2].v 1328 pF[icur].edge[3].v = pF[icur].edge[2].v;
1329 pF[icur].edge[2].v = vcheck; 1329 pF[icur].edge[2].v = vcheck;
1330 } 1330 }
1331 else 1331 else
1332 { 1332 {
1333 pF[icur].edge[3].v = vcheck; 1333 pF[icur].edge[3].v = vcheck;
1334 } 1334 }
1335 } 1335 }
1336 if (njoin == 0) return; // no joined facets 1336 if (njoin == 0) return; // no joined facets
1337 1337
1338 // restructure facets 1338 // restructure facets
1339 G4int nnew = 0; 1339 G4int nnew = 0;
1340 for (G4int icur = 1; icur <= nface; ++icur) 1340 for (G4int icur = 1; icur <= nface; ++icur)
1341 { 1341 {
1342 if (pF[icur].edge[0].v == 0) continue; 1342 if (pF[icur].edge[0].v == 0) continue;
1343 nnew++; 1343 nnew++;
1344 pF[nnew].edge[0].v = pF[icur].edge[0].v; 1344 pF[nnew].edge[0].v = pF[icur].edge[0].v;
1345 pF[nnew].edge[1].v = pF[icur].edge[1].v; 1345 pF[nnew].edge[1].v = pF[icur].edge[1].v;
1346 pF[nnew].edge[2].v = pF[icur].edge[2].v; 1346 pF[nnew].edge[2].v = pF[icur].edge[2].v;
1347 pF[nnew].edge[3].v = pF[icur].edge[3].v; 1347 pF[nnew].edge[3].v = pF[icur].edge[3].v;
1348 } 1348 }
1349 nface = nnew; 1349 nface = nnew;
1350 SetReferences(); 1350 SetReferences();
1351 } 1351 }
1352 1352
1353 void HepPolyhedron::InvertFacets() 1353 void HepPolyhedron::InvertFacets()
1354 /******************************************** 1354 /***********************************************************************
1355 * 1355 * *
1356 * Name: HepPolyhedron::InvertFacets 1356 * Name: HepPolyhedron::InvertFacets Date: 01.12.99 *
1357 * Author: E.Chernyaev 1357 * Author: E.Chernyaev Revised: *
1358 * 1358 * *
1359 * Function: Invert the order of the nodes in 1359 * Function: Invert the order of the nodes in the facets *
1360 * 1360 * *
1361 ******************************************** 1361 ***********************************************************************/
1362 { 1362 {
1363 if (nface <= 0) return; 1363 if (nface <= 0) return;
1364 G4int i, k, nnode, v[4],f[4]; 1364 G4int i, k, nnode, v[4],f[4];
1365 for (i=1; i<=nface; i++) { 1365 for (i=1; i<=nface; i++) {
1366 nnode = (pF[i].edge[3].v == 0) ? 3 : 4; 1366 nnode = (pF[i].edge[3].v == 0) ? 3 : 4;
1367 for (k=0; k<nnode; k++) { 1367 for (k=0; k<nnode; k++) {
1368 v[k] = (k+1 == nnode) ? pF[i].edge[0].v 1368 v[k] = (k+1 == nnode) ? pF[i].edge[0].v : pF[i].edge[k+1].v;
1369 if (v[k] * pF[i].edge[k].v < 0) v[k] = 1369 if (v[k] * pF[i].edge[k].v < 0) v[k] = -v[k];
1370 f[k] = pF[i].edge[k].f; 1370 f[k] = pF[i].edge[k].f;
1371 } 1371 }
1372 for (k=0; k<nnode; k++) { 1372 for (k=0; k<nnode; k++) {
1373 pF[i].edge[nnode-1-k].v = v[k]; 1373 pF[i].edge[nnode-1-k].v = v[k];
1374 pF[i].edge[nnode-1-k].f = f[k]; 1374 pF[i].edge[nnode-1-k].f = f[k];
1375 } 1375 }
1376 } 1376 }
1377 } 1377 }
1378 1378
1379 HepPolyhedron & HepPolyhedron::Transform(cons 1379 HepPolyhedron & HepPolyhedron::Transform(const G4Transform3D &t)
1380 /******************************************** 1380 /***********************************************************************
1381 * 1381 * *
1382 * Name: HepPolyhedron::Transform 1382 * Name: HepPolyhedron::Transform Date: 01.12.99 *
1383 * Author: E.Chernyaev 1383 * Author: E.Chernyaev Revised: *
1384 * 1384 * *
1385 * Function: Make transformation of the polyh 1385 * Function: Make transformation of the polyhedron *
1386 * 1386 * *
1387 ******************************************** 1387 ***********************************************************************/
1388 { 1388 {
1389 if (nvert > 0) { 1389 if (nvert > 0) {
1390 for (G4int i=1; i<=nvert; i++) { pV[i] = 1390 for (G4int i=1; i<=nvert; i++) { pV[i] = t * pV[i]; }
1391 1391
1392 // C H E C K D E T E R M I N A N T A 1392 // C H E C K D E T E R M I N A N T A N D
1393 // I N V E R T F A C E T S I F I T 1393 // I N V E R T F A C E T S I F I T I S N E G A T I V E
1394 1394
1395 G4Vector3D d = t * G4Vector3D(0,0,0); 1395 G4Vector3D d = t * G4Vector3D(0,0,0);
1396 G4Vector3D x = t * G4Vector3D(1,0,0) - d; 1396 G4Vector3D x = t * G4Vector3D(1,0,0) - d;
1397 G4Vector3D y = t * G4Vector3D(0,1,0) - d; 1397 G4Vector3D y = t * G4Vector3D(0,1,0) - d;
1398 G4Vector3D z = t * G4Vector3D(0,0,1) - d; 1398 G4Vector3D z = t * G4Vector3D(0,0,1) - d;
1399 if ((x.cross(y))*z < 0) InvertFacets(); 1399 if ((x.cross(y))*z < 0) InvertFacets();
1400 } 1400 }
1401 return *this; 1401 return *this;
1402 } 1402 }
1403 1403
1404 G4bool HepPolyhedron::GetNextVertexIndex(G4in 1404 G4bool HepPolyhedron::GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
1405 /******************************************** 1405 /***********************************************************************
1406 * 1406 * *
1407 * Name: HepPolyhedron::GetNextVertexIndex 1407 * Name: HepPolyhedron::GetNextVertexIndex Date: 03.09.96 *
1408 * Author: Yasuhide Sawada 1408 * Author: Yasuhide Sawada Revised: *
1409 * 1409 * *
1410 * Function: 1410 * Function: *
1411 * 1411 * *
1412 ******************************************** 1412 ***********************************************************************/
1413 { 1413 {
1414 static G4ThreadLocal G4int iFace = 1; 1414 static G4ThreadLocal G4int iFace = 1;
1415 static G4ThreadLocal G4int iQVertex = 0; 1415 static G4ThreadLocal G4int iQVertex = 0;
1416 G4int vIndex = pF[iFace].edge[iQVertex].v; 1416 G4int vIndex = pF[iFace].edge[iQVertex].v;
1417 1417
1418 edgeFlag = (vIndex > 0) ? 1 : 0; 1418 edgeFlag = (vIndex > 0) ? 1 : 0;
1419 index = std::abs(vIndex); 1419 index = std::abs(vIndex);
1420 1420
1421 if (iQVertex >= 3 || pF[iFace].edge[iQVerte 1421 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1422 iQVertex = 0; 1422 iQVertex = 0;
1423 if (++iFace > nface) iFace = 1; 1423 if (++iFace > nface) iFace = 1;
1424 return false; // Last Edge 1424 return false; // Last Edge
1425 } 1425 }
1426 1426
1427 ++iQVertex; 1427 ++iQVertex;
1428 return true; // not Last Edge 1428 return true; // not Last Edge
1429 } 1429 }
1430 1430
1431 G4Point3D HepPolyhedron::GetVertex(G4int inde 1431 G4Point3D HepPolyhedron::GetVertex(G4int index) const
1432 /******************************************** 1432 /***********************************************************************
1433 * 1433 * *
1434 * Name: HepPolyhedron::GetVertex 1434 * Name: HepPolyhedron::GetVertex Date: 03.09.96 *
1435 * Author: Yasuhide Sawada 1435 * Author: Yasuhide Sawada Revised: 17.11.99 *
1436 * 1436 * *
1437 * Function: Get vertex of the index. 1437 * Function: Get vertex of the index. *
1438 * 1438 * *
1439 ******************************************** 1439 ***********************************************************************/
1440 { 1440 {
1441 if (index <= 0 || index > nvert) { 1441 if (index <= 0 || index > nvert) {
1442 std::cerr 1442 std::cerr
1443 << "HepPolyhedron::GetVertex: irrelevan 1443 << "HepPolyhedron::GetVertex: irrelevant index " << index
1444 << std::endl; 1444 << std::endl;
1445 return G4Point3D(); 1445 return G4Point3D();
1446 } 1446 }
1447 return pV[index]; 1447 return pV[index];
1448 } 1448 }
1449 1449
1450 G4bool 1450 G4bool
1451 HepPolyhedron::GetNextVertex(G4Point3D &verte 1451 HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag) const
1452 /******************************************** 1452 /***********************************************************************
1453 * 1453 * *
1454 * Name: HepPolyhedron::GetNextVertex 1454 * Name: HepPolyhedron::GetNextVertex Date: 22.07.96 *
1455 * Author: John Allison 1455 * Author: John Allison Revised: *
1456 * 1456 * *
1457 * Function: Get vertices of the quadrilatera 1457 * Function: Get vertices of the quadrilaterals in order for each *
1458 * face in face order. Returns fal 1458 * face in face order. Returns false when finished each *
1459 * face. 1459 * face. *
1460 * 1460 * *
1461 ******************************************** 1461 ***********************************************************************/
1462 { 1462 {
1463 G4int index; 1463 G4int index;
1464 G4bool rep = GetNextVertexIndex(index, edge 1464 G4bool rep = GetNextVertexIndex(index, edgeFlag);
1465 vertex = pV[index]; 1465 vertex = pV[index];
1466 return rep; 1466 return rep;
1467 } 1467 }
1468 1468
1469 G4bool HepPolyhedron::GetNextVertex(G4Point3D 1469 G4bool HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag,
1470 G4Normal3D 1470 G4Normal3D &normal) const
1471 /******************************************** 1471 /***********************************************************************
1472 * 1472 * *
1473 * Name: HepPolyhedron::GetNextVertex 1473 * Name: HepPolyhedron::GetNextVertex Date: 26.11.99 *
1474 * Author: E.Chernyaev 1474 * Author: E.Chernyaev Revised: *
1475 * 1475 * *
1476 * Function: Get vertices with normals of the 1476 * Function: Get vertices with normals of the quadrilaterals in order *
1477 * for each face in face order. 1477 * for each face in face order. *
1478 * Returns false when finished each 1478 * Returns false when finished each face. *
1479 * 1479 * *
1480 ******************************************** 1480 ***********************************************************************/
1481 { 1481 {
1482 static G4ThreadLocal G4int iFace = 1; 1482 static G4ThreadLocal G4int iFace = 1;
1483 static G4ThreadLocal G4int iNode = 0; 1483 static G4ThreadLocal G4int iNode = 0;
1484 1484
1485 if (nface == 0) return false; // empty pol 1485 if (nface == 0) return false; // empty polyhedron
1486 1486
1487 G4int k = pF[iFace].edge[iNode].v; 1487 G4int k = pF[iFace].edge[iNode].v;
1488 if (k > 0) { edgeFlag = 1; } else { edgeFla 1488 if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
1489 vertex = pV[k]; 1489 vertex = pV[k];
1490 normal = FindNodeNormal(iFace,k); 1490 normal = FindNodeNormal(iFace,k);
1491 if (iNode >= 3 || pF[iFace].edge[iNode+1].v 1491 if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
1492 iNode = 0; 1492 iNode = 0;
1493 if (++iFace > nface) iFace = 1; 1493 if (++iFace > nface) iFace = 1;
1494 return false; // last node 1494 return false; // last node
1495 } 1495 }
1496 ++iNode; 1496 ++iNode;
1497 return true; // not last no 1497 return true; // not last node
1498 } 1498 }
1499 1499
1500 G4bool HepPolyhedron::GetNextEdgeIndices(G4in 1500 G4bool HepPolyhedron::GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag,
1501 G4int 1501 G4int &iface1, G4int &iface2) const
1502 /******************************************** 1502 /***********************************************************************
1503 * 1503 * *
1504 * Name: HepPolyhedron::GetNextEdgeIndices 1504 * Name: HepPolyhedron::GetNextEdgeIndices Date: 30.09.96 *
1505 * Author: E.Chernyaev 1505 * Author: E.Chernyaev Revised: 17.11.99 *
1506 * 1506 * *
1507 * Function: Get indices of the next edge tog 1507 * Function: Get indices of the next edge together with indices of *
1508 * of the faces which share the edg 1508 * of the faces which share the edge. *
1509 * Returns false when the last edge 1509 * Returns false when the last edge. *
1510 * 1510 * *
1511 ******************************************** 1511 ***********************************************************************/
1512 { 1512 {
1513 static G4ThreadLocal G4int iFace = 1; 1513 static G4ThreadLocal G4int iFace = 1;
1514 static G4ThreadLocal G4int iQVertex = 0; 1514 static G4ThreadLocal G4int iQVertex = 0;
1515 static G4ThreadLocal G4int iOrder = 1; 1515 static G4ThreadLocal G4int iOrder = 1;
1516 G4int k1, k2, kflag, kface1, kface2; 1516 G4int k1, k2, kflag, kface1, kface2;
1517 1517
1518 if (iFace == 1 && iQVertex == 0) { 1518 if (iFace == 1 && iQVertex == 0) {
1519 k2 = pF[nface].edge[0].v; 1519 k2 = pF[nface].edge[0].v;
1520 k1 = pF[nface].edge[3].v; 1520 k1 = pF[nface].edge[3].v;
1521 if (k1 == 0) k1 = pF[nface].edge[2].v; 1521 if (k1 == 0) k1 = pF[nface].edge[2].v;
1522 if (std::abs(k1) > std::abs(k2)) iOrder = 1522 if (std::abs(k1) > std::abs(k2)) iOrder = -1;
1523 } 1523 }
1524 1524
1525 do { 1525 do {
1526 k1 = pF[iFace].edge[iQVertex].v; 1526 k1 = pF[iFace].edge[iQVertex].v;
1527 kflag = k1; 1527 kflag = k1;
1528 k1 = std::abs(k1); 1528 k1 = std::abs(k1);
1529 kface1 = iFace; 1529 kface1 = iFace;
1530 kface2 = pF[iFace].edge[iQVertex].f; 1530 kface2 = pF[iFace].edge[iQVertex].f;
1531 if (iQVertex >= 3 || pF[iFace].edge[iQVer 1531 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1532 iQVertex = 0; 1532 iQVertex = 0;
1533 k2 = std::abs(pF[iFace].edge[iQVertex]. 1533 k2 = std::abs(pF[iFace].edge[iQVertex].v);
1534 iFace++; 1534 iFace++;
1535 }else{ 1535 }else{
1536 iQVertex++; 1536 iQVertex++;
1537 k2 = std::abs(pF[iFace].edge[iQVertex]. 1537 k2 = std::abs(pF[iFace].edge[iQVertex].v);
1538 } 1538 }
1539 } while (iOrder*k1 > iOrder*k2); 1539 } while (iOrder*k1 > iOrder*k2);
1540 1540
1541 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1541 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1542 iface1 = kface1; iface2 = kface2; 1542 iface1 = kface1; iface2 = kface2;
1543 1543
1544 if (iFace > nface) { 1544 if (iFace > nface) {
1545 iFace = 1; iOrder = 1; 1545 iFace = 1; iOrder = 1;
1546 return false; 1546 return false;
1547 } 1547 }
1548 1548
1549 return true; 1549 return true;
1550 } 1550 }
1551 1551
1552 G4bool 1552 G4bool
1553 HepPolyhedron::GetNextEdgeIndices(G4int &i1, 1553 HepPolyhedron::GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag) const
1554 /******************************************** 1554 /***********************************************************************
1555 * 1555 * *
1556 * Name: HepPolyhedron::GetNextEdgeIndices 1556 * Name: HepPolyhedron::GetNextEdgeIndices Date: 17.11.99 *
1557 * Author: E.Chernyaev 1557 * Author: E.Chernyaev Revised: *
1558 * 1558 * *
1559 * Function: Get indices of the next edge. 1559 * Function: Get indices of the next edge. *
1560 * Returns false when the last edge 1560 * Returns false when the last edge. *
1561 * 1561 * *
1562 ******************************************** 1562 ***********************************************************************/
1563 { 1563 {
1564 G4int kface1, kface2; 1564 G4int kface1, kface2;
1565 return GetNextEdgeIndices(i1, i2, edgeFlag, 1565 return GetNextEdgeIndices(i1, i2, edgeFlag, kface1, kface2);
1566 } 1566 }
1567 1567
1568 G4bool 1568 G4bool
1569 HepPolyhedron::GetNextEdge(G4Point3D &p1, 1569 HepPolyhedron::GetNextEdge(G4Point3D &p1,
1570 G4Point3D &p2, 1570 G4Point3D &p2,
1571 G4int &edgeFlag) c 1571 G4int &edgeFlag) const
1572 /******************************************** 1572 /***********************************************************************
1573 * 1573 * *
1574 * Name: HepPolyhedron::GetNextEdge 1574 * Name: HepPolyhedron::GetNextEdge Date: 30.09.96 *
1575 * Author: E.Chernyaev 1575 * Author: E.Chernyaev Revised: *
1576 * 1576 * *
1577 * Function: Get next edge. 1577 * Function: Get next edge. *
1578 * Returns false when the last edge 1578 * Returns false when the last edge. *
1579 * 1579 * *
1580 ******************************************** 1580 ***********************************************************************/
1581 { 1581 {
1582 G4int i1,i2; 1582 G4int i1,i2;
1583 G4bool rep = GetNextEdgeIndices(i1,i2,edgeF 1583 G4bool rep = GetNextEdgeIndices(i1,i2,edgeFlag);
1584 p1 = pV[i1]; 1584 p1 = pV[i1];
1585 p2 = pV[i2]; 1585 p2 = pV[i2];
1586 return rep; 1586 return rep;
1587 } 1587 }
1588 1588
1589 G4bool 1589 G4bool
1590 HepPolyhedron::GetNextEdge(G4Point3D &p1, G4P 1590 HepPolyhedron::GetNextEdge(G4Point3D &p1, G4Point3D &p2,
1591 G4int &edgeFlag, G4 1591 G4int &edgeFlag, G4int &iface1, G4int &iface2) const
1592 /******************************************** 1592 /***********************************************************************
1593 * 1593 * *
1594 * Name: HepPolyhedron::GetNextEdge 1594 * Name: HepPolyhedron::GetNextEdge Date: 17.11.99 *
1595 * Author: E.Chernyaev 1595 * Author: E.Chernyaev Revised: *
1596 * 1596 * *
1597 * Function: Get next edge with indices of th 1597 * Function: Get next edge with indices of the faces which share *
1598 * the edge. 1598 * the edge. *
1599 * Returns false when the last edge 1599 * Returns false when the last edge. *
1600 * 1600 * *
1601 ******************************************** 1601 ***********************************************************************/
1602 { 1602 {
1603 G4int i1,i2; 1603 G4int i1,i2;
1604 G4bool rep = GetNextEdgeIndices(i1,i2,edgeF 1604 G4bool rep = GetNextEdgeIndices(i1,i2,edgeFlag,iface1,iface2);
1605 p1 = pV[i1]; 1605 p1 = pV[i1];
1606 p2 = pV[i2]; 1606 p2 = pV[i2];
1607 return rep; 1607 return rep;
1608 } 1608 }
1609 1609
1610 void HepPolyhedron::GetFacet(G4int iFace, G4i 1610 void HepPolyhedron::GetFacet(G4int iFace, G4int &n, G4int *iNodes,
1611 G4int *edgeFlags, 1611 G4int *edgeFlags, G4int *iFaces) const
1612 /******************************************** 1612 /***********************************************************************
1613 * 1613 * *
1614 * Name: HepPolyhedron::GetFacet 1614 * Name: HepPolyhedron::GetFacet Date: 15.12.99 *
1615 * Author: E.Chernyaev 1615 * Author: E.Chernyaev Revised: *
1616 * 1616 * *
1617 * Function: Get face by index 1617 * Function: Get face by index *
1618 * 1618 * *
1619 ******************************************** 1619 ***********************************************************************/
1620 { 1620 {
1621 if (iFace < 1 || iFace > nface) { 1621 if (iFace < 1 || iFace > nface) {
1622 std::cerr 1622 std::cerr
1623 << "HepPolyhedron::GetFacet: irrelevant 1623 << "HepPolyhedron::GetFacet: irrelevant index " << iFace
1624 << std::endl; 1624 << std::endl;
1625 n = 0; 1625 n = 0;
1626 }else{ 1626 }else{
1627 G4int i, k; 1627 G4int i, k;
1628 for (i=0; i<4; i++) { 1628 for (i=0; i<4; i++) {
1629 k = pF[iFace].edge[i].v; 1629 k = pF[iFace].edge[i].v;
1630 if (k == 0) break; 1630 if (k == 0) break;
1631 if (iFaces != nullptr) iFaces[i] = pF[i 1631 if (iFaces != nullptr) iFaces[i] = pF[iFace].edge[i].f;
1632 if (k > 0) { 1632 if (k > 0) {
1633 iNodes[i] = k; 1633 iNodes[i] = k;
1634 if (edgeFlags != nullptr) edgeFlags[i 1634 if (edgeFlags != nullptr) edgeFlags[i] = 1;
1635 }else{ 1635 }else{
1636 iNodes[i] = -k; 1636 iNodes[i] = -k;
1637 if (edgeFlags != nullptr) edgeFlags[i 1637 if (edgeFlags != nullptr) edgeFlags[i] = -1;
1638 } 1638 }
1639 } 1639 }
1640 n = i; 1640 n = i;
1641 } 1641 }
1642 } 1642 }
1643 1643
1644 void HepPolyhedron::GetFacet(G4int index, G4i 1644 void HepPolyhedron::GetFacet(G4int index, G4int &n, G4Point3D *nodes,
1645 G4int *edgeFlags 1645 G4int *edgeFlags, G4Normal3D *normals) const
1646 /******************************************** 1646 /***********************************************************************
1647 * 1647 * *
1648 * Name: HepPolyhedron::GetFacet 1648 * Name: HepPolyhedron::GetFacet Date: 17.11.99 *
1649 * Author: E.Chernyaev 1649 * Author: E.Chernyaev Revised: *
1650 * 1650 * *
1651 * Function: Get face by index 1651 * Function: Get face by index *
1652 * 1652 * *
1653 ******************************************** 1653 ***********************************************************************/
1654 { 1654 {
1655 G4int iNodes[4]; 1655 G4int iNodes[4];
1656 GetFacet(index, n, iNodes, edgeFlags); 1656 GetFacet(index, n, iNodes, edgeFlags);
1657 if (n != 0) { 1657 if (n != 0) {
1658 for (G4int i=0; i<n; i++) { 1658 for (G4int i=0; i<n; i++) {
1659 nodes[i] = pV[iNodes[i]]; 1659 nodes[i] = pV[iNodes[i]];
1660 if (normals != nullptr) normals[i] = Fi 1660 if (normals != nullptr) normals[i] = FindNodeNormal(index,iNodes[i]);
1661 } 1661 }
1662 } 1662 }
1663 } 1663 }
1664 1664
1665 G4bool 1665 G4bool
1666 HepPolyhedron::GetNextFacet(G4int &n, G4Point 1666 HepPolyhedron::GetNextFacet(G4int &n, G4Point3D *nodes,
1667 G4int *edgeFlags, 1667 G4int *edgeFlags, G4Normal3D *normals) const
1668 /******************************************** 1668 /***********************************************************************
1669 * 1669 * *
1670 * Name: HepPolyhedron::GetNextFacet 1670 * Name: HepPolyhedron::GetNextFacet Date: 19.11.99 *
1671 * Author: E.Chernyaev 1671 * Author: E.Chernyaev Revised: *
1672 * 1672 * *
1673 * Function: Get next face with normals of un 1673 * Function: Get next face with normals of unit length at the nodes. *
1674 * Returns false when finished all 1674 * Returns false when finished all faces. *
1675 * 1675 * *
1676 ******************************************** 1676 ***********************************************************************/
1677 { 1677 {
1678 static G4ThreadLocal G4int iFace = 1; 1678 static G4ThreadLocal G4int iFace = 1;
1679 1679
1680 if (edgeFlags == nullptr) { 1680 if (edgeFlags == nullptr) {
1681 GetFacet(iFace, n, nodes); 1681 GetFacet(iFace, n, nodes);
1682 }else if (normals == nullptr) { 1682 }else if (normals == nullptr) {
1683 GetFacet(iFace, n, nodes, edgeFlags); 1683 GetFacet(iFace, n, nodes, edgeFlags);
1684 }else{ 1684 }else{
1685 GetFacet(iFace, n, nodes, edgeFlags, norm 1685 GetFacet(iFace, n, nodes, edgeFlags, normals);
1686 } 1686 }
1687 1687
1688 if (++iFace > nface) { 1688 if (++iFace > nface) {
1689 iFace = 1; 1689 iFace = 1;
1690 return false; 1690 return false;
1691 } 1691 }
1692 1692
1693 return true; 1693 return true;
1694 } 1694 }
1695 1695
1696 G4Normal3D HepPolyhedron::GetNormal(G4int iFa 1696 G4Normal3D HepPolyhedron::GetNormal(G4int iFace) const
1697 /******************************************** 1697 /***********************************************************************
1698 * 1698 * *
1699 * Name: HepPolyhedron::GetNormal 1699 * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1700 * Author: E.Chernyaev 1700 * Author: E.Chernyaev Revised: *
1701 * 1701 * *
1702 * Function: Get normal of the face given by 1702 * Function: Get normal of the face given by index *
1703 * 1703 * *
1704 ******************************************** 1704 ***********************************************************************/
1705 { 1705 {
1706 if (iFace < 1 || iFace > nface) { 1706 if (iFace < 1 || iFace > nface) {
1707 std::cerr 1707 std::cerr
1708 << "HepPolyhedron::GetNormal: irrelevan 1708 << "HepPolyhedron::GetNormal: irrelevant index " << iFace
1709 << std::endl; 1709 << std::endl;
1710 return G4Normal3D(); 1710 return G4Normal3D();
1711 } 1711 }
1712 1712
1713 G4int i0 = std::abs(pF[iFace].edge[0].v); 1713 G4int i0 = std::abs(pF[iFace].edge[0].v);
1714 G4int i1 = std::abs(pF[iFace].edge[1].v); 1714 G4int i1 = std::abs(pF[iFace].edge[1].v);
1715 G4int i2 = std::abs(pF[iFace].edge[2].v); 1715 G4int i2 = std::abs(pF[iFace].edge[2].v);
1716 G4int i3 = std::abs(pF[iFace].edge[3].v); 1716 G4int i3 = std::abs(pF[iFace].edge[3].v);
1717 if (i3 == 0) i3 = i0; 1717 if (i3 == 0) i3 = i0;
1718 return (pV[i2] - pV[i0]).cross(pV[i3] - pV[ 1718 return (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1]);
1719 } 1719 }
1720 1720
1721 G4Normal3D HepPolyhedron::GetUnitNormal(G4int 1721 G4Normal3D HepPolyhedron::GetUnitNormal(G4int iFace) const
1722 /******************************************** 1722 /***********************************************************************
1723 * 1723 * *
1724 * Name: HepPolyhedron::GetNormal 1724 * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1725 * Author: E.Chernyaev 1725 * Author: E.Chernyaev Revised: *
1726 * 1726 * *
1727 * Function: Get unit normal of the face give 1727 * Function: Get unit normal of the face given by index *
1728 * 1728 * *
1729 ******************************************** 1729 ***********************************************************************/
1730 { 1730 {
1731 if (iFace < 1 || iFace > nface) { 1731 if (iFace < 1 || iFace > nface) {
1732 std::cerr 1732 std::cerr
1733 << "HepPolyhedron::GetUnitNormal: irrel 1733 << "HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
1734 << std::endl; 1734 << std::endl;
1735 return G4Normal3D(); 1735 return G4Normal3D();
1736 } 1736 }
1737 1737
1738 G4int i0 = std::abs(pF[iFace].edge[0].v); 1738 G4int i0 = std::abs(pF[iFace].edge[0].v);
1739 G4int i1 = std::abs(pF[iFace].edge[1].v); 1739 G4int i1 = std::abs(pF[iFace].edge[1].v);
1740 G4int i2 = std::abs(pF[iFace].edge[2].v); 1740 G4int i2 = std::abs(pF[iFace].edge[2].v);
1741 G4int i3 = std::abs(pF[iFace].edge[3].v); 1741 G4int i3 = std::abs(pF[iFace].edge[3].v);
1742 if (i3 == 0) i3 = i0; 1742 if (i3 == 0) i3 = i0;
1743 return ((pV[i2] - pV[i0]).cross(pV[i3] - pV 1743 return ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).unit();
1744 } 1744 }
1745 1745
1746 G4bool HepPolyhedron::GetNextNormal(G4Normal3 1746 G4bool HepPolyhedron::GetNextNormal(G4Normal3D &normal) const
1747 /******************************************** 1747 /***********************************************************************
1748 * 1748 * *
1749 * Name: HepPolyhedron::GetNextNormal 1749 * Name: HepPolyhedron::GetNextNormal Date: 22.07.96 *
1750 * Author: John Allison 1750 * Author: John Allison Revised: 19.11.99 *
1751 * 1751 * *
1752 * Function: Get normals of each face in face 1752 * Function: Get normals of each face in face order. Returns false *
1753 * when finished all faces. 1753 * when finished all faces. *
1754 * 1754 * *
1755 ******************************************** 1755 ***********************************************************************/
1756 { 1756 {
1757 static G4ThreadLocal G4int iFace = 1; 1757 static G4ThreadLocal G4int iFace = 1;
1758 normal = GetNormal(iFace); 1758 normal = GetNormal(iFace);
1759 if (++iFace > nface) { 1759 if (++iFace > nface) {
1760 iFace = 1; 1760 iFace = 1;
1761 return false; 1761 return false;
1762 } 1762 }
1763 return true; 1763 return true;
1764 } 1764 }
1765 1765
1766 G4bool HepPolyhedron::GetNextUnitNormal(G4Nor 1766 G4bool HepPolyhedron::GetNextUnitNormal(G4Normal3D &normal) const
1767 /******************************************** 1767 /***********************************************************************
1768 * 1768 * *
1769 * Name: HepPolyhedron::GetNextUnitNormal 1769 * Name: HepPolyhedron::GetNextUnitNormal Date: 16.09.96 *
1770 * Author: E.Chernyaev 1770 * Author: E.Chernyaev Revised: *
1771 * 1771 * *
1772 * Function: Get normals of unit length of ea 1772 * Function: Get normals of unit length of each face in face order. *
1773 * Returns false when finished all 1773 * Returns false when finished all faces. *
1774 * 1774 * *
1775 ******************************************** 1775 ***********************************************************************/
1776 { 1776 {
1777 G4bool rep = GetNextNormal(normal); 1777 G4bool rep = GetNextNormal(normal);
1778 normal = normal.unit(); 1778 normal = normal.unit();
1779 return rep; 1779 return rep;
1780 } 1780 }
1781 1781
1782 G4double HepPolyhedron::GetSurfaceArea() cons 1782 G4double HepPolyhedron::GetSurfaceArea() const
1783 /******************************************** 1783 /***********************************************************************
1784 * 1784 * *
1785 * Name: HepPolyhedron::GetSurfaceArea 1785 * Name: HepPolyhedron::GetSurfaceArea Date: 25.05.01 *
1786 * Author: E.Chernyaev 1786 * Author: E.Chernyaev Revised: *
1787 * 1787 * *
1788 * Function: Returns area of the surface of t 1788 * Function: Returns area of the surface of the polyhedron. *
1789 * 1789 * *
1790 ******************************************** 1790 ***********************************************************************/
1791 { 1791 {
1792 G4double srf = 0.; 1792 G4double srf = 0.;
1793 for (G4int iFace=1; iFace<=nface; iFace++) 1793 for (G4int iFace=1; iFace<=nface; iFace++) {
1794 G4int i0 = std::abs(pF[iFace].edge[0].v); 1794 G4int i0 = std::abs(pF[iFace].edge[0].v);
1795 G4int i1 = std::abs(pF[iFace].edge[1].v); 1795 G4int i1 = std::abs(pF[iFace].edge[1].v);
1796 G4int i2 = std::abs(pF[iFace].edge[2].v); 1796 G4int i2 = std::abs(pF[iFace].edge[2].v);
1797 G4int i3 = std::abs(pF[iFace].edge[3].v); 1797 G4int i3 = std::abs(pF[iFace].edge[3].v);
1798 if (i3 == 0) i3 = i0; 1798 if (i3 == 0) i3 = i0;
1799 srf += ((pV[i2] - pV[i0]).cross(pV[i3] - 1799 srf += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).mag();
1800 } 1800 }
1801 return srf/2.; 1801 return srf/2.;
1802 } 1802 }
1803 1803
1804 G4double HepPolyhedron::GetVolume() const 1804 G4double HepPolyhedron::GetVolume() const
1805 /******************************************** 1805 /***********************************************************************
1806 * 1806 * *
1807 * Name: HepPolyhedron::GetVolume 1807 * Name: HepPolyhedron::GetVolume Date: 25.05.01 *
1808 * Author: E.Chernyaev 1808 * Author: E.Chernyaev Revised: *
1809 * 1809 * *
1810 * Function: Returns volume of the polyhedron 1810 * Function: Returns volume of the polyhedron. *
1811 * 1811 * *
1812 ******************************************** 1812 ***********************************************************************/
1813 { 1813 {
1814 G4double v = 0.; 1814 G4double v = 0.;
1815 for (G4int iFace=1; iFace<=nface; iFace++) 1815 for (G4int iFace=1; iFace<=nface; iFace++) {
1816 G4int i0 = std::abs(pF[iFace].edge[0].v); 1816 G4int i0 = std::abs(pF[iFace].edge[0].v);
1817 G4int i1 = std::abs(pF[iFace].edge[1].v); 1817 G4int i1 = std::abs(pF[iFace].edge[1].v);
1818 G4int i2 = std::abs(pF[iFace].edge[2].v); 1818 G4int i2 = std::abs(pF[iFace].edge[2].v);
1819 G4int i3 = std::abs(pF[iFace].edge[3].v); 1819 G4int i3 = std::abs(pF[iFace].edge[3].v);
1820 G4Point3D pt; 1820 G4Point3D pt;
1821 if (i3 == 0) { 1821 if (i3 == 0) {
1822 i3 = i0; 1822 i3 = i0;
1823 pt = (pV[i0]+pV[i1]+pV[i2]) * (1./3.); 1823 pt = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);
1824 }else{ 1824 }else{
1825 pt = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0. 1825 pt = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0.25;
1826 } 1826 }
1827 v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV 1827 v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).dot(pt);
1828 } 1828 }
1829 return v/6.; 1829 return v/6.;
1830 } 1830 }
1831 1831
1832 G4int 1832 G4int
1833 HepPolyhedron::createTwistedTrap(G4double Dz, 1833 HepPolyhedron::createTwistedTrap(G4double Dz,
1834 const G4doub 1834 const G4double xy1[][2],
1835 const G4doub 1835 const G4double xy2[][2])
1836 /******************************************** 1836 /***********************************************************************
1837 * 1837 * *
1838 * Name: createTwistedTrap 1838 * Name: createTwistedTrap Date: 05.11.02 *
1839 * Author: E.Chernyaev (IHEP/Protvino) 1839 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1840 * 1840 * *
1841 * Function: Creates polyhedron for twisted t 1841 * Function: Creates polyhedron for twisted trapezoid *
1842 * 1842 * *
1843 * Input: Dz - half-length along Z 1843 * Input: Dz - half-length along Z 8----7 *
1844 * xy1[2,4] - quadrilateral at Z=-Dz 1844 * xy1[2,4] - quadrilateral at Z=-Dz 5----6 ! *
1845 * xy2[2,4] - quadrilateral at Z=+Dz 1845 * xy2[2,4] - quadrilateral at Z=+Dz ! 4-!--3 *
1846 * 1846 * 1----2 *
1847 * 1847 * *
1848 ******************************************** 1848 ***********************************************************************/
1849 { 1849 {
1850 AllocateMemory(12,18); 1850 AllocateMemory(12,18);
1851 1851
1852 pV[ 1] = G4Point3D(xy1[0][0],xy1[0][1],-Dz) 1852 pV[ 1] = G4Point3D(xy1[0][0],xy1[0][1],-Dz);
1853 pV[ 2] = G4Point3D(xy1[1][0],xy1[1][1],-Dz) 1853 pV[ 2] = G4Point3D(xy1[1][0],xy1[1][1],-Dz);
1854 pV[ 3] = G4Point3D(xy1[2][0],xy1[2][1],-Dz) 1854 pV[ 3] = G4Point3D(xy1[2][0],xy1[2][1],-Dz);
1855 pV[ 4] = G4Point3D(xy1[3][0],xy1[3][1],-Dz) 1855 pV[ 4] = G4Point3D(xy1[3][0],xy1[3][1],-Dz);
1856 1856
1857 pV[ 5] = G4Point3D(xy2[0][0],xy2[0][1], Dz) 1857 pV[ 5] = G4Point3D(xy2[0][0],xy2[0][1], Dz);
1858 pV[ 6] = G4Point3D(xy2[1][0],xy2[1][1], Dz) 1858 pV[ 6] = G4Point3D(xy2[1][0],xy2[1][1], Dz);
1859 pV[ 7] = G4Point3D(xy2[2][0],xy2[2][1], Dz) 1859 pV[ 7] = G4Point3D(xy2[2][0],xy2[2][1], Dz);
1860 pV[ 8] = G4Point3D(xy2[3][0],xy2[3][1], Dz) 1860 pV[ 8] = G4Point3D(xy2[3][0],xy2[3][1], Dz);
1861 1861
1862 pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.; 1862 pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;
1863 pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.; 1863 pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;
1864 pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.; 1864 pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;
1865 pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.; 1865 pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;
1866 1866
1867 enum {DUMMY, BOTTOM, 1867 enum {DUMMY, BOTTOM,
1868 LEFT_BOTTOM, LEFT_FRONT, LEFT_TOP, 1868 LEFT_BOTTOM, LEFT_FRONT, LEFT_TOP, LEFT_BACK,
1869 BACK_BOTTOM, BACK_LEFT, BACK_TOP, 1869 BACK_BOTTOM, BACK_LEFT, BACK_TOP, BACK_RIGHT,
1870 RIGHT_BOTTOM, RIGHT_BACK, RIGHT_TOP 1870 RIGHT_BOTTOM, RIGHT_BACK, RIGHT_TOP, RIGHT_FRONT,
1871 FRONT_BOTTOM, FRONT_RIGHT, FRONT_TOP 1871 FRONT_BOTTOM, FRONT_RIGHT, FRONT_TOP, FRONT_LEFT,
1872 TOP}; 1872 TOP};
1873 1873
1874 pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM 1874 pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
1875 1875
1876 pF[ 2]=G4Facet(4,BOTTOM, -1,LEFT_FRONT, 1876 pF[ 2]=G4Facet(4,BOTTOM, -1,LEFT_FRONT, -12,LEFT_BACK, 0,0);
1877 pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP, 1877 pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP, -12,LEFT_BOTTOM, 0,0);
1878 pF[ 4]=G4Facet(5,TOP, -8,LEFT_BACK, 1878 pF[ 4]=G4Facet(5,TOP, -8,LEFT_BACK, -12,LEFT_FRONT, 0,0);
1879 pF[ 5]=G4Facet(8,BACK_LEFT, -4,LEFT_BOTTOM 1879 pF[ 5]=G4Facet(8,BACK_LEFT, -4,LEFT_BOTTOM, -12,LEFT_TOP, 0,0);
1880 1880
1881 pF[ 6]=G4Facet(3,BOTTOM, -4,BACK_LEFT, 1881 pF[ 6]=G4Facet(3,BOTTOM, -4,BACK_LEFT, -11,BACK_RIGHT, 0,0);
1882 pF[ 7]=G4Facet(4,LEFT_BACK, -8,BACK_TOP, 1882 pF[ 7]=G4Facet(4,LEFT_BACK, -8,BACK_TOP, -11,BACK_BOTTOM, 0,0);
1883 pF[ 8]=G4Facet(8,TOP, -7,BACK_RIGHT, 1883 pF[ 8]=G4Facet(8,TOP, -7,BACK_RIGHT, -11,BACK_LEFT, 0,0);
1884 pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM 1884 pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP, 0,0);
1885 1885
1886 pF[10]=G4Facet(2,BOTTOM, -3,RIGHT_BACK, 1886 pF[10]=G4Facet(2,BOTTOM, -3,RIGHT_BACK, -10,RIGHT_FRONT, 0,0);
1887 pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP, 1887 pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP, -10,RIGHT_BOTTOM, 0,0);
1888 pF[12]=G4Facet(7,TOP, -6,RIGHT_FRONT 1888 pF[12]=G4Facet(7,TOP, -6,RIGHT_FRONT, -10,RIGHT_BACK, 0,0);
1889 pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTO 1889 pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP, 0,0);
1890 1890
1891 pF[14]=G4Facet(1,BOTTOM, -2,FRONT_RIGHT 1891 pF[14]=G4Facet(1,BOTTOM, -2,FRONT_RIGHT, -9,FRONT_LEFT, 0,0);
1892 pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP, 1892 pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP, -9,FRONT_BOTTOM, 0,0);
1893 pF[16]=G4Facet(6,TOP, -5,FRONT_LEFT, 1893 pF[16]=G4Facet(6,TOP, -5,FRONT_LEFT, -9,FRONT_RIGHT, 0,0);
1894 pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTO 1894 pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP, 0,0);
1895 1895
1896 pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7, 1896 pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
1897 1897
1898 return 0; 1898 return 0;
1899 } 1899 }
1900 1900
1901 G4int 1901 G4int
1902 HepPolyhedron::createPolyhedron(G4int Nnodes, 1902 HepPolyhedron::createPolyhedron(G4int Nnodes, G4int Nfaces,
1903 const G4doubl 1903 const G4double xyz[][3],
1904 const G4int 1904 const G4int faces[][4])
1905 /******************************************** 1905 /***********************************************************************
1906 * 1906 * *
1907 * Name: createPolyhedron 1907 * Name: createPolyhedron Date: 05.11.02 *
1908 * Author: E.Chernyaev (IHEP/Protvino) 1908 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1909 * 1909 * *
1910 * Function: Creates user defined polyhedron 1910 * Function: Creates user defined polyhedron *
1911 * 1911 * *
1912 * Input: Nnodes - number of nodes 1912 * Input: Nnodes - number of nodes *
1913 * Nfaces - number of faces 1913 * Nfaces - number of faces *
1914 * nodes[][3] - node coordinates 1914 * nodes[][3] - node coordinates *
1915 * faces[][4] - faces 1915 * faces[][4] - faces *
1916 * 1916 * *
1917 ******************************************** 1917 ***********************************************************************/
1918 { 1918 {
1919 AllocateMemory(Nnodes, Nfaces); 1919 AllocateMemory(Nnodes, Nfaces);
1920 if (nvert == 0) return 1; 1920 if (nvert == 0) return 1;
1921 1921
1922 for (G4int i=0; i<Nnodes; i++) { 1922 for (G4int i=0; i<Nnodes; i++) {
1923 pV[i+1] = G4Point3D(xyz[i][0], xyz[i][1], 1923 pV[i+1] = G4Point3D(xyz[i][0], xyz[i][1], xyz[i][2]);
1924 } 1924 }
1925 for (G4int k=0; k<Nfaces; k++) { 1925 for (G4int k=0; k<Nfaces; k++) {
1926 pF[k+1] = G4Facet(faces[k][0],0,faces[k][ 1926 pF[k+1] = G4Facet(faces[k][0],0,faces[k][1],0,faces[k][2],0,faces[k][3],0);
1927 } 1927 }
1928 SetReferences(); 1928 SetReferences();
1929 return 0; 1929 return 0;
1930 } 1930 }
1931 1931
1932 G4Point3D HepPolyhedron::vertexUnweightedMean 1932 G4Point3D HepPolyhedron::vertexUnweightedMean() const {
1933 /****************************************** 1933 /***********************************************************************
1934 * 1934 * *
1935 * Name: vertexUnweightedMean 1935 * Name: vertexUnweightedMean Date: 23.10.23 *
1936 * Author: S. Boogert (Manchester) 1936 * Author: S. Boogert (Manchester) Revised: *
1937 * 1937 * *
1938 * Function: Calculate the unweighted mean 1938 * Function: Calculate the unweighted mean of all the vertices *
1939 * in the polyhedron. Not to be confused wi 1939 * in the polyhedron. Not to be confused with the polyhedron centre or *
1940 * centre of mass 1940 * centre of mass *
1941 ****************************************** 1941 ***********************************************************************/
1942 1942
1943 auto centre = G4Point3D(); 1943 auto centre = G4Point3D();
1944 for(int i=1;i<=nvert;i++) { 1944 for(int i=1;i<=nvert;i++) {
1945 centre += pV[i]; 1945 centre += pV[i];
1946 } 1946 }
1947 centre = centre/nvert; 1947 centre = centre/nvert;
1948 return centre; 1948 return centre;
1949 } 1949 }
1950 1950
1951 HepPolyhedronTrd2::HepPolyhedronTrd2(G4double 1951 HepPolyhedronTrd2::HepPolyhedronTrd2(G4double Dx1, G4double Dx2,
1952 G4double 1952 G4double Dy1, G4double Dy2,
1953 G4double 1953 G4double Dz)
1954 /******************************************** 1954 /***********************************************************************
1955 * 1955 * *
1956 * Name: HepPolyhedronTrd2 1956 * Name: HepPolyhedronTrd2 Date: 22.07.96 *
1957 * Author: E.Chernyaev (IHEP/Protvino) 1957 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1958 * 1958 * *
1959 * Function: Create GEANT4 TRD2-trapezoid 1959 * Function: Create GEANT4 TRD2-trapezoid *
1960 * 1960 * *
1961 * Input: Dx1 - half-length along X at -Dz 1961 * Input: Dx1 - half-length along X at -Dz 8----7 *
1962 * Dx2 - half-length along X ay +Dz 1962 * Dx2 - half-length along X ay +Dz 5----6 ! *
1963 * Dy1 - half-length along Y ay -Dz 1963 * Dy1 - half-length along Y ay -Dz ! 4-!--3 *
1964 * Dy2 - half-length along Y ay +Dz 1964 * Dy2 - half-length along Y ay +Dz 1----2 *
1965 * Dz - half-length along Z 1965 * Dz - half-length along Z *
1966 * 1966 * *
1967 ******************************************** 1967 ***********************************************************************/
1968 { 1968 {
1969 AllocateMemory(8,6); 1969 AllocateMemory(8,6);
1970 1970
1971 pV[1] = G4Point3D(-Dx1,-Dy1,-Dz); 1971 pV[1] = G4Point3D(-Dx1,-Dy1,-Dz);
1972 pV[2] = G4Point3D( Dx1,-Dy1,-Dz); 1972 pV[2] = G4Point3D( Dx1,-Dy1,-Dz);
1973 pV[3] = G4Point3D( Dx1, Dy1,-Dz); 1973 pV[3] = G4Point3D( Dx1, Dy1,-Dz);
1974 pV[4] = G4Point3D(-Dx1, Dy1,-Dz); 1974 pV[4] = G4Point3D(-Dx1, Dy1,-Dz);
1975 pV[5] = G4Point3D(-Dx2,-Dy2, Dz); 1975 pV[5] = G4Point3D(-Dx2,-Dy2, Dz);
1976 pV[6] = G4Point3D( Dx2,-Dy2, Dz); 1976 pV[6] = G4Point3D( Dx2,-Dy2, Dz);
1977 pV[7] = G4Point3D( Dx2, Dy2, Dz); 1977 pV[7] = G4Point3D( Dx2, Dy2, Dz);
1978 pV[8] = G4Point3D(-Dx2, Dy2, Dz); 1978 pV[8] = G4Point3D(-Dx2, Dy2, Dz);
1979 1979
1980 CreatePrism(); 1980 CreatePrism();
1981 } 1981 }
1982 1982
1983 HepPolyhedronTrd2::~HepPolyhedronTrd2() = def 1983 HepPolyhedronTrd2::~HepPolyhedronTrd2() = default;
1984 1984
1985 HepPolyhedronTrd1::HepPolyhedronTrd1(G4double 1985 HepPolyhedronTrd1::HepPolyhedronTrd1(G4double Dx1, G4double Dx2,
1986 G4double 1986 G4double Dy, G4double Dz)
1987 : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) { 1987 : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1988 1988
1989 HepPolyhedronTrd1::~HepPolyhedronTrd1() = def 1989 HepPolyhedronTrd1::~HepPolyhedronTrd1() = default;
1990 1990
1991 HepPolyhedronBox::HepPolyhedronBox(G4double D 1991 HepPolyhedronBox::HepPolyhedronBox(G4double Dx, G4double Dy, G4double Dz)
1992 : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {} 1992 : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1993 1993
1994 HepPolyhedronBox::~HepPolyhedronBox() = defau 1994 HepPolyhedronBox::~HepPolyhedronBox() = default;
1995 1995
1996 HepPolyhedronTrap::HepPolyhedronTrap(G4double 1996 HepPolyhedronTrap::HepPolyhedronTrap(G4double Dz,
1997 G4double 1997 G4double Theta,
1998 G4double 1998 G4double Phi,
1999 G4double 1999 G4double Dy1,
2000 G4double 2000 G4double Dx1,
2001 G4double 2001 G4double Dx2,
2002 G4double 2002 G4double Alp1,
2003 G4double 2003 G4double Dy2,
2004 G4double 2004 G4double Dx3,
2005 G4double 2005 G4double Dx4,
2006 G4double 2006 G4double Alp2)
2007 /******************************************** 2007 /***********************************************************************
2008 * 2008 * *
2009 * Name: HepPolyhedronTrap 2009 * Name: HepPolyhedronTrap Date: 20.11.96 *
2010 * Author: E.Chernyaev 2010 * Author: E.Chernyaev Revised: *
2011 * 2011 * *
2012 * Function: Create GEANT4 TRAP-trapezoid 2012 * Function: Create GEANT4 TRAP-trapezoid *
2013 * 2013 * *
2014 * Input: DZ - half-length in Z 2014 * Input: DZ - half-length in Z *
2015 * Theta,Phi - polar angles of the lin 2015 * Theta,Phi - polar angles of the line joining centres of the *
2016 * faces at Z=-Dz and Z=+D 2016 * faces at Z=-Dz and Z=+Dz *
2017 * Dy1 - half-length in Y of the face 2017 * Dy1 - half-length in Y of the face at Z=-Dz *
2018 * Dx1 - half-length in X of low edge 2018 * Dx1 - half-length in X of low edge of the face at Z=-Dz *
2019 * Dx2 - half-length in X of top edge 2019 * Dx2 - half-length in X of top edge of the face at Z=-Dz *
2020 * Alp1 - angle between Y-axis and the 2020 * Alp1 - angle between Y-axis and the median joining top and *
2021 * low edges of the face at Z=- 2021 * low edges of the face at Z=-Dz *
2022 * Dy2 - half-length in Y of the face 2022 * Dy2 - half-length in Y of the face at Z=+Dz *
2023 * Dx3 - half-length in X of low edge 2023 * Dx3 - half-length in X of low edge of the face at Z=+Dz *
2024 * Dx4 - half-length in X of top edge 2024 * Dx4 - half-length in X of top edge of the face at Z=+Dz *
2025 * Alp2 - angle between Y-axis and the 2025 * Alp2 - angle between Y-axis and the median joining top and *
2026 * low edges of the face at Z=+ 2026 * low edges of the face at Z=+Dz *
2027 * 2027 * *
2028 ******************************************** 2028 ***********************************************************************/
2029 { 2029 {
2030 G4double DzTthetaCphi = Dz*std::tan(Theta)* 2030 G4double DzTthetaCphi = Dz*std::tan(Theta)*std::cos(Phi);
2031 G4double DzTthetaSphi = Dz*std::tan(Theta)* 2031 G4double DzTthetaSphi = Dz*std::tan(Theta)*std::sin(Phi);
2032 G4double Dy1Talp1 = Dy1*std::tan(Alp1); 2032 G4double Dy1Talp1 = Dy1*std::tan(Alp1);
2033 G4double Dy2Talp2 = Dy2*std::tan(Alp2); 2033 G4double Dy2Talp2 = Dy2*std::tan(Alp2);
2034 2034
2035 AllocateMemory(8,6); 2035 AllocateMemory(8,6);
2036 2036
2037 pV[1] = G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx 2037 pV[1] = G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
2038 pV[2] = G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx 2038 pV[2] = G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
2039 pV[3] = G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx 2039 pV[3] = G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
2040 pV[4] = G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx 2040 pV[4] = G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
2041 pV[5] = G4Point3D( DzTthetaCphi-Dy2Talp2-Dx 2041 pV[5] = G4Point3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
2042 pV[6] = G4Point3D( DzTthetaCphi-Dy2Talp2+Dx 2042 pV[6] = G4Point3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
2043 pV[7] = G4Point3D( DzTthetaCphi+Dy2Talp2+Dx 2043 pV[7] = G4Point3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
2044 pV[8] = G4Point3D( DzTthetaCphi+Dy2Talp2-Dx 2044 pV[8] = G4Point3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
2045 2045
2046 CreatePrism(); 2046 CreatePrism();
2047 } 2047 }
2048 2048
2049 HepPolyhedronTrap::~HepPolyhedronTrap() = def 2049 HepPolyhedronTrap::~HepPolyhedronTrap() = default;
2050 2050
2051 HepPolyhedronPara::HepPolyhedronPara(G4double 2051 HepPolyhedronPara::HepPolyhedronPara(G4double Dx, G4double Dy, G4double Dz,
2052 G4double 2052 G4double Alpha, G4double Theta,
2053 G4double 2053 G4double Phi)
2054 : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, 2054 : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
2055 2055
2056 HepPolyhedronPara::~HepPolyhedronPara() = def 2056 HepPolyhedronPara::~HepPolyhedronPara() = default;
2057 2057
2058 HepPolyhedronParaboloid::HepPolyhedronParabol 2058 HepPolyhedronParaboloid::HepPolyhedronParaboloid(G4double r1,
2059 2059 G4double r2,
2060 2060 G4double dz,
2061 2061 G4double sPhi,
2062 2062 G4double dPhi)
2063 /******************************************** 2063 /***********************************************************************
2064 * 2064 * *
2065 * Name: HepPolyhedronParaboloid 2065 * Name: HepPolyhedronParaboloid Date: 28.06.07 *
2066 * Author: L.Lindroos, T.Nikitina (CERN), Jul 2066 * Author: L.Lindroos, T.Nikitina (CERN), July 2007 Revised: 28.06.07 *
2067 * 2067 * *
2068 * Function: Constructor for paraboloid 2068 * Function: Constructor for paraboloid *
2069 * 2069 * *
2070 * Input: r1 - inside and outside radiuses 2070 * Input: r1 - inside and outside radiuses at -Dz *
2071 * r2 - inside and outside radiuses 2071 * r2 - inside and outside radiuses at +Dz *
2072 * dz - half length in Z 2072 * dz - half length in Z *
2073 * sPhi - starting angle of the segme 2073 * sPhi - starting angle of the segment *
2074 * dPhi - segment range 2074 * dPhi - segment range *
2075 * 2075 * *
2076 ******************************************** 2076 ***********************************************************************/
2077 { 2077 {
2078 static const G4double wholeCircle=twopi; 2078 static const G4double wholeCircle=twopi;
2079 2079
2080 // C H E C K I N P U T P A R A M E T 2080 // C H E C K I N P U T P A R A M E T E R S
2081 2081
2082 G4int k = 0; 2082 G4int k = 0;
2083 if (r1 < 0. || r2 <= 0.) k = 1; 2083 if (r1 < 0. || r2 <= 0.) k = 1;
2084 2084
2085 if (dz <= 0.) k += 2; 2085 if (dz <= 0.) k += 2;
2086 2086
2087 G4double phi1, phi2, dphi; 2087 G4double phi1, phi2, dphi;
2088 2088
2089 if(dPhi < 0.) 2089 if(dPhi < 0.)
2090 { 2090 {
2091 phi2 = sPhi; phi1 = phi2 + dPhi; 2091 phi2 = sPhi; phi1 = phi2 + dPhi;
2092 } 2092 }
2093 else if(dPhi == 0.) 2093 else if(dPhi == 0.)
2094 { 2094 {
2095 phi1 = sPhi; phi2 = phi1 + wholeCircle; 2095 phi1 = sPhi; phi2 = phi1 + wholeCircle;
2096 } 2096 }
2097 else 2097 else
2098 { 2098 {
2099 phi1 = sPhi; phi2 = phi1 + dPhi; 2099 phi1 = sPhi; phi2 = phi1 + dPhi;
2100 } 2100 }
2101 dphi = phi2 - phi1; 2101 dphi = phi2 - phi1;
2102 2102
2103 if (std::abs(dphi-wholeCircle) < perMillion 2103 if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
2104 if (dphi > wholeCircle) k += 4; 2104 if (dphi > wholeCircle) k += 4;
2105 2105
2106 if (k != 0) { 2106 if (k != 0) {
2107 std::cerr << "HepPolyhedronParaboloid: er 2107 std::cerr << "HepPolyhedronParaboloid: error in input parameters";
2108 if ((k & 1) != 0) std::cerr << " (radiuse 2108 if ((k & 1) != 0) std::cerr << " (radiuses)";
2109 if ((k & 2) != 0) std::cerr << " (half-le 2109 if ((k & 2) != 0) std::cerr << " (half-length)";
2110 if ((k & 4) != 0) std::cerr << " (angles) 2110 if ((k & 4) != 0) std::cerr << " (angles)";
2111 std::cerr << std::endl; 2111 std::cerr << std::endl;
2112 std::cerr << " r1=" << r1; 2112 std::cerr << " r1=" << r1;
2113 std::cerr << " r2=" << r2; 2113 std::cerr << " r2=" << r2;
2114 std::cerr << " dz=" << dz << " sPhi=" << 2114 std::cerr << " dz=" << dz << " sPhi=" << sPhi << " dPhi=" << dPhi
2115 << std::endl; 2115 << std::endl;
2116 return; 2116 return;
2117 } 2117 }
2118 2118
2119 // P R E P A R E T W O P O L Y L I N 2119 // P R E P A R E T W O P O L Y L I N E S
2120 2120
2121 G4int n = GetNumberOfRotationSteps(); 2121 G4int n = GetNumberOfRotationSteps();
2122 G4double dl = (r2 - r1) / n; 2122 G4double dl = (r2 - r1) / n;
2123 G4double k1 = (r2*r2 - r1*r1) / 2 / dz; 2123 G4double k1 = (r2*r2 - r1*r1) / 2 / dz;
2124 G4double k2 = (r2*r2 + r1*r1) / 2; 2124 G4double k2 = (r2*r2 + r1*r1) / 2;
2125 2125
2126 auto zz = new G4double[n + 2], rr = new G4d 2126 auto zz = new G4double[n + 2], rr = new G4double[n + 2];
2127 2127
2128 zz[0] = dz; 2128 zz[0] = dz;
2129 rr[0] = r2; 2129 rr[0] = r2;
2130 2130
2131 for(G4int i = 1; i < n - 1; i++) 2131 for(G4int i = 1; i < n - 1; i++)
2132 { 2132 {
2133 rr[i] = rr[i-1] - dl; 2133 rr[i] = rr[i-1] - dl;
2134 zz[i] = (rr[i]*rr[i] - k2) / k1; 2134 zz[i] = (rr[i]*rr[i] - k2) / k1;
2135 if(rr[i] < 0) 2135 if(rr[i] < 0)
2136 { 2136 {
2137 rr[i] = 0; 2137 rr[i] = 0;
2138 zz[i] = 0; 2138 zz[i] = 0;
2139 } 2139 }
2140 } 2140 }
2141 2141
2142 zz[n-1] = -dz; 2142 zz[n-1] = -dz;
2143 rr[n-1] = r1; 2143 rr[n-1] = r1;
2144 2144
2145 zz[n] = dz; 2145 zz[n] = dz;
2146 rr[n] = 0; 2146 rr[n] = 0;
2147 2147
2148 zz[n+1] = -dz; 2148 zz[n+1] = -dz;
2149 rr[n+1] = 0; 2149 rr[n+1] = 0;
2150 2150
2151 // R O T A T E P O L Y L I N E S 2151 // R O T A T E P O L Y L I N E S
2152 2152
2153 RotateAroundZ(0, phi1, dphi, n, 2, zz, rr, 2153 RotateAroundZ(0, phi1, dphi, n, 2, zz, rr, -1, -1);
2154 SetReferences(); 2154 SetReferences();
2155 2155
2156 delete [] zz; 2156 delete [] zz;
2157 delete [] rr; 2157 delete [] rr;
2158 } 2158 }
2159 2159
2160 HepPolyhedronParaboloid::~HepPolyhedronParabo 2160 HepPolyhedronParaboloid::~HepPolyhedronParaboloid() = default;
2161 2161
2162 HepPolyhedronHype::HepPolyhedronHype(G4double 2162 HepPolyhedronHype::HepPolyhedronHype(G4double r1,
2163 G4double 2163 G4double r2,
2164 G4double 2164 G4double sqrtan1,
2165 G4double 2165 G4double sqrtan2,
2166 G4double 2166 G4double halfZ)
2167 /******************************************** 2167 /***********************************************************************
2168 * 2168 * *
2169 * Name: HepPolyhedronHype 2169 * Name: HepPolyhedronHype Date: 14.04.08 *
2170 * Author: Tatiana Nikitina (CERN) 2170 * Author: Tatiana Nikitina (CERN) Revised: 14.04.08 *
2171 * Evgueni Tcherniaev 2171 * Evgueni Tcherniaev 01.12.20 *
2172 * 2172 * *
2173 * Function: Constructor for Hype 2173 * Function: Constructor for Hype *
2174 * 2174 * *
2175 * Input: r1 - inside radius at z=0 2175 * Input: r1 - inside radius at z=0 *
2176 * r2 - outside radiuses at z=0 2176 * r2 - outside radiuses at z=0 *
2177 * sqrtan1 - sqr of tan of Inner Ster 2177 * sqrtan1 - sqr of tan of Inner Stereo Angle *
2178 * sqrtan2 - sqr of tan of Outer Ster 2178 * sqrtan2 - sqr of tan of Outer Stereo Angle *
2179 * halfZ - half length in Z 2179 * halfZ - half length in Z *
2180 * 2180 * *
2181 ******************************************** 2181 ***********************************************************************/
2182 { 2182 {
2183 static const G4double wholeCircle = twopi; 2183 static const G4double wholeCircle = twopi;
2184 2184
2185 // C H E C K I N P U T P A R A M E T 2185 // C H E C K I N P U T P A R A M E T E R S
2186 2186
2187 G4int k = 0; 2187 G4int k = 0;
2188 if (r1 < 0. || r2 < 0. || r1 >= r2) k = 1; 2188 if (r1 < 0. || r2 < 0. || r1 >= r2) k = 1;
2189 if (halfZ <= 0.) k += 2; 2189 if (halfZ <= 0.) k += 2;
2190 if (sqrtan1 < 0.|| sqrtan2 < 0.) k += 4; 2190 if (sqrtan1 < 0.|| sqrtan2 < 0.) k += 4;
2191 2191
2192 if (k != 0) 2192 if (k != 0)
2193 { 2193 {
2194 std::cerr << "HepPolyhedronHype: error in 2194 std::cerr << "HepPolyhedronHype: error in input parameters";
2195 if ((k & 1) != 0) std::cerr << " (radiuse 2195 if ((k & 1) != 0) std::cerr << " (radiuses)";
2196 if ((k & 2) != 0) std::cerr << " (half-le 2196 if ((k & 2) != 0) std::cerr << " (half-length)";
2197 if ((k & 4) != 0) std::cerr << " (angles) 2197 if ((k & 4) != 0) std::cerr << " (angles)";
2198 std::cerr << std::endl; 2198 std::cerr << std::endl;
2199 std::cerr << " r1=" << r1 << " r2=" << r2 2199 std::cerr << " r1=" << r1 << " r2=" << r2;
2200 std::cerr << " halfZ=" << halfZ << " sqrT 2200 std::cerr << " halfZ=" << halfZ << " sqrTan1=" << sqrtan1
2201 << " sqrTan2=" << sqrtan2 2201 << " sqrTan2=" << sqrtan2
2202 << std::endl; 2202 << std::endl;
2203 return; 2203 return;
2204 } 2204 }
2205 2205
2206 // P R E P A R E T W O P O L Y L I N 2206 // P R E P A R E T W O P O L Y L I N E S
2207 2207
2208 G4int ns = std::max(3, GetNumberOfRotationS 2208 G4int ns = std::max(3, GetNumberOfRotationSteps()/4);
2209 G4int nz1 = (sqrtan1 == 0.) ? 2 : ns + 1; 2209 G4int nz1 = (sqrtan1 == 0.) ? 2 : ns + 1;
2210 G4int nz2 = (sqrtan2 == 0.) ? 2 : ns + 1; 2210 G4int nz2 = (sqrtan2 == 0.) ? 2 : ns + 1;
2211 auto zz = new G4double[nz1 + nz2]; 2211 auto zz = new G4double[nz1 + nz2];
2212 auto rr = new G4double[nz1 + nz2]; 2212 auto rr = new G4double[nz1 + nz2];
2213 2213
2214 // external polyline 2214 // external polyline
2215 G4double dz2 = 2.*halfZ/(nz2 - 1); 2215 G4double dz2 = 2.*halfZ/(nz2 - 1);
2216 for(G4int i = 0; i < nz2; ++i) 2216 for(G4int i = 0; i < nz2; ++i)
2217 { 2217 {
2218 zz[i] = halfZ - dz2*i; 2218 zz[i] = halfZ - dz2*i;
2219 rr[i] = std::sqrt(sqrtan2*zz[i]*zz[i] + r 2219 rr[i] = std::sqrt(sqrtan2*zz[i]*zz[i] + r2*r2);
2220 } 2220 }
2221 2221
2222 // internal polyline 2222 // internal polyline
2223 G4double dz1 = 2.*halfZ/(nz1 - 1); 2223 G4double dz1 = 2.*halfZ/(nz1 - 1);
2224 for(G4int i = 0; i < nz1; ++i) 2224 for(G4int i = 0; i < nz1; ++i)
2225 { 2225 {
2226 G4int j = nz2 + i; 2226 G4int j = nz2 + i;
2227 zz[j] = halfZ - dz1*i; 2227 zz[j] = halfZ - dz1*i;
2228 rr[j] = std::sqrt(sqrtan1*zz[j]*zz[j] + r 2228 rr[j] = std::sqrt(sqrtan1*zz[j]*zz[j] + r1*r1);
2229 } 2229 }
2230 2230
2231 // R O T A T E P O L Y L I N E S 2231 // R O T A T E P O L Y L I N E S
2232 2232
2233 RotateAroundZ(0, 0., wholeCircle, nz2, nz1, 2233 RotateAroundZ(0, 0., wholeCircle, nz2, nz1, zz, rr, -1, -1);
2234 SetReferences(); 2234 SetReferences();
2235 2235
2236 delete [] zz; 2236 delete [] zz;
2237 delete [] rr; 2237 delete [] rr;
2238 } 2238 }
2239 2239
2240 HepPolyhedronHype::~HepPolyhedronHype() = def 2240 HepPolyhedronHype::~HepPolyhedronHype() = default;
2241 2241
2242 HepPolyhedronCons::HepPolyhedronCons(G4double 2242 HepPolyhedronCons::HepPolyhedronCons(G4double Rmn1,
2243 G4double 2243 G4double Rmx1,
2244 G4double 2244 G4double Rmn2,
2245 G4double 2245 G4double Rmx2,
2246 G4double 2246 G4double Dz,
2247 G4double 2247 G4double Phi1,
2248 G4double 2248 G4double Dphi)
2249 /******************************************** 2249 /***********************************************************************
2250 * 2250 * *
2251 * Name: HepPolyhedronCons::HepPolyhedronCons 2251 * Name: HepPolyhedronCons::HepPolyhedronCons Date: 15.12.96 *
2252 * Author: E.Chernyaev (IHEP/Protvino) 2252 * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
2253 * 2253 * *
2254 * Function: Constructor for CONS, TUBS, CONE 2254 * Function: Constructor for CONS, TUBS, CONE, TUBE *
2255 * 2255 * *
2256 * Input: Rmn1, Rmx1 - inside and outside rad 2256 * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
2257 * Rmn2, Rmx2 - inside and outside rad 2257 * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
2258 * Dz - half length in Z 2258 * Dz - half length in Z *
2259 * Phi1 - starting angle of the 2259 * Phi1 - starting angle of the segment *
2260 * Dphi - segment range 2260 * Dphi - segment range *
2261 * 2261 * *
2262 ******************************************** 2262 ***********************************************************************/
2263 { 2263 {
2264 static const G4double wholeCircle=twopi; 2264 static const G4double wholeCircle=twopi;
2265 2265
2266 // C H E C K I N P U T P A R A M E T 2266 // C H E C K I N P U T P A R A M E T E R S
2267 2267
2268 G4int k = 0; 2268 G4int k = 0;
2269 if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || 2269 if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
2270 if (Rmn1 > Rmx1 || Rmn2 > Rmx2) 2270 if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
2271 if (Rmn1 == Rmx1 && Rmn2 == Rmx2) 2271 if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
2272 2272
2273 if (Dz <= 0.) k += 2; 2273 if (Dz <= 0.) k += 2;
2274 2274
2275 G4double phi1, phi2, dphi; 2275 G4double phi1, phi2, dphi;
2276 if (Dphi < 0.) { 2276 if (Dphi < 0.) {
2277 phi2 = Phi1; phi1 = phi2 - Dphi; 2277 phi2 = Phi1; phi1 = phi2 - Dphi;
2278 }else if (Dphi == 0.) { 2278 }else if (Dphi == 0.) {
2279 phi1 = Phi1; phi2 = phi1 + wholeCircle; 2279 phi1 = Phi1; phi2 = phi1 + wholeCircle;
2280 }else{ 2280 }else{
2281 phi1 = Phi1; phi2 = phi1 + Dphi; 2281 phi1 = Phi1; phi2 = phi1 + Dphi;
2282 } 2282 }
2283 dphi = phi2 - phi1; 2283 dphi = phi2 - phi1;
2284 if (std::abs(dphi-wholeCircle) < perMillion 2284 if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
2285 if (dphi > wholeCircle) k += 4; 2285 if (dphi > wholeCircle) k += 4;
2286 2286
2287 if (k != 0) { 2287 if (k != 0) {
2288 std::cerr << "HepPolyhedronCone(s)/Tube(s 2288 std::cerr << "HepPolyhedronCone(s)/Tube(s): error in input parameters";
2289 if ((k & 1) != 0) std::cerr << " (radiuse 2289 if ((k & 1) != 0) std::cerr << " (radiuses)";
2290 if ((k & 2) != 0) std::cerr << " (half-le 2290 if ((k & 2) != 0) std::cerr << " (half-length)";
2291 if ((k & 4) != 0) std::cerr << " (angles) 2291 if ((k & 4) != 0) std::cerr << " (angles)";
2292 std::cerr << std::endl; 2292 std::cerr << std::endl;
2293 std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" 2293 std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
2294 std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" 2294 std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
2295 std::cerr << " Dz=" << Dz << " Phi1=" << 2295 std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
2296 << std::endl; 2296 << std::endl;
2297 return; 2297 return;
2298 } 2298 }
2299 2299
2300 // P R E P A R E T W O P O L Y L I N 2300 // P R E P A R E T W O P O L Y L I N E S
2301 2301
2302 G4double zz[4], rr[4]; 2302 G4double zz[4], rr[4];
2303 zz[0] = Dz; 2303 zz[0] = Dz;
2304 zz[1] = -Dz; 2304 zz[1] = -Dz;
2305 zz[2] = Dz; 2305 zz[2] = Dz;
2306 zz[3] = -Dz; 2306 zz[3] = -Dz;
2307 rr[0] = Rmx2; 2307 rr[0] = Rmx2;
2308 rr[1] = Rmx1; 2308 rr[1] = Rmx1;
2309 rr[2] = Rmn2; 2309 rr[2] = Rmn2;
2310 rr[3] = Rmn1; 2310 rr[3] = Rmn1;
2311 2311
2312 // R O T A T E P O L Y L I N E S 2312 // R O T A T E P O L Y L I N E S
2313 2313
2314 RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, 2314 RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1);
2315 SetReferences(); 2315 SetReferences();
2316 } 2316 }
2317 2317
2318 HepPolyhedronCons::~HepPolyhedronCons() = def 2318 HepPolyhedronCons::~HepPolyhedronCons() = default;
2319 2319
2320 HepPolyhedronCone::HepPolyhedronCone(G4double 2320 HepPolyhedronCone::HepPolyhedronCone(G4double Rmn1, G4double Rmx1,
2321 G4double 2321 G4double Rmn2, G4double Rmx2,
2322 G4double 2322 G4double Dz) :
2323 HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, D 2323 HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
2324 2324
2325 HepPolyhedronCone::~HepPolyhedronCone() = def 2325 HepPolyhedronCone::~HepPolyhedronCone() = default;
2326 2326
2327 HepPolyhedronTubs::HepPolyhedronTubs(G4double 2327 HepPolyhedronTubs::HepPolyhedronTubs(G4double Rmin, G4double Rmax,
2328 G4double 2328 G4double Dz,
2329 G4double 2329 G4double Phi1, G4double Dphi)
2330 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rma 2330 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
2331 2331
2332 HepPolyhedronTubs::~HepPolyhedronTubs() = def 2332 HepPolyhedronTubs::~HepPolyhedronTubs() = default;
2333 2333
2334 HepPolyhedronTube::HepPolyhedronTube (G4doubl 2334 HepPolyhedronTube::HepPolyhedronTube (G4double Rmin, G4double Rmax,
2335 G4doubl 2335 G4double Dz)
2336 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, 2336 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
2337 2337
2338 HepPolyhedronTube::~HepPolyhedronTube () = de 2338 HepPolyhedronTube::~HepPolyhedronTube () = default;
2339 2339
2340 HepPolyhedronPgon::HepPolyhedronPgon(G4double 2340 HepPolyhedronPgon::HepPolyhedronPgon(G4double phi,
2341 G4double 2341 G4double dphi,
2342 G4int np 2342 G4int npdv,
2343 G4int nz 2343 G4int nz,
2344 const G4 2344 const G4double *z,
2345 const G4 2345 const G4double *rmin,
2346 const G4 2346 const G4double *rmax)
2347 /******************************************** 2347 /***********************************************************************
2348 * 2348 * *
2349 * Name: HepPolyhedronPgon 2349 * Name: HepPolyhedronPgon Date: 09.12.96 *
2350 * Author: E.Chernyaev 2350 * Author: E.Chernyaev Revised: *
2351 * 2351 * *
2352 * Function: Constructor of polyhedron for PG 2352 * Function: Constructor of polyhedron for PGON, PCON *
2353 * 2353 * *
2354 * Input: phi - initial phi 2354 * Input: phi - initial phi *
2355 * dphi - delta phi 2355 * dphi - delta phi *
2356 * npdv - number of steps along phi 2356 * npdv - number of steps along phi *
2357 * nz - number of z-planes (at least 2357 * nz - number of z-planes (at least two) *
2358 * z[] - z coordinates of the slices 2358 * z[] - z coordinates of the slices *
2359 * rmin[] - smaller r at the slices 2359 * rmin[] - smaller r at the slices *
2360 * rmax[] - bigger r at the slices 2360 * rmax[] - bigger r at the slices *
2361 * 2361 * *
2362 ******************************************** 2362 ***********************************************************************/
2363 { 2363 {
2364 // C H E C K I N P U T P A R A M E T 2364 // C H E C K I N P U T P A R A M E T E R S
2365 2365
2366 if (dphi <= 0. || dphi > twopi) { 2366 if (dphi <= 0. || dphi > twopi) {
2367 std::cerr 2367 std::cerr
2368 << "HepPolyhedronPgon/Pcon: wrong delta 2368 << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
2369 << std::endl; 2369 << std::endl;
2370 return; 2370 return;
2371 } 2371 }
2372 2372
2373 if (nz < 2) { 2373 if (nz < 2) {
2374 std::cerr 2374 std::cerr
2375 << "HepPolyhedronPgon/Pcon: number of z 2375 << "HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
2376 << std::endl; 2376 << std::endl;
2377 return; 2377 return;
2378 } 2378 }
2379 2379
2380 if (npdv < 0) { 2380 if (npdv < 0) {
2381 std::cerr 2381 std::cerr
2382 << "HepPolyhedronPgon/Pcon: error in nu 2382 << "HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
2383 << std::endl; 2383 << std::endl;
2384 return; 2384 return;
2385 } 2385 }
2386 2386
2387 G4int i; 2387 G4int i;
2388 for (i=0; i<nz; i++) { 2388 for (i=0; i<nz; i++) {
2389 if (rmin[i] < 0. || rmax[i] < 0. || rmin[ 2389 if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
2390 std::cerr 2390 std::cerr
2391 << "HepPolyhedronPgon: error in radiu 2391 << "HepPolyhedronPgon: error in radiuses rmin[" << i << "]="
2392 << rmin[i] << " rmax[" << i << "]=" < 2392 << rmin[i] << " rmax[" << i << "]=" << rmax[i]
2393 << std::endl; 2393 << std::endl;
2394 return; 2394 return;
2395 } 2395 }
2396 } 2396 }
2397 2397
2398 // P R E P A R E T W O P O L Y L I N 2398 // P R E P A R E T W O P O L Y L I N E S
2399 2399
2400 G4double *zz, *rr; 2400 G4double *zz, *rr;
2401 zz = new G4double[2*nz]; 2401 zz = new G4double[2*nz];
2402 rr = new G4double[2*nz]; 2402 rr = new G4double[2*nz];
2403 2403
2404 if (z[0] > z[nz-1]) { 2404 if (z[0] > z[nz-1]) {
2405 for (i=0; i<nz; i++) { 2405 for (i=0; i<nz; i++) {
2406 zz[i] = z[i]; 2406 zz[i] = z[i];
2407 rr[i] = rmax[i]; 2407 rr[i] = rmax[i];
2408 zz[i+nz] = z[i]; 2408 zz[i+nz] = z[i];
2409 rr[i+nz] = rmin[i]; 2409 rr[i+nz] = rmin[i];
2410 } 2410 }
2411 }else{ 2411 }else{
2412 for (i=0; i<nz; i++) { 2412 for (i=0; i<nz; i++) {
2413 zz[i] = z[nz-i-1]; 2413 zz[i] = z[nz-i-1];
2414 rr[i] = rmax[nz-i-1]; 2414 rr[i] = rmax[nz-i-1];
2415 zz[i+nz] = z[nz-i-1]; 2415 zz[i+nz] = z[nz-i-1];
2416 rr[i+nz] = rmin[nz-i-1]; 2416 rr[i+nz] = rmin[nz-i-1];
2417 } 2417 }
2418 } 2418 }
2419 2419
2420 // R O T A T E P O L Y L I N E S 2420 // R O T A T E P O L Y L I N E S
2421 2421
2422 G4int nodeVis = 1; 2422 G4int nodeVis = 1;
2423 G4int edgeVis = (npdv == 0) ? -1 : 1; 2423 G4int edgeVis = (npdv == 0) ? -1 : 1;
2424 RotateAroundZ(npdv, phi, dphi, nz, nz, zz, 2424 RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, nodeVis, edgeVis);
2425 SetReferences(); 2425 SetReferences();
2426 2426
2427 delete [] zz; 2427 delete [] zz;
2428 delete [] rr; 2428 delete [] rr;
2429 } 2429 }
2430 2430
2431 HepPolyhedronPgon::HepPolyhedronPgon(G4double 2431 HepPolyhedronPgon::HepPolyhedronPgon(G4double phi,
2432 G4double 2432 G4double dphi,
2433 G4int np 2433 G4int npdv,
2434 const st 2434 const std::vector<G4TwoVector> &rz)
2435 /******************************************** 2435 /***********************************************************************
2436 * 2436 * *
2437 * Name: HepPolyhedronPgon 2437 * Name: HepPolyhedronPgon Date: 12.05.21 *
2438 * Author: E.Tcherniaev (E.Chernyaev) 2438 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2439 * 2439 * *
2440 * Function: Constructor of polyhedron for PG 2440 * Function: Constructor of polyhedron for PGON, PCON *
2441 * 2441 * *
2442 * Input: phi - initial phi 2442 * Input: phi - initial phi *
2443 * dphi - delta phi 2443 * dphi - delta phi *
2444 * npdv - number of steps along phi 2444 * npdv - number of steps along phi *
2445 * rz - rz-contour 2445 * rz - rz-contour *
2446 * 2446 * *
2447 ******************************************** 2447 ***********************************************************************/
2448 { 2448 {
2449 // C H E C K I N P U T P A R A M E T 2449 // C H E C K I N P U T P A R A M E T E R S
2450 2450
2451 if (dphi <= 0. || dphi > twopi) { 2451 if (dphi <= 0. || dphi > twopi) {
2452 std::cerr 2452 std::cerr
2453 << "HepPolyhedronPgon/Pcon: wrong delta 2453 << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
2454 << std::endl; 2454 << std::endl;
2455 return; 2455 return;
2456 } 2456 }
2457 2457
2458 if (npdv < 0) { 2458 if (npdv < 0) {
2459 std::cerr 2459 std::cerr
2460 << "HepPolyhedronPgon/Pcon: error in nu 2460 << "HepPolyhedronPgon/Pcon: error in number of phi-steps = " << npdv
2461 << std::endl; 2461 << std::endl;
2462 return; 2462 return;
2463 } 2463 }
2464 2464
2465 G4int nrz = (G4int)rz.size(); 2465 G4int nrz = (G4int)rz.size();
2466 if (nrz < 3) { 2466 if (nrz < 3) {
2467 std::cerr 2467 std::cerr
2468 << "HepPolyhedronPgon/Pcon: invalid num 2468 << "HepPolyhedronPgon/Pcon: invalid number of nodes in rz-contour = " << nrz
2469 << std::endl; 2469 << std::endl;
2470 return; 2470 return;
2471 } 2471 }
2472 2472
2473 // R O T A T E P O L Y L I N E 2473 // R O T A T E P O L Y L I N E
2474 2474
2475 G4int nodeVis = 1; 2475 G4int nodeVis = 1;
2476 G4int edgeVis = (npdv == 0) ? -1 : 1; 2476 G4int edgeVis = (npdv == 0) ? -1 : 1;
2477 RotateContourAroundZ(npdv, phi, dphi, rz, n 2477 RotateContourAroundZ(npdv, phi, dphi, rz, nodeVis, edgeVis);
2478 SetReferences(); 2478 SetReferences();
2479 } 2479 }
2480 2480
2481 HepPolyhedronPgon::~HepPolyhedronPgon() = def 2481 HepPolyhedronPgon::~HepPolyhedronPgon() = default;
2482 2482
2483 HepPolyhedronPcon::HepPolyhedronPcon(G4double 2483 HepPolyhedronPcon::HepPolyhedronPcon(G4double phi, G4double dphi, G4int nz,
2484 const G4 2484 const G4double *z,
2485 const G4 2485 const G4double *rmin,
2486 const G4 2486 const G4double *rmax)
2487 : HepPolyhedronPgon(phi, dphi, 0, nz, z, rm 2487 : HepPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
2488 2488
2489 HepPolyhedronPcon::HepPolyhedronPcon(G4double 2489 HepPolyhedronPcon::HepPolyhedronPcon(G4double phi, G4double dphi,
2490 const st 2490 const std::vector<G4TwoVector> &rz)
2491 : HepPolyhedronPgon(phi, dphi, 0, rz) {} 2491 : HepPolyhedronPgon(phi, dphi, 0, rz) {}
2492 2492
2493 HepPolyhedronPcon::~HepPolyhedronPcon() = def 2493 HepPolyhedronPcon::~HepPolyhedronPcon() = default;
2494 2494
2495 HepPolyhedronSphere::HepPolyhedronSphere(G4do 2495 HepPolyhedronSphere::HepPolyhedronSphere(G4double rmin, G4double rmax,
2496 G4do 2496 G4double phi, G4double dphi,
2497 G4do 2497 G4double the, G4double dthe)
2498 /******************************************** 2498 /***********************************************************************
2499 * 2499 * *
2500 * Name: HepPolyhedronSphere 2500 * Name: HepPolyhedronSphere Date: 11.12.96 *
2501 * Author: E.Chernyaev (IHEP/Protvino) 2501 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
2502 * 2502 * *
2503 * Function: Constructor of polyhedron for SP 2503 * Function: Constructor of polyhedron for SPHERE *
2504 * 2504 * *
2505 * Input: rmin - internal radius 2505 * Input: rmin - internal radius *
2506 * rmax - external radius 2506 * rmax - external radius *
2507 * phi - initial phi 2507 * phi - initial phi *
2508 * dphi - delta phi 2508 * dphi - delta phi *
2509 * the - initial theta 2509 * the - initial theta *
2510 * dthe - delta theta 2510 * dthe - delta theta *
2511 * 2511 * *
2512 ******************************************** 2512 ***********************************************************************/
2513 { 2513 {
2514 // C H E C K I N P U T P A R A M E T 2514 // C H E C K I N P U T P A R A M E T E R S
2515 2515
2516 if (dphi <= 0. || dphi > twopi) { 2516 if (dphi <= 0. || dphi > twopi) {
2517 std::cerr 2517 std::cerr
2518 << "HepPolyhedronSphere: wrong delta ph 2518 << "HepPolyhedronSphere: wrong delta phi = " << dphi
2519 << std::endl; 2519 << std::endl;
2520 return; 2520 return;
2521 } 2521 }
2522 2522
2523 if (the < 0. || the > pi) { 2523 if (the < 0. || the > pi) {
2524 std::cerr 2524 std::cerr
2525 << "HepPolyhedronSphere: wrong theta = 2525 << "HepPolyhedronSphere: wrong theta = " << the
2526 << std::endl; 2526 << std::endl;
2527 return; 2527 return;
2528 } 2528 }
2529 2529
2530 if (dthe <= 0. || dthe > pi) { 2530 if (dthe <= 0. || dthe > pi) {
2531 std::cerr 2531 std::cerr
2532 << "HepPolyhedronSphere: wrong delta th 2532 << "HepPolyhedronSphere: wrong delta theta = " << dthe
2533 << std::endl; 2533 << std::endl;
2534 return; 2534 return;
2535 } 2535 }
2536 2536
2537 if (the+dthe > pi) { 2537 if (the+dthe > pi) {
2538 std::cerr 2538 std::cerr
2539 << "HepPolyhedronSphere: wrong theta + 2539 << "HepPolyhedronSphere: wrong theta + delta theta = "
2540 << the << " " << dthe 2540 << the << " " << dthe
2541 << std::endl; 2541 << std::endl;
2542 return; 2542 return;
2543 } 2543 }
2544 2544
2545 if (rmin < 0. || rmin >= rmax) { 2545 if (rmin < 0. || rmin >= rmax) {
2546 std::cerr 2546 std::cerr
2547 << "HepPolyhedronSphere: error in radiu 2547 << "HepPolyhedronSphere: error in radiuses"
2548 << " rmin=" << rmin << " rmax=" << rmax 2548 << " rmin=" << rmin << " rmax=" << rmax
2549 << std::endl; 2549 << std::endl;
2550 return; 2550 return;
2551 } 2551 }
2552 2552
2553 // P R E P A R E T W O P O L Y L I N 2553 // P R E P A R E T W O P O L Y L I N E S
2554 2554
2555 G4int nds = (GetNumberOfRotationSteps() + 1 2555 G4int nds = (GetNumberOfRotationSteps() + 1) / 2;
2556 G4int np1 = G4int(dthe*nds/pi+.5) + 1; 2556 G4int np1 = G4int(dthe*nds/pi+.5) + 1;
2557 if (np1 <= 1) np1 = 2; 2557 if (np1 <= 1) np1 = 2;
2558 G4int np2 = rmin < spatialTolerance ? 1 : n 2558 G4int np2 = rmin < spatialTolerance ? 1 : np1;
2559 2559
2560 G4double *zz, *rr; 2560 G4double *zz, *rr;
2561 zz = new G4double[np1+np2]; 2561 zz = new G4double[np1+np2];
2562 rr = new G4double[np1+np2]; 2562 rr = new G4double[np1+np2];
2563 2563
2564 G4double a = dthe/(np1-1); 2564 G4double a = dthe/(np1-1);
2565 G4double cosa, sina; 2565 G4double cosa, sina;
2566 for (G4int i=0; i<np1; i++) { 2566 for (G4int i=0; i<np1; i++) {
2567 cosa = std::cos(the+i*a); 2567 cosa = std::cos(the+i*a);
2568 sina = std::sin(the+i*a); 2568 sina = std::sin(the+i*a);
2569 zz[i] = rmax*cosa; 2569 zz[i] = rmax*cosa;
2570 rr[i] = rmax*sina; 2570 rr[i] = rmax*sina;
2571 if (np2 > 1) { 2571 if (np2 > 1) {
2572 zz[i+np1] = rmin*cosa; 2572 zz[i+np1] = rmin*cosa;
2573 rr[i+np1] = rmin*sina; 2573 rr[i+np1] = rmin*sina;
2574 } 2574 }
2575 } 2575 }
2576 if (np2 == 1) { 2576 if (np2 == 1) {
2577 zz[np1] = 0.; 2577 zz[np1] = 0.;
2578 rr[np1] = 0.; 2578 rr[np1] = 0.;
2579 } 2579 }
2580 2580
2581 // R O T A T E P O L Y L I N E S 2581 // R O T A T E P O L Y L I N E S
2582 2582
2583 RotateAroundZ(0, phi, dphi, np1, np2, zz, r 2583 RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
2584 SetReferences(); 2584 SetReferences();
2585 2585
2586 delete [] zz; 2586 delete [] zz;
2587 delete [] rr; 2587 delete [] rr;
2588 } 2588 }
2589 2589
2590 HepPolyhedronSphere::~HepPolyhedronSphere() = 2590 HepPolyhedronSphere::~HepPolyhedronSphere() = default;
2591 2591
2592 HepPolyhedronTorus::HepPolyhedronTorus(G4doub 2592 HepPolyhedronTorus::HepPolyhedronTorus(G4double rmin,
2593 G4doub 2593 G4double rmax,
2594 G4doub 2594 G4double rtor,
2595 G4doub 2595 G4double phi,
2596 G4doub 2596 G4double dphi)
2597 /******************************************** 2597 /***********************************************************************
2598 * 2598 * *
2599 * Name: HepPolyhedronTorus 2599 * Name: HepPolyhedronTorus Date: 11.12.96 *
2600 * Author: E.Chernyaev (IHEP/Protvino) 2600 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
2601 * 2601 * *
2602 * Function: Constructor of polyhedron for TO 2602 * Function: Constructor of polyhedron for TORUS *
2603 * 2603 * *
2604 * Input: rmin - internal radius 2604 * Input: rmin - internal radius *
2605 * rmax - external radius 2605 * rmax - external radius *
2606 * rtor - radius of torus 2606 * rtor - radius of torus *
2607 * phi - initial phi 2607 * phi - initial phi *
2608 * dphi - delta phi 2608 * dphi - delta phi *
2609 * 2609 * *
2610 ******************************************** 2610 ***********************************************************************/
2611 { 2611 {
2612 // C H E C K I N P U T P A R A M E T 2612 // C H E C K I N P U T P A R A M E T E R S
2613 2613
2614 if (dphi <= 0. || dphi > twopi) { 2614 if (dphi <= 0. || dphi > twopi) {
2615 std::cerr 2615 std::cerr
2616 << "HepPolyhedronTorus: wrong delta phi 2616 << "HepPolyhedronTorus: wrong delta phi = " << dphi
2617 << std::endl; 2617 << std::endl;
2618 return; 2618 return;
2619 } 2619 }
2620 2620
2621 if (rmin < 0. || rmin >= rmax || rmax >= rt 2621 if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2622 std::cerr 2622 std::cerr
2623 << "HepPolyhedronTorus: error in radius 2623 << "HepPolyhedronTorus: error in radiuses"
2624 << " rmin=" << rmin << " rmax=" << rmax 2624 << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
2625 << std::endl; 2625 << std::endl;
2626 return; 2626 return;
2627 } 2627 }
2628 2628
2629 // P R E P A R E T W O P O L Y L I N 2629 // P R E P A R E T W O P O L Y L I N E S
2630 2630
2631 G4int np1 = GetNumberOfRotationSteps(); 2631 G4int np1 = GetNumberOfRotationSteps();
2632 G4int np2 = rmin < spatialTolerance ? 1 : n 2632 G4int np2 = rmin < spatialTolerance ? 1 : np1;
2633 2633
2634 G4double *zz, *rr; 2634 G4double *zz, *rr;
2635 zz = new G4double[np1+np2]; 2635 zz = new G4double[np1+np2];
2636 rr = new G4double[np1+np2]; 2636 rr = new G4double[np1+np2];
2637 2637
2638 G4double a = twopi/np1; 2638 G4double a = twopi/np1;
2639 G4double cosa, sina; 2639 G4double cosa, sina;
2640 for (G4int i=0; i<np1; i++) { 2640 for (G4int i=0; i<np1; i++) {
2641 cosa = std::cos(i*a); 2641 cosa = std::cos(i*a);
2642 sina = std::sin(i*a); 2642 sina = std::sin(i*a);
2643 zz[i] = rmax*cosa; 2643 zz[i] = rmax*cosa;
2644 rr[i] = rtor+rmax*sina; 2644 rr[i] = rtor+rmax*sina;
2645 if (np2 > 1) { 2645 if (np2 > 1) {
2646 zz[i+np1] = rmin*cosa; 2646 zz[i+np1] = rmin*cosa;
2647 rr[i+np1] = rtor+rmin*sina; 2647 rr[i+np1] = rtor+rmin*sina;
2648 } 2648 }
2649 } 2649 }
2650 if (np2 == 1) { 2650 if (np2 == 1) {
2651 zz[np1] = 0.; 2651 zz[np1] = 0.;
2652 rr[np1] = rtor; 2652 rr[np1] = rtor;
2653 np2 = -1; 2653 np2 = -1;
2654 } 2654 }
2655 2655
2656 // R O T A T E P O L Y L I N E S 2656 // R O T A T E P O L Y L I N E S
2657 2657
2658 RotateAroundZ(0, phi, dphi, -np1, -np2, zz, 2658 RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1);
2659 SetReferences(); 2659 SetReferences();
2660 2660
2661 delete [] zz; 2661 delete [] zz;
2662 delete [] rr; 2662 delete [] rr;
2663 } 2663 }
2664 2664
2665 HepPolyhedronTorus::~HepPolyhedronTorus() = d 2665 HepPolyhedronTorus::~HepPolyhedronTorus() = default;
2666 2666
2667 HepPolyhedronTet::HepPolyhedronTet(const G4do 2667 HepPolyhedronTet::HepPolyhedronTet(const G4double p0[3],
2668 const G4do 2668 const G4double p1[3],
2669 const G4do 2669 const G4double p2[3],
2670 const G4do 2670 const G4double p3[3])
2671 /******************************************** 2671 /***********************************************************************
2672 * 2672 * *
2673 * Name: HepPolyhedronTet 2673 * Name: HepPolyhedronTet Date: 21.02.2020 *
2674 * Author: E.Tcherniaev (E.Chernyaev) 2674 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2675 * 2675 * *
2676 * Function: Constructor of polyhedron for TE 2676 * Function: Constructor of polyhedron for TETrahedron *
2677 * 2677 * *
2678 * Input: p0,p1,p2,p3 - vertices 2678 * Input: p0,p1,p2,p3 - vertices *
2679 * 2679 * *
2680 ******************************************** 2680 ***********************************************************************/
2681 { 2681 {
2682 AllocateMemory(4,4); 2682 AllocateMemory(4,4);
2683 2683
2684 pV[1].set(p0[0], p0[1], p0[2]); 2684 pV[1].set(p0[0], p0[1], p0[2]);
2685 pV[2].set(p1[0], p1[1], p1[2]); 2685 pV[2].set(p1[0], p1[1], p1[2]);
2686 pV[3].set(p2[0], p2[1], p2[2]); 2686 pV[3].set(p2[0], p2[1], p2[2]);
2687 pV[4].set(p3[0], p3[1], p3[2]); 2687 pV[4].set(p3[0], p3[1], p3[2]);
2688 2688
2689 G4Vector3D v1(pV[2] - pV[1]); 2689 G4Vector3D v1(pV[2] - pV[1]);
2690 G4Vector3D v2(pV[3] - pV[1]); 2690 G4Vector3D v2(pV[3] - pV[1]);
2691 G4Vector3D v3(pV[4] - pV[1]); 2691 G4Vector3D v3(pV[4] - pV[1]);
2692 2692
2693 if (v1.cross(v2).dot(v3) < 0.) 2693 if (v1.cross(v2).dot(v3) < 0.)
2694 { 2694 {
2695 pV[3].set(p3[0], p3[1], p3[2]); 2695 pV[3].set(p3[0], p3[1], p3[2]);
2696 pV[4].set(p2[0], p2[1], p2[2]); 2696 pV[4].set(p2[0], p2[1], p2[2]);
2697 } 2697 }
2698 2698
2699 pF[1] = G4Facet(1,2, 3,4, 2,3); 2699 pF[1] = G4Facet(1,2, 3,4, 2,3);
2700 pF[2] = G4Facet(1,3, 4,4, 3,1); 2700 pF[2] = G4Facet(1,3, 4,4, 3,1);
2701 pF[3] = G4Facet(1,1, 2,4, 4,2); 2701 pF[3] = G4Facet(1,1, 2,4, 4,2);
2702 pF[4] = G4Facet(2,1, 3,2, 4,3); 2702 pF[4] = G4Facet(2,1, 3,2, 4,3);
2703 } 2703 }
2704 2704
2705 HepPolyhedronTet::~HepPolyhedronTet() = defau 2705 HepPolyhedronTet::~HepPolyhedronTet() = default;
2706 2706
2707 HepPolyhedronEllipsoid::HepPolyhedronEllipsoi 2707 HepPolyhedronEllipsoid::HepPolyhedronEllipsoid(G4double ax, G4double by,
2708 2708 G4double cz, G4double zCut1,
2709 2709 G4double zCut2)
2710 /******************************************** 2710 /***********************************************************************
2711 * 2711 * *
2712 * Name: HepPolyhedronEllipsoid 2712 * Name: HepPolyhedronEllipsoid Date: 25.02.05 *
2713 * Author: G.Guerrieri 2713 * Author: G.Guerrieri Revised: *
2714 * Evgueni Tcherniaev 2714 * Evgueni Tcherniaev 20.01.22 *
2715 * 2715 * *
2716 * Function: Constructor of polyhedron for EL 2716 * Function: Constructor of polyhedron for ELLIPSOID *
2717 * 2717 * *
2718 * Input: ax - semiaxis x 2718 * Input: ax - semiaxis x *
2719 * by - semiaxis y 2719 * by - semiaxis y *
2720 * cz - semiaxis z 2720 * cz - semiaxis z *
2721 * zCut1 - lower cut plane level (soli 2721 * zCut1 - lower cut plane level (solid lies above this plane) *
2722 * zCut2 - upper cut plane level (soli 2722 * zCut2 - upper cut plane level (solid lies below this plane) *
2723 * 2723 * *
2724 ******************************************** 2724 ***********************************************************************/
2725 { 2725 {
2726 // C H E C K I N P U T P A R A M E T 2726 // C H E C K I N P U T P A R A M E T E R S
2727 2727
2728 if (zCut1 >= cz || zCut2 <= -cz || zCut1 > 2728 if (zCut1 >= cz || zCut2 <= -cz || zCut1 > zCut2) {
2729 std::cerr << "HepPolyhedronEllipsoid: wro 2729 std::cerr << "HepPolyhedronEllipsoid: wrong zCut1 = " << zCut1
2730 << " zCut2 = " << zCut2 2730 << " zCut2 = " << zCut2
2731 << " for given cz = " << cz << std 2731 << " for given cz = " << cz << std::endl;
2732 return; 2732 return;
2733 } 2733 }
2734 if (cz <= 0.0) { 2734 if (cz <= 0.0) {
2735 std::cerr << "HepPolyhedronEllipsoid: bad 2735 std::cerr << "HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
2736 << std::endl; 2736 << std::endl;
2737 return; 2737 return;
2738 } 2738 }
2739 2739
2740 // P R E P A R E T W O P O L Y L I N 2740 // P R E P A R E T W O P O L Y L I N E S
2741 // generate sphere of radius cz first, th 2741 // generate sphere of radius cz first, then rescale x and y later
2742 2742
2743 G4double sthe = std::acos(zCut2/cz); 2743 G4double sthe = std::acos(zCut2/cz);
2744 G4double dthe = std::acos(zCut1/cz) - sthe; 2744 G4double dthe = std::acos(zCut1/cz) - sthe;
2745 G4int nds = (GetNumberOfRotationSteps() + 1 2745 G4int nds = (GetNumberOfRotationSteps() + 1)/2;
2746 G4int np1 = G4int(dthe*nds/pi + 0.5) + 1; 2746 G4int np1 = G4int(dthe*nds/pi + 0.5) + 1;
2747 if (np1 <= 1) np1 = 2; 2747 if (np1 <= 1) np1 = 2;
2748 G4int np2 = 2; 2748 G4int np2 = 2;
2749 2749
2750 G4double *zz, *rr; 2750 G4double *zz, *rr;
2751 zz = new G4double[np1 + np2]; 2751 zz = new G4double[np1 + np2];
2752 rr = new G4double[np1 + np2]; 2752 rr = new G4double[np1 + np2];
2753 if ((zz == nullptr) || (rr == nullptr)) 2753 if ((zz == nullptr) || (rr == nullptr))
2754 { 2754 {
2755 G4Exception("HepPolyhedronEllipsoid::HepP 2755 G4Exception("HepPolyhedronEllipsoid::HepPolyhedronEllipsoid",
2756 "greps1002", FatalException, 2756 "greps1002", FatalException, "Out of memory");
2757 } 2757 }
2758 2758
2759 G4double a = dthe/(np1 - 1); 2759 G4double a = dthe/(np1 - 1);
2760 G4double cosa, sina; 2760 G4double cosa, sina;
2761 for (G4int i = 0; i < np1; ++i) 2761 for (G4int i = 0; i < np1; ++i)
2762 { 2762 {
2763 cosa = std::cos(sthe + i*a); 2763 cosa = std::cos(sthe + i*a);
2764 sina = std::sin(sthe + i*a); 2764 sina = std::sin(sthe + i*a);
2765 zz[i] = cz*cosa; 2765 zz[i] = cz*cosa;
2766 rr[i] = cz*sina; 2766 rr[i] = cz*sina;
2767 } 2767 }
2768 zz[np1 + 0] = zCut2; 2768 zz[np1 + 0] = zCut2;
2769 rr[np1 + 0] = 0.; 2769 rr[np1 + 0] = 0.;
2770 zz[np1 + 1] = zCut1; 2770 zz[np1 + 1] = zCut1;
2771 rr[np1 + 1] = 0.; 2771 rr[np1 + 1] = 0.;
2772 2772
2773 // R O T A T E P O L Y L I N E S 2773 // R O T A T E P O L Y L I N E S
2774 2774
2775 RotateAroundZ(0, 0., twopi, np1, np2, zz, r 2775 RotateAroundZ(0, 0., twopi, np1, np2, zz, rr, -1, -1);
2776 SetReferences(); 2776 SetReferences();
2777 2777
2778 delete [] zz; 2778 delete [] zz;
2779 delete [] rr; 2779 delete [] rr;
2780 2780
2781 // rescale x and y vertex coordinates 2781 // rescale x and y vertex coordinates
2782 G4double kx = ax/cz; 2782 G4double kx = ax/cz;
2783 G4double ky = by/cz; 2783 G4double ky = by/cz;
2784 G4Point3D* p = pV; 2784 G4Point3D* p = pV;
2785 for (G4int i = 0; i < nvert; ++i, ++p) 2785 for (G4int i = 0; i < nvert; ++i, ++p)
2786 { 2786 {
2787 p->setX(p->x()*kx); 2787 p->setX(p->x()*kx);
2788 p->setY(p->y()*ky); 2788 p->setY(p->y()*ky);
2789 } 2789 }
2790 } 2790 }
2791 2791
2792 HepPolyhedronEllipsoid::~HepPolyhedronEllipso 2792 HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() = default;
2793 2793
2794 HepPolyhedronEllipticalCone::HepPolyhedronEll 2794 HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone(G4double ax,
2795 2795 G4double ay,
2796 2796 G4double h,
2797 2797 G4double zTopCut)
2798 /******************************************** 2798 /***********************************************************************
2799 * 2799 * *
2800 * Name: HepPolyhedronEllipticalCone 2800 * Name: HepPolyhedronEllipticalCone Date: 8.9.2005 *
2801 * Author: D.Anninos 2801 * Author: D.Anninos Revised: 9.9.2005 *
2802 * 2802 * *
2803 * Function: Constructor for EllipticalCone 2803 * Function: Constructor for EllipticalCone *
2804 * 2804 * *
2805 * Input: ax, ay - X & Y semi axes at z = 2805 * Input: ax, ay - X & Y semi axes at z = 0 *
2806 * h - height of full cone 2806 * h - height of full cone *
2807 * zTopCut - Top Cut in Z Axis 2807 * zTopCut - Top Cut in Z Axis *
2808 * 2808 * *
2809 ******************************************** 2809 ***********************************************************************/
2810 { 2810 {
2811 // C H E C K I N P U T P A R A M E T 2811 // C H E C K I N P U T P A R A M E T E R S
2812 2812
2813 G4int k = 0; 2813 G4int k = 0;
2814 if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) 2814 if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
2815 2815
2816 if (k != 0) { 2816 if (k != 0) {
2817 std::cerr << "HepPolyhedronCone: error in 2817 std::cerr << "HepPolyhedronCone: error in input parameters";
2818 std::cerr << std::endl; 2818 std::cerr << std::endl;
2819 return; 2819 return;
2820 } 2820 }
2821 2821
2822 // P R E P A R E T W O P O L Y L I N 2822 // P R E P A R E T W O P O L Y L I N E S
2823 2823
2824 zTopCut = (h >= zTopCut ? zTopCut : h); 2824 zTopCut = (h >= zTopCut ? zTopCut : h);
2825 2825
2826 G4double *zz, *rr; 2826 G4double *zz, *rr;
2827 zz = new G4double[4]; 2827 zz = new G4double[4];
2828 rr = new G4double[4]; 2828 rr = new G4double[4];
2829 zz[0] = zTopCut; 2829 zz[0] = zTopCut;
2830 zz[1] = -zTopCut; 2830 zz[1] = -zTopCut;
2831 zz[2] = zTopCut; 2831 zz[2] = zTopCut;
2832 zz[3] = -zTopCut; 2832 zz[3] = -zTopCut;
2833 rr[0] = (h-zTopCut); 2833 rr[0] = (h-zTopCut);
2834 rr[1] = (h+zTopCut); 2834 rr[1] = (h+zTopCut);
2835 rr[2] = 0.; 2835 rr[2] = 0.;
2836 rr[3] = 0.; 2836 rr[3] = 0.;
2837 2837
2838 // R O T A T E P O L Y L I N E S 2838 // R O T A T E P O L Y L I N E S
2839 2839
2840 RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, - 2840 RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, -1, -1);
2841 SetReferences(); 2841 SetReferences();
2842 2842
2843 delete [] zz; 2843 delete [] zz;
2844 delete [] rr; 2844 delete [] rr;
2845 2845
2846 // rescale x and y vertex coordinates 2846 // rescale x and y vertex coordinates
2847 { 2847 {
2848 G4Point3D * p= pV; 2848 G4Point3D * p= pV;
2849 for (G4int i=0; i<nvert; i++, p++) { 2849 for (G4int i=0; i<nvert; i++, p++) {
2850 p->setX( p->x() * ax ); 2850 p->setX( p->x() * ax );
2851 p->setY( p->y() * ay ); 2851 p->setY( p->y() * ay );
2852 } 2852 }
2853 } 2853 }
2854 } 2854 }
2855 2855
2856 HepPolyhedronEllipticalCone::~HepPolyhedronEl 2856 HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() = default;
2857 2857
2858 HepPolyhedronHyperbolicMirror::HepPolyhedronH 2858 HepPolyhedronHyperbolicMirror::HepPolyhedronHyperbolicMirror(G4double a,
2859 2859 G4double h,
2860 2860 G4double r)
2861 /******************************************** 2861 /***********************************************************************
2862 * 2862 * *
2863 * Name: HepPolyhedronHyperbolicMirror 2863 * Name: HepPolyhedronHyperbolicMirror Date: 22.02.2020 *
2864 * Author: E.Tcherniaev (E.Chernyaev) 2864 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2865 * 2865 * *
2866 * Function: Create polyhedron for Hyperbolic 2866 * Function: Create polyhedron for Hyperbolic mirror *
2867 * 2867 * *
2868 * Input: a - half-separation 2868 * Input: a - half-separation *
2869 * h - height 2869 * h - height *
2870 * r - radius 2870 * r - radius *
2871 * 2871 * *
2872 ******************************************** 2872 ***********************************************************************/
2873 { 2873 {
2874 G4double H = std::abs(h); 2874 G4double H = std::abs(h);
2875 G4double R = std::abs(r); 2875 G4double R = std::abs(r);
2876 G4double A = std::abs(a); 2876 G4double A = std::abs(a);
2877 G4double B = A*R/std::sqrt(2*A*H + H*H); 2877 G4double B = A*R/std::sqrt(2*A*H + H*H);
2878 2878
2879 // P R E P A R E T W O P O L Y L I N 2879 // P R E P A R E T W O P O L Y L I N E S
2880 2880
2881 G4int np1 = (A == 0.) ? 2 : std::max(3, Get 2881 G4int np1 = (A == 0.) ? 2 : std::max(3, GetNumberOfRotationSteps()/4) + 1;
2882 G4int np2 = 2; 2882 G4int np2 = 2;
2883 G4double maxAng = (A == 0.) ? 0. : std::aco 2883 G4double maxAng = (A == 0.) ? 0. : std::acosh(1. + H/A);
2884 G4double delAng = maxAng/(np1 - 1); 2884 G4double delAng = maxAng/(np1 - 1);
2885 2885
2886 auto zz = new G4double[np1 + np2]; 2886 auto zz = new G4double[np1 + np2];
2887 auto rr = new G4double[np1 + np2]; 2887 auto rr = new G4double[np1 + np2];
2888 2888
2889 // 1st polyline 2889 // 1st polyline
2890 zz[0] = H; 2890 zz[0] = H;
2891 rr[0] = R; 2891 rr[0] = R;
2892 for (G4int iz = 1; iz < np1 - 1; ++iz) 2892 for (G4int iz = 1; iz < np1 - 1; ++iz)
2893 { 2893 {
2894 G4double ang = maxAng - iz*delAng; 2894 G4double ang = maxAng - iz*delAng;
2895 zz[iz] = A*std::cosh(ang) - A; 2895 zz[iz] = A*std::cosh(ang) - A;
2896 rr[iz] = B*std::sinh(ang); 2896 rr[iz] = B*std::sinh(ang);
2897 } 2897 }
2898 zz[np1 - 1] = 0.; 2898 zz[np1 - 1] = 0.;
2899 rr[np1 - 1] = 0.; 2899 rr[np1 - 1] = 0.;
2900 2900
2901 // 2nd polyline 2901 // 2nd polyline
2902 zz[np1] = H; 2902 zz[np1] = H;
2903 rr[np1] = 0.; 2903 rr[np1] = 0.;
2904 zz[np1 + 1] = 0.; 2904 zz[np1 + 1] = 0.;
2905 rr[np1 + 1] = 0.; 2905 rr[np1 + 1] = 0.;
2906 2906
2907 // R O T A T E P O L Y L I N E S 2907 // R O T A T E P O L Y L I N E S
2908 2908
2909 G4double phi = 0.; 2909 G4double phi = 0.;
2910 G4double dphi = CLHEP::twopi; 2910 G4double dphi = CLHEP::twopi;
2911 RotateAroundZ(0, phi, dphi, np1, np2, zz, r 2911 RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
2912 SetReferences(); 2912 SetReferences();
2913 2913
2914 delete [] zz; 2914 delete [] zz;
2915 delete [] rr; 2915 delete [] rr;
2916 } 2916 }
2917 2917
2918 HepPolyhedronHyperbolicMirror::~HepPolyhedron 2918 HepPolyhedronHyperbolicMirror::~HepPolyhedronHyperbolicMirror() = default;
2919 2919
2920 HepPolyhedronTetMesh:: 2920 HepPolyhedronTetMesh::
2921 HepPolyhedronTetMesh(const std::vector<G4Thre 2921 HepPolyhedronTetMesh(const std::vector<G4ThreeVector>& tetrahedra)
2922 /******************************************** 2922 /***********************************************************************
2923 * 2923 * *
2924 * Name: HepPolyhedronTetMesh 2924 * Name: HepPolyhedronTetMesh Date: 26.03.2022 *
2925 * Author: E.Tcherniaev (E.Chernyaev) 2925 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2926 * 2926 * *
2927 * Function: Create polyhedron for tetrahedro 2927 * Function: Create polyhedron for tetrahedron mesh *
2928 * 2928 * *
2929 * Input: tetrahedra - array of tetrahedron v 2929 * Input: tetrahedra - array of tetrahedron vertices, four vertices *
2930 * per tetrahedron 2930 * per tetrahedron *
2931 * 2931 * *
2932 ******************************************** 2932 ***********************************************************************/
2933 { 2933 {
2934 // Check size of input vector 2934 // Check size of input vector
2935 G4int nnodes = (G4int)tetrahedra.size(); 2935 G4int nnodes = (G4int)tetrahedra.size();
2936 if (nnodes == 0) 2936 if (nnodes == 0)
2937 { 2937 {
2938 std::cerr 2938 std::cerr
2939 << "HepPolyhedronTetMesh: Empty tetrahe 2939 << "HepPolyhedronTetMesh: Empty tetrahedron mesh" << std::endl;
2940 return; 2940 return;
2941 } 2941 }
2942 G4int ntet = nnodes/4; 2942 G4int ntet = nnodes/4;
2943 if (nnodes != ntet*4) 2943 if (nnodes != ntet*4)
2944 { 2944 {
2945 std::cerr << "HepPolyhedronTetMesh: Numbe 2945 std::cerr << "HepPolyhedronTetMesh: Number of nodes = " << nnodes
2946 << " in tetrahedron mesh is NOT 2946 << " in tetrahedron mesh is NOT multiple of 4"
2947 << std::endl; 2947 << std::endl;
2948 return; 2948 return;
2949 } 2949 }
2950 2950
2951 // Find coincident vertices using hash tabl 2951 // Find coincident vertices using hash table techniques.
2952 // This could be done using std::unordered_ 2952 // This could be done using std::unordered_map, but the code
2953 // below runs faster. 2953 // below runs faster.
2954 std::vector<G4int> iheads(nnodes, -1); 2954 std::vector<G4int> iheads(nnodes, -1);
2955 std::vector<std::pair<G4int,G4int>> ipairs( 2955 std::vector<std::pair<G4int,G4int>> ipairs(nnodes,std::pair(-1,-1));
2956 for (G4int i = 0; i < nnodes; ++i) 2956 for (G4int i = 0; i < nnodes; ++i)
2957 { 2957 {
2958 // Generate hash key 2958 // Generate hash key
2959 G4ThreeVector point = tetrahedra[i]; 2959 G4ThreeVector point = tetrahedra[i];
2960 auto key = std::hash<G4double>()(point.x( 2960 auto key = std::hash<G4double>()(point.x());
2961 key ^= std::hash<G4double>()(point.y()); 2961 key ^= std::hash<G4double>()(point.y());
2962 key ^= std::hash<G4double>()(point.z()); 2962 key ^= std::hash<G4double>()(point.z());
2963 key %= nnodes; 2963 key %= nnodes;
2964 // Check head of the list 2964 // Check head of the list
2965 if (iheads[key] < 0) 2965 if (iheads[key] < 0)
2966 { 2966 {
2967 iheads[key] = i; 2967 iheads[key] = i;
2968 ipairs[i].first = i; 2968 ipairs[i].first = i;
2969 continue; 2969 continue;
2970 } 2970 }
2971 // Loop along the list 2971 // Loop along the list
2972 for (G4int icur = iheads[key], iprev = 0; 2972 for (G4int icur = iheads[key], iprev = 0;;)
2973 { 2973 {
2974 G4int icheck = ipairs[icur].first; 2974 G4int icheck = ipairs[icur].first;
2975 if (tetrahedra[icheck] == point) 2975 if (tetrahedra[icheck] == point)
2976 { 2976 {
2977 ipairs[i].first = icheck; // coincide 2977 ipairs[i].first = icheck; // coincident vertex
2978 break; 2978 break;
2979 } 2979 }
2980 iprev = icur; 2980 iprev = icur;
2981 icur = ipairs[icur].second; 2981 icur = ipairs[icur].second;
2982 // Append vertex to the list 2982 // Append vertex to the list
2983 if (icur < 0) 2983 if (icur < 0)
2984 { 2984 {
2985 ipairs[i].first = i; 2985 ipairs[i].first = i;
2986 ipairs[iprev].second = i; 2986 ipairs[iprev].second = i;
2987 break; 2987 break;
2988 } 2988 }
2989 } 2989 }
2990 } 2990 }
2991 2991
2992 // Create vector of original facets 2992 // Create vector of original facets
2993 struct facet 2993 struct facet
2994 { 2994 {
2995 G4int i1, i2, i3; 2995 G4int i1, i2, i3;
2996 facet() : i1(0), i2(0), i3(0) {}; 2996 facet() : i1(0), i2(0), i3(0) {};
2997 facet(G4int k1, G4int k2, G4int k3) : i1( 2997 facet(G4int k1, G4int k2, G4int k3) : i1(k1), i2(k2), i3(k3) {};
2998 }; 2998 };
2999 G4int nfacets = nnodes; 2999 G4int nfacets = nnodes;
3000 std::vector<facet> ifacets(nfacets); 3000 std::vector<facet> ifacets(nfacets);
3001 for (G4int i = 0; i < nfacets; i += 4) 3001 for (G4int i = 0; i < nfacets; i += 4)
3002 { 3002 {
3003 G4int i0 = ipairs[i + 0].first; 3003 G4int i0 = ipairs[i + 0].first;
3004 G4int i1 = ipairs[i + 1].first; 3004 G4int i1 = ipairs[i + 1].first;
3005 G4int i2 = ipairs[i + 2].first; 3005 G4int i2 = ipairs[i + 2].first;
3006 G4int i3 = ipairs[i + 3].first; 3006 G4int i3 = ipairs[i + 3].first;
3007 if (i0 > i1) std::swap(i0, i1); 3007 if (i0 > i1) std::swap(i0, i1);
3008 if (i0 > i2) std::swap(i0, i2); 3008 if (i0 > i2) std::swap(i0, i2);
3009 if (i0 > i3) std::swap(i0, i3); 3009 if (i0 > i3) std::swap(i0, i3);
3010 if (i1 > i2) std::swap(i1, i2); 3010 if (i1 > i2) std::swap(i1, i2);
3011 if (i1 > i3) std::swap(i1, i3); 3011 if (i1 > i3) std::swap(i1, i3);
3012 G4ThreeVector e1 = tetrahedra[i1] - tetra 3012 G4ThreeVector e1 = tetrahedra[i1] - tetrahedra[i0];
3013 G4ThreeVector e2 = tetrahedra[i2] - tetra 3013 G4ThreeVector e2 = tetrahedra[i2] - tetrahedra[i0];
3014 G4ThreeVector e3 = tetrahedra[i3] - tetra 3014 G4ThreeVector e3 = tetrahedra[i3] - tetrahedra[i0];
3015 G4double volume = (e1.cross(e2)).dot(e3); 3015 G4double volume = (e1.cross(e2)).dot(e3);
3016 if (volume > 0.) std::swap(i2, i3); 3016 if (volume > 0.) std::swap(i2, i3);
3017 ifacets[i + 0] = facet(i0, i1, i2); 3017 ifacets[i + 0] = facet(i0, i1, i2);
3018 ifacets[i + 1] = facet(i0, i2, i3); 3018 ifacets[i + 1] = facet(i0, i2, i3);
3019 ifacets[i + 2] = facet(i0, i3, i1); 3019 ifacets[i + 2] = facet(i0, i3, i1);
3020 ifacets[i + 3] = facet(i1, i3, i2); 3020 ifacets[i + 3] = facet(i1, i3, i2);
3021 } 3021 }
3022 3022
3023 // Find shared facets 3023 // Find shared facets
3024 std::fill(iheads.begin(), iheads.end(), -1) 3024 std::fill(iheads.begin(), iheads.end(), -1);
3025 std::fill(ipairs.begin(), ipairs.end(), std 3025 std::fill(ipairs.begin(), ipairs.end(), std::pair(-1,-1));
3026 for (G4int i = 0; i < nfacets; ++i) 3026 for (G4int i = 0; i < nfacets; ++i)
3027 { 3027 {
3028 // Check head of the list 3028 // Check head of the list
3029 G4int key = ifacets[i].i1; 3029 G4int key = ifacets[i].i1;
3030 if (iheads[key] < 0) 3030 if (iheads[key] < 0)
3031 { 3031 {
3032 iheads[key] = i; 3032 iheads[key] = i;
3033 ipairs[i].first = i; 3033 ipairs[i].first = i;
3034 continue; 3034 continue;
3035 } 3035 }
3036 // Loop along the list 3036 // Loop along the list
3037 G4int i2 = ifacets[i].i2, i3 = ifacets[i] 3037 G4int i2 = ifacets[i].i2, i3 = ifacets[i].i3;
3038 for (G4int icur = iheads[key], iprev = -1 3038 for (G4int icur = iheads[key], iprev = -1;;)
3039 { 3039 {
3040 G4int icheck = ipairs[icur].first; 3040 G4int icheck = ipairs[icur].first;
3041 if (ifacets[icheck].i2 == i3 && ifacets 3041 if (ifacets[icheck].i2 == i3 && ifacets[icheck].i3 == i2)
3042 { 3042 {
3043 if (iprev < 0) 3043 if (iprev < 0)
3044 { 3044 {
3045 iheads[key] = ipairs[icur].second; 3045 iheads[key] = ipairs[icur].second;
3046 } 3046 }
3047 else 3047 else
3048 { 3048 {
3049 ipairs[iprev].second = ipairs[icur] 3049 ipairs[iprev].second = ipairs[icur].second;
3050 } 3050 }
3051 ipairs[icur].first = -1; // shared fa 3051 ipairs[icur].first = -1; // shared facet
3052 ipairs[icur].second = -1; 3052 ipairs[icur].second = -1;
3053 break; 3053 break;
3054 } 3054 }
3055 iprev = icur; 3055 iprev = icur;
3056 icur = ipairs[icur].second; 3056 icur = ipairs[icur].second;
3057 // Append facet to the list 3057 // Append facet to the list
3058 if (icur < 0) 3058 if (icur < 0)
3059 { 3059 {
3060 ipairs[i].first = i; 3060 ipairs[i].first = i;
3061 ipairs[iprev].second = i; 3061 ipairs[iprev].second = i;
3062 break; 3062 break;
3063 } 3063 }
3064 } 3064 }
3065 } 3065 }
3066 3066
3067 // Count vertices and facets skipping share 3067 // Count vertices and facets skipping shared facets
3068 std::fill(iheads.begin(), iheads.end(), -1) 3068 std::fill(iheads.begin(), iheads.end(), -1);
3069 G4int nver = 0, nfac = 0; 3069 G4int nver = 0, nfac = 0;
3070 for (G4int i = 0; i < nfacets; ++i) 3070 for (G4int i = 0; i < nfacets; ++i)
3071 { 3071 {
3072 if (ipairs[i].first < 0) continue; 3072 if (ipairs[i].first < 0) continue;
3073 G4int i1 = ifacets[i].i1; 3073 G4int i1 = ifacets[i].i1;
3074 G4int i2 = ifacets[i].i2; 3074 G4int i2 = ifacets[i].i2;
3075 G4int i3 = ifacets[i].i3; 3075 G4int i3 = ifacets[i].i3;
3076 if (iheads[i1] < 0) iheads[i1] = nver++; 3076 if (iheads[i1] < 0) iheads[i1] = nver++;
3077 if (iheads[i2] < 0) iheads[i2] = nver++; 3077 if (iheads[i2] < 0) iheads[i2] = nver++;
3078 if (iheads[i3] < 0) iheads[i3] = nver++; 3078 if (iheads[i3] < 0) iheads[i3] = nver++;
3079 nfac++; 3079 nfac++;
3080 } 3080 }
3081 3081
3082 // Construct polyhedron 3082 // Construct polyhedron
3083 AllocateMemory(nver, nfac); 3083 AllocateMemory(nver, nfac);
3084 for (G4int i = 0; i < nnodes; ++i) 3084 for (G4int i = 0; i < nnodes; ++i)
3085 { 3085 {
3086 G4int k = iheads[i]; 3086 G4int k = iheads[i];
3087 if (k >= 0) SetVertex(k + 1, tetrahedra[i 3087 if (k >= 0) SetVertex(k + 1, tetrahedra[i]);
3088 } 3088 }
3089 for (G4int i = 0, k = 0; i < nfacets; ++i) 3089 for (G4int i = 0, k = 0; i < nfacets; ++i)
3090 { 3090 {
3091 if (ipairs[i].first < 0) continue; 3091 if (ipairs[i].first < 0) continue;
3092 G4int i1 = iheads[ifacets[i].i1] + 1; 3092 G4int i1 = iheads[ifacets[i].i1] + 1;
3093 G4int i2 = iheads[ifacets[i].i2] + 1; 3093 G4int i2 = iheads[ifacets[i].i2] + 1;
3094 G4int i3 = iheads[ifacets[i].i3] + 1; 3094 G4int i3 = iheads[ifacets[i].i3] + 1;
3095 SetFacet(++k, i1, i2, i3); 3095 SetFacet(++k, i1, i2, i3);
3096 } 3096 }
3097 SetReferences(); 3097 SetReferences();
3098 } 3098 }
3099 3099
3100 HepPolyhedronTetMesh::~HepPolyhedronTetMesh() 3100 HepPolyhedronTetMesh::~HepPolyhedronTetMesh() = default;
3101 3101
3102 HepPolyhedronBoxMesh:: 3102 HepPolyhedronBoxMesh::
3103 HepPolyhedronBoxMesh(G4double sizeX, G4double 3103 HepPolyhedronBoxMesh(G4double sizeX, G4double sizeY, G4double sizeZ,
3104 const std::vector<G4Thre 3104 const std::vector<G4ThreeVector>& positions)
3105 /******************************************** 3105 /***********************************************************************
3106 * 3106 * *
3107 * Name: HepPolyhedronBoxMesh 3107 * Name: HepPolyhedronBoxMesh Date: 07.04.2022 *
3108 * Author: E.Tcherniaev (E.Chernyaev) 3108 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
3109 * 3109 * *
3110 * Function: Create polyhedron for box mesh 3110 * Function: Create polyhedron for box mesh *
3111 * 3111 * *
3112 * Input: sizeX, sizeY, sizeZ - dimensions of 3112 * Input: sizeX, sizeY, sizeZ - dimensions of the mesh cell *
3113 * positions - vector of cell centres 3113 * positions - vector of cell centres *
3114 * 3114 * *
3115 ******************************************** 3115 ***********************************************************************/
3116 { 3116 {
3117 G4int nbox = (G4int)positions.size(); 3117 G4int nbox = (G4int)positions.size();
3118 if (nbox == 0) 3118 if (nbox == 0)
3119 { 3119 {
3120 std::cerr << "HepPolyhedronBoxMesh: Empty 3120 std::cerr << "HepPolyhedronBoxMesh: Empty box mesh" << std::endl;
3121 return; 3121 return;
3122 } 3122 }
3123 // compute inverse dimensions 3123 // compute inverse dimensions
3124 G4double invx = 1./sizeX, invy = 1./sizeY, 3124 G4double invx = 1./sizeX, invy = 1./sizeY, invz = 1./sizeZ;
3125 // find mesh bounding box 3125 // find mesh bounding box
3126 G4ThreeVector pmin = positions[0], pmax = p 3126 G4ThreeVector pmin = positions[0], pmax = positions[0];
3127 for (const auto& p: positions) 3127 for (const auto& p: positions)
3128 { 3128 {
3129 if (pmin.x() > p.x()) pmin.setX(p.x()); 3129 if (pmin.x() > p.x()) pmin.setX(p.x());
3130 if (pmin.y() > p.y()) pmin.setY(p.y()); 3130 if (pmin.y() > p.y()) pmin.setY(p.y());
3131 if (pmin.z() > p.z()) pmin.setZ(p.z()); 3131 if (pmin.z() > p.z()) pmin.setZ(p.z());
3132 if (pmax.x() < p.x()) pmax.setX(p.x()); 3132 if (pmax.x() < p.x()) pmax.setX(p.x());
3133 if (pmax.y() < p.y()) pmax.setY(p.y()); 3133 if (pmax.y() < p.y()) pmax.setY(p.y());
3134 if (pmax.z() < p.z()) pmax.setZ(p.z()); 3134 if (pmax.z() < p.z()) pmax.setZ(p.z());
3135 } 3135 }
3136 // find number of voxels 3136 // find number of voxels
3137 G4int nx = (pmax.x() - pmin.x())*invx + 1.5 3137 G4int nx = (pmax.x() - pmin.x())*invx + 1.5;
3138 G4int ny = (pmax.y() - pmin.y())*invy + 1.5 3138 G4int ny = (pmax.y() - pmin.y())*invy + 1.5;
3139 G4int nz = (pmax.z() - pmin.z())*invz + 1.5 3139 G4int nz = (pmax.z() - pmin.z())*invz + 1.5;
3140 // create structures for voxels and node in 3140 // create structures for voxels and node indices
3141 std::vector<char> voxels(nx*ny*nz, 0); 3141 std::vector<char> voxels(nx*ny*nz, 0);
3142 std::vector<G4int> indices((nx+1)*(ny+1)*(n 3142 std::vector<G4int> indices((nx+1)*(ny+1)*(nz+1), 0);
3143 // mark voxels listed in positions 3143 // mark voxels listed in positions
3144 G4int kx = ny*nz, ky = nz; 3144 G4int kx = ny*nz, ky = nz;
3145 for (const auto& p: positions) 3145 for (const auto& p: positions)
3146 { 3146 {
3147 G4int ix = (p.x() - pmin.x())*invx + 0.5; 3147 G4int ix = (p.x() - pmin.x())*invx + 0.5;
3148 G4int iy = (p.y() - pmin.y())*invy + 0.5; 3148 G4int iy = (p.y() - pmin.y())*invy + 0.5;
3149 G4int iz = (p.z() - pmin.z())*invz + 0.5; 3149 G4int iz = (p.z() - pmin.z())*invz + 0.5;
3150 G4int i = ix*kx + iy*ky + iz; 3150 G4int i = ix*kx + iy*ky + iz;
3151 voxels[i] = 1; 3151 voxels[i] = 1;
3152 } 3152 }
3153 // count number of vertices and facets 3153 // count number of vertices and facets
3154 // set indices 3154 // set indices
3155 G4int kvx = (ny + 1)*(nz + 1), kvy = nz + 1 3155 G4int kvx = (ny + 1)*(nz + 1), kvy = nz + 1;
3156 G4int nver = 0, nfac = 0; 3156 G4int nver = 0, nfac = 0;
3157 for (const auto& p: positions) 3157 for (const auto& p: positions)
3158 { 3158 {
3159 G4int ix = (p.x() - pmin.x())*invx + 0.5; 3159 G4int ix = (p.x() - pmin.x())*invx + 0.5;
3160 G4int iy = (p.y() - pmin.y())*invy + 0.5; 3160 G4int iy = (p.y() - pmin.y())*invy + 0.5;
3161 G4int iz = (p.z() - pmin.z())*invz + 0.5; 3161 G4int iz = (p.z() - pmin.z())*invz + 0.5;
3162 // 3162 //
3163 // 011 111 3163 // 011 111
3164 // +---–---+ 3164 // +---–---+
3165 // | 001 | 101 3165 // | 001 | 101
3166 // | +---–---+ 3166 // | +---–---+
3167 // | | | | 3167 // | | | |
3168 // +---|---+ | 3168 // +---|---+ |
3169 // 010 | 110 | 3169 // 010 | 110 |
3170 // +-------+ 3170 // +-------+
3171 // 000 100 3171 // 000 100
3172 // 3172 //
3173 G4int vcheck = 0; 3173 G4int vcheck = 0;
3174 // check (ix - 1) side 3174 // check (ix - 1) side
3175 vcheck = (ix == 0) ? 0 : voxels[(ix-1)*kx 3175 vcheck = (ix == 0) ? 0 : voxels[(ix-1)*kx + iy*ky + iz];
3176 if (vcheck == 0) 3176 if (vcheck == 0)
3177 { 3177 {
3178 nfac++; 3178 nfac++;
3179 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i 3179 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3180 G4int i2 = (ix+0)*kvx + (iy+0)*kvy + (i 3180 G4int i2 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3181 G4int i3 = (ix+0)*kvx + (iy+1)*kvy + (i 3181 G4int i3 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3182 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i 3182 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3183 if (indices[i1] == 0) indices[i1] = ++n 3183 if (indices[i1] == 0) indices[i1] = ++nver;
3184 if (indices[i2] == 0) indices[i2] = ++n 3184 if (indices[i2] == 0) indices[i2] = ++nver;
3185 if (indices[i3] == 0) indices[i3] = ++n 3185 if (indices[i3] == 0) indices[i3] = ++nver;
3186 if (indices[i4] == 0) indices[i4] = ++n 3186 if (indices[i4] == 0) indices[i4] = ++nver;
3187 } 3187 }
3188 // check (ix + 1) side 3188 // check (ix + 1) side
3189 vcheck = (ix == nx - 1) ? 0 : voxels[(ix+ 3189 vcheck = (ix == nx - 1) ? 0 : voxels[(ix+1)*kx + iy*ky + iz];
3190 if (vcheck == 0) 3190 if (vcheck == 0)
3191 { 3191 {
3192 nfac++; 3192 nfac++;
3193 G4int i1 = (ix+1)*kvx + (iy+1)*kvy + (i 3193 G4int i1 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3194 G4int i2 = (ix+1)*kvx + (iy+1)*kvy + (i 3194 G4int i2 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3195 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i 3195 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3196 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i 3196 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3197 if (indices[i1] == 0) indices[i1] = ++n 3197 if (indices[i1] == 0) indices[i1] = ++nver;
3198 if (indices[i2] == 0) indices[i2] = ++n 3198 if (indices[i2] == 0) indices[i2] = ++nver;
3199 if (indices[i3] == 0) indices[i3] = ++n 3199 if (indices[i3] == 0) indices[i3] = ++nver;
3200 if (indices[i4] == 0) indices[i4] = ++n 3200 if (indices[i4] == 0) indices[i4] = ++nver;
3201 } 3201 }
3202 // check (iy - 1) side 3202 // check (iy - 1) side
3203 vcheck = (iy == 0) ? 0 : voxels[ix*kx + ( 3203 vcheck = (iy == 0) ? 0 : voxels[ix*kx + (iy-1)*ky + iz];
3204 if (vcheck == 0) 3204 if (vcheck == 0)
3205 { 3205 {
3206 nfac++; 3206 nfac++;
3207 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i 3207 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3208 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i 3208 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3209 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i 3209 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3210 G4int i4 = (ix+0)*kvx + (iy+0)*kvy + (i 3210 G4int i4 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3211 if (indices[i1] == 0) indices[i1] = ++n 3211 if (indices[i1] == 0) indices[i1] = ++nver;
3212 if (indices[i2] == 0) indices[i2] = ++n 3212 if (indices[i2] == 0) indices[i2] = ++nver;
3213 if (indices[i3] == 0) indices[i3] = ++n 3213 if (indices[i3] == 0) indices[i3] = ++nver;
3214 if (indices[i4] == 0) indices[i4] = ++n 3214 if (indices[i4] == 0) indices[i4] = ++nver;
3215 } 3215 }
3216 // check (iy + 1) side 3216 // check (iy + 1) side
3217 vcheck = (iy == ny - 1) ? 0 : voxels[ix*k 3217 vcheck = (iy == ny - 1) ? 0 : voxels[ix*kx + (iy+1)*ky + iz];
3218 if (vcheck == 0) 3218 if (vcheck == 0)
3219 { 3219 {
3220 nfac++; 3220 nfac++;
3221 G4int i1 = (ix+0)*kvx + (iy+1)*kvy + (i 3221 G4int i1 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3222 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i 3222 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3223 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i 3223 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3224 G4int i4 = (ix+1)*kvx + (iy+1)*kvy + (i 3224 G4int i4 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3225 if (indices[i1] == 0) indices[i1] = ++n 3225 if (indices[i1] == 0) indices[i1] = ++nver;
3226 if (indices[i2] == 0) indices[i2] = ++n 3226 if (indices[i2] == 0) indices[i2] = ++nver;
3227 if (indices[i3] == 0) indices[i3] = ++n 3227 if (indices[i3] == 0) indices[i3] = ++nver;
3228 if (indices[i4] == 0) indices[i4] = ++n 3228 if (indices[i4] == 0) indices[i4] = ++nver;
3229 } 3229 }
3230 // check (iz - 1) side 3230 // check (iz - 1) side
3231 vcheck = (iz == 0) ? 0 : voxels[ix*kx + i 3231 vcheck = (iz == 0) ? 0 : voxels[ix*kx + iy*ky + (iz-1)];
3232 if (vcheck == 0) 3232 if (vcheck == 0)
3233 { 3233 {
3234 nfac++; 3234 nfac++;
3235 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i 3235 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3236 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i 3236 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3237 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i 3237 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3238 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i 3238 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3239 if (indices[i1] == 0) indices[i1] = ++n 3239 if (indices[i1] == 0) indices[i1] = ++nver;
3240 if (indices[i2] == 0) indices[i2] = ++n 3240 if (indices[i2] == 0) indices[i2] = ++nver;
3241 if (indices[i3] == 0) indices[i3] = ++n 3241 if (indices[i3] == 0) indices[i3] = ++nver;
3242 if (indices[i4] == 0) indices[i4] = ++n 3242 if (indices[i4] == 0) indices[i4] = ++nver;
3243 } 3243 }
3244 // check (iz + 1) side 3244 // check (iz + 1) side
3245 vcheck = (iz == nz - 1) ? 0 : voxels[ix*k 3245 vcheck = (iz == nz - 1) ? 0 : voxels[ix*kx + iy*ky + (iz+1)];
3246 if (vcheck == 0) 3246 if (vcheck == 0)
3247 { 3247 {
3248 nfac++; 3248 nfac++;
3249 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i 3249 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3250 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i 3250 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3251 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i 3251 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3252 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i 3252 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3253 if (indices[i1] == 0) indices[i1] = ++n 3253 if (indices[i1] == 0) indices[i1] = ++nver;
3254 if (indices[i2] == 0) indices[i2] = ++n 3254 if (indices[i2] == 0) indices[i2] = ++nver;
3255 if (indices[i3] == 0) indices[i3] = ++n 3255 if (indices[i3] == 0) indices[i3] = ++nver;
3256 if (indices[i4] == 0) indices[i4] = ++n 3256 if (indices[i4] == 0) indices[i4] = ++nver;
3257 } 3257 }
3258 } 3258 }
3259 // Construct polyhedron 3259 // Construct polyhedron
3260 AllocateMemory(nver, nfac); 3260 AllocateMemory(nver, nfac);
3261 G4ThreeVector p0(pmin.x() - 0.5*sizeX, pmin 3261 G4ThreeVector p0(pmin.x() - 0.5*sizeX, pmin.y() - 0.5*sizeY, pmin.z() - 0.5*sizeZ);
3262 for (G4int ix = 0; ix <= nx; ++ix) 3262 for (G4int ix = 0; ix <= nx; ++ix)
3263 { 3263 {
3264 for (G4int iy = 0; iy <= ny; ++iy) 3264 for (G4int iy = 0; iy <= ny; ++iy)
3265 { 3265 {
3266 for (G4int iz = 0; iz <= nz; ++iz) 3266 for (G4int iz = 0; iz <= nz; ++iz)
3267 { 3267 {
3268 G4int i = ix*kvx + iy*kvy + iz; 3268 G4int i = ix*kvx + iy*kvy + iz;
3269 if (indices[i] == 0) continue; 3269 if (indices[i] == 0) continue;
3270 SetVertex(indices[i], p0 + G4ThreeVector(ix 3270 SetVertex(indices[i], p0 + G4ThreeVector(ix*sizeX, iy*sizeY, iz*sizeZ));
3271 } 3271 }
3272 } 3272 }
3273 } 3273 }
3274 nfac = 0; 3274 nfac = 0;
3275 for (const auto& p: positions) 3275 for (const auto& p: positions)
3276 { 3276 {
3277 G4int ix = (p.x() - pmin.x())*invx + 0.5; 3277 G4int ix = (p.x() - pmin.x())*invx + 0.5;
3278 G4int iy = (p.y() - pmin.y())*invy + 0.5; 3278 G4int iy = (p.y() - pmin.y())*invy + 0.5;
3279 G4int iz = (p.z() - pmin.z())*invz + 0.5; 3279 G4int iz = (p.z() - pmin.z())*invz + 0.5;
3280 G4int vcheck = 0; 3280 G4int vcheck = 0;
3281 // check (ix - 1) side 3281 // check (ix - 1) side
3282 vcheck = (ix == 0) ? 0 : voxels[(ix-1)*kx 3282 vcheck = (ix == 0) ? 0 : voxels[(ix-1)*kx + iy*ky + iz];
3283 if (vcheck == 0) 3283 if (vcheck == 0)
3284 { 3284 {
3285 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i 3285 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3286 G4int i2 = (ix+0)*kvx + (iy+0)*kvy + (i 3286 G4int i2 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3287 G4int i3 = (ix+0)*kvx + (iy+1)*kvy + (i 3287 G4int i3 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3288 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i 3288 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3289 SetFacet(++nfac, indices[i1], indices[i 3289 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3290 } 3290 }
3291 // check (ix + 1) side 3291 // check (ix + 1) side
3292 vcheck = (ix == nx - 1) ? 0 : voxels[(ix+ 3292 vcheck = (ix == nx - 1) ? 0 : voxels[(ix+1)*kx + iy*ky + iz];
3293 if (vcheck == 0) 3293 if (vcheck == 0)
3294 { 3294 {
3295 G4int i1 = (ix+1)*kvx + (iy+1)*kvy + (i 3295 G4int i1 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3296 G4int i2 = (ix+1)*kvx + (iy+1)*kvy + (i 3296 G4int i2 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3297 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i 3297 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3298 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i 3298 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3299 SetFacet(++nfac, indices[i1], indices[i 3299 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3300 3300
3301 } 3301 }
3302 // check (iy - 1) side 3302 // check (iy - 1) side
3303 vcheck = (iy == 0) ? 0 : voxels[ix*kx + ( 3303 vcheck = (iy == 0) ? 0 : voxels[ix*kx + (iy-1)*ky + iz];
3304 if (vcheck == 0) 3304 if (vcheck == 0)
3305 { 3305 {
3306 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i 3306 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3307 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i 3307 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3308 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i 3308 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3309 G4int i4 = (ix+0)*kvx + (iy+0)*kvy + (i 3309 G4int i4 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3310 SetFacet(++nfac, indices[i1], indices[i 3310 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3311 } 3311 }
3312 // check (iy + 1) side 3312 // check (iy + 1) side
3313 vcheck = (iy == ny - 1) ? 0 : voxels[ix*k 3313 vcheck = (iy == ny - 1) ? 0 : voxels[ix*kx + (iy+1)*ky + iz];
3314 if (vcheck == 0) 3314 if (vcheck == 0)
3315 { 3315 {
3316 G4int i1 = (ix+0)*kvx + (iy+1)*kvy + (i 3316 G4int i1 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3317 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i 3317 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3318 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i 3318 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3319 G4int i4 = (ix+1)*kvx + (iy+1)*kvy + (i 3319 G4int i4 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3320 SetFacet(++nfac, indices[i1], indices[i 3320 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3321 } 3321 }
3322 // check (iz - 1) side 3322 // check (iz - 1) side
3323 vcheck = (iz == 0) ? 0 : voxels[ix*kx + i 3323 vcheck = (iz == 0) ? 0 : voxels[ix*kx + iy*ky + (iz-1)];
3324 if (vcheck == 0) 3324 if (vcheck == 0)
3325 { 3325 {
3326 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i 3326 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3327 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i 3327 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3328 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i 3328 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3329 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i 3329 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3330 SetFacet(++nfac, indices[i1], indices[i 3330 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3331 } 3331 }
3332 // check (iz + 1) side 3332 // check (iz + 1) side
3333 vcheck = (iz == nz - 1) ? 0 : voxels[ix*k 3333 vcheck = (iz == nz - 1) ? 0 : voxels[ix*kx + iy*ky + (iz+1)];
3334 if (vcheck == 0) 3334 if (vcheck == 0)
3335 { 3335 {
3336 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i 3336 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3337 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i 3337 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3338 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i 3338 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3339 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i 3339 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3340 SetFacet(++nfac, indices[i1], indices[i 3340 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3341 } 3341 }
3342 } 3342 }
3343 SetReferences(); 3343 SetReferences();
3344 } 3344 }
3345 3345
3346 HepPolyhedronBoxMesh::~HepPolyhedronBoxMesh() 3346 HepPolyhedronBoxMesh::~HepPolyhedronBoxMesh() = default;
3347 3347
3348 G4ThreadLocal 3348 G4ThreadLocal
3349 G4int HepPolyhedron::fNumberOfRotationSteps = 3349 G4int HepPolyhedron::fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
3350 /******************************************** 3350 /***********************************************************************
3351 * 3351 * *
3352 * Name: HepPolyhedron::fNumberOfRotationStep 3352 * Name: HepPolyhedron::fNumberOfRotationSteps Date: 24.06.97 *
3353 * Author: J.Allison (Manchester University) 3353 * Author: J.Allison (Manchester University) Revised: *
3354 * 3354 * *
3355 * Function: Number of steps for whole circle 3355 * Function: Number of steps for whole circle *
3356 * 3356 * *
3357 ******************************************** 3357 ***********************************************************************/
3358 3358
3359 #include "BooleanProcessor.src" 3359 #include "BooleanProcessor.src"
3360 3360
3361 HepPolyhedron HepPolyhedron::add(const HepPol 3361 HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const
3362 /******************************************** 3362 /***********************************************************************
3363 * 3363 * *
3364 * Name: HepPolyhedron::add 3364 * Name: HepPolyhedron::add Date: 19.03.00 *
3365 * Author: E.Chernyaev 3365 * Author: E.Chernyaev Revised: *
3366 * 3366 * *
3367 * Function: Boolean "union" of two polyhedra 3367 * Function: Boolean "union" of two polyhedra *
3368 * 3368 * *
3369 ******************************************** 3369 ***********************************************************************/
3370 { 3370 {
3371 G4int ierr; 3371 G4int ierr;
3372 BooleanProcessor processor; 3372 BooleanProcessor processor;
3373 return processor.execute(OP_UNION, *this, p 3373 return processor.execute(OP_UNION, *this, p,ierr);
3374 } 3374 }
3375 3375
3376 HepPolyhedron HepPolyhedron::intersect(const 3376 HepPolyhedron HepPolyhedron::intersect(const HepPolyhedron & p) const
3377 /******************************************** 3377 /***********************************************************************
3378 * 3378 * *
3379 * Name: HepPolyhedron::intersect 3379 * Name: HepPolyhedron::intersect Date: 19.03.00 *
3380 * Author: E.Chernyaev 3380 * Author: E.Chernyaev Revised: *
3381 * 3381 * *
3382 * Function: Boolean "intersection" of two po 3382 * Function: Boolean "intersection" of two polyhedra *
3383 * 3383 * *
3384 ******************************************** 3384 ***********************************************************************/
3385 { 3385 {
3386 G4int ierr; 3386 G4int ierr;
3387 BooleanProcessor processor; 3387 BooleanProcessor processor;
3388 return processor.execute(OP_INTERSECTION, * 3388 return processor.execute(OP_INTERSECTION, *this, p,ierr);
3389 } 3389 }
3390 3390
3391 HepPolyhedron HepPolyhedron::subtract(const H 3391 HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const
3392 /******************************************** 3392 /***********************************************************************
3393 * 3393 * *
3394 * Name: HepPolyhedron::add 3394 * Name: HepPolyhedron::add Date: 19.03.00 *
3395 * Author: E.Chernyaev 3395 * Author: E.Chernyaev Revised: *
3396 * 3396 * *
3397 * Function: Boolean "subtraction" of "p" fro 3397 * Function: Boolean "subtraction" of "p" from "this" *
3398 * 3398 * *
3399 ******************************************** 3399 ***********************************************************************/
3400 { 3400 {
3401 G4int ierr; 3401 G4int ierr;
3402 BooleanProcessor processor; 3402 BooleanProcessor processor;
3403 return processor.execute(OP_SUBTRACTION, *t 3403 return processor.execute(OP_SUBTRACTION, *this, p,ierr);
3404 } 3404 }
3405 3405
3406 //NOTE : include the code of HepPolyhedronPro 3406 //NOTE : include the code of HepPolyhedronProcessor here
3407 // since there is no BooleanProcessor.h 3407 // since there is no BooleanProcessor.h
3408 3408
3409 #undef INTERSECTION 3409 #undef INTERSECTION
3410 3410
3411 #include "HepPolyhedronProcessor.src" 3411 #include "HepPolyhedronProcessor.src"
3412 3412