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