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