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