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