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