| 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 // Implementation of G4UGenericPolycone wrappe <<
27 // 26 //
28 // 30.10.13 G.Cosmo, CERN << 27 // $Id:$
>> 28 //
>> 29 //
>> 30 // Implementation of G4UGenericPolycone wrapper class
29 // ------------------------------------------- 31 // --------------------------------------------------------------------
30 32
31 #include "G4GenericPolycone.hh" 33 #include "G4GenericPolycone.hh"
32 #include "G4UGenericPolycone.hh" 34 #include "G4UGenericPolycone.hh"
33 <<
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G <<
35 <<
36 #include "G4GeomTools.hh" <<
37 #include "G4AffineTransform.hh" <<
38 #include "G4VPVParameterisation.hh" <<
39 #include "G4BoundingEnvelope.hh" <<
40 <<
41 #include "G4Polyhedron.hh" 35 #include "G4Polyhedron.hh"
42 36
43 using namespace CLHEP; << 37 using CLHEP::twopi;
44 38
45 ////////////////////////////////////////////// 39 ////////////////////////////////////////////////////////////////////////
46 // 40 //
47 // Constructor (generic parameters) 41 // Constructor (generic parameters)
48 // 42 //
49 G4UGenericPolycone::G4UGenericPolycone(const G << 43 G4UGenericPolycone::G4UGenericPolycone(const G4String& name,
50 G 44 G4double phiStart,
51 G 45 G4double phiTotal,
52 G 46 G4int numRZ,
53 const G 47 const G4double r[],
54 const G 48 const G4double z[] )
55 : Base_t(name, phiStart, phiTotal, numRZ, r, << 49 : G4USolid(name, new UGenericPolycone(name, phiStart, phiTotal, numRZ, r, z))
56 { << 50 {
57 wrStart = phiStart; while (wrStart < 0) wrSt <<
58 wrDelta = phiTotal; <<
59 if (wrDelta <= 0 || wrDelta >= twopi*(1-DBL_ <<
60 { <<
61 wrStart = 0; <<
62 wrDelta = twopi; <<
63 } <<
64 rzcorners.resize(0); <<
65 for (G4int i=0; i<numRZ; ++i) <<
66 { <<
67 rzcorners.emplace_back(r[i],z[i]); <<
68 } <<
69 std::vector<G4int> iout; <<
70 G4GeomTools::RemoveRedundantVertices(rzcorne <<
71 } 51 }
72 52
73 53
74 ////////////////////////////////////////////// 54 ////////////////////////////////////////////////////////////////////////
75 // 55 //
76 // Fake default constructor - sets only member 56 // Fake default constructor - sets only member data and allocates memory
77 // for usage restri 57 // for usage restricted to object persistency.
78 // 58 //
79 G4UGenericPolycone::G4UGenericPolycone(__void_ 59 G4UGenericPolycone::G4UGenericPolycone(__void__& a)
80 : Base_t(a) << 60 : G4USolid(a)
81 { 61 {
82 } 62 }
83 63
84 64
85 ////////////////////////////////////////////// 65 //////////////////////////////////////////////////////////////////////////
86 // 66 //
87 // Destructor 67 // Destructor
88 // 68 //
89 G4UGenericPolycone::~G4UGenericPolycone() = de << 69 G4UGenericPolycone::~G4UGenericPolycone()
>> 70 {
>> 71 }
90 72
91 73
92 ////////////////////////////////////////////// 74 //////////////////////////////////////////////////////////////////////////
93 // 75 //
94 // Copy constructor 76 // Copy constructor
95 // 77 //
96 G4UGenericPolycone::G4UGenericPolycone(const G << 78 G4UGenericPolycone::G4UGenericPolycone(const G4UGenericPolycone &source)
97 : Base_t(source) << 79 : G4USolid(source)
98 { 80 {
99 wrStart = source.wrStart; <<
100 wrDelta = source.wrDelta; <<
101 rzcorners = source.rzcorners; <<
102 } 81 }
103 82
104 83
105 ////////////////////////////////////////////// 84 //////////////////////////////////////////////////////////////////////////
106 // 85 //
107 // Assignment operator 86 // Assignment operator
108 // 87 //
109 G4UGenericPolycone& 88 G4UGenericPolycone&
110 G4UGenericPolycone::operator=(const G4UGeneric << 89 G4UGenericPolycone::operator=(const G4UGenericPolycone &source)
111 { 90 {
112 if (this == &source) return *this; 91 if (this == &source) return *this;
113 << 92
114 Base_t::operator=( source ); << 93 G4USolid::operator=( source );
115 wrStart = source.wrStart; << 94
116 wrDelta = source.wrDelta; <<
117 rzcorners = source.rzcorners; <<
118 <<
119 return *this; 95 return *this;
120 } 96 }
121 97
122 G4double G4UGenericPolycone::GetStartPhi() con << 98 G4Polyhedron* G4UGenericPolycone::CreatePolyhedron() const
123 { <<
124 return wrStart; <<
125 } <<
126 G4double G4UGenericPolycone::GetEndPhi() const <<
127 { <<
128 return (wrStart + wrDelta); <<
129 } <<
130 G4double G4UGenericPolycone::GetSinStartPhi() <<
131 { <<
132 if (IsOpen()) return 0.; <<
133 G4double phi = GetStartPhi(); <<
134 return std::sin(phi); <<
135 } <<
136 G4double G4UGenericPolycone::GetCosStartPhi() <<
137 { <<
138 if (IsOpen()) return 1.; <<
139 G4double phi = GetStartPhi(); <<
140 return std::cos(phi); <<
141 } <<
142 G4double G4UGenericPolycone::GetSinEndPhi() co <<
143 { <<
144 if (IsOpen()) return 0.; <<
145 G4double phi = GetEndPhi(); <<
146 return std::sin(phi); <<
147 } <<
148 G4double G4UGenericPolycone::GetCosEndPhi() co <<
149 { <<
150 if (IsOpen()) return 1.; <<
151 G4double phi = GetEndPhi(); <<
152 return std::cos(phi); <<
153 } <<
154 G4bool G4UGenericPolycone::IsOpen() const <<
155 { <<
156 return (wrDelta < twopi); <<
157 } <<
158 G4int G4UGenericPolycone::GetNumRZCorner() con <<
159 { <<
160 return rzcorners.size(); <<
161 } <<
162 G4PolyconeSideRZ G4UGenericPolycone::GetCorner <<
163 { <<
164 G4TwoVector rz = rzcorners.at(index); <<
165 G4PolyconeSideRZ psiderz = { rz.x(), rz.y() <<
166 <<
167 return psiderz; <<
168 } <<
169 <<
170 ////////////////////////////////////////////// <<
171 // <<
172 // Make a clone of the object <<
173 <<
174 G4VSolid* G4UGenericPolycone::Clone() const <<
175 { 99 {
176 return new G4UGenericPolycone(*this); <<
177 } <<
178 100
179 ////////////////////////////////////////////// <<
180 // <<
181 // Get bounding box <<
182 101
183 void << 102 // The following code prepares for:
184 G4UGenericPolycone::BoundingLimits(G4ThreeVect << 103 // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
185 G4ThreeVect << 104 // const double xyz[][3],
186 { << 105 // const int faces_vec[][4])
187 G4double rmin = kInfinity, rmax = -kInfinity << 106 // Here is an extract from the header file HepPolyhedron.h:
188 G4double zmin = kInfinity, zmax = -kInfinity << 107 /**
189 << 108 * Creates user defined polyhedron.
190 for (G4int i=0; i<GetNumRZCorner(); ++i) << 109 * This function allows to the user to define arbitrary polyhedron.
191 { << 110 * The faces of the polyhedron should be either triangles or planar
192 G4PolyconeSideRZ corner = GetCorner(i); << 111 * quadrilateral. Nodes of a face are defined by indexes pointing to
193 if (corner.r < rmin) rmin = corner.r; << 112 * the elements in the xyz array. Numeration of the elements in the
194 if (corner.r > rmax) rmax = corner.r; << 113 * array starts from 1 (like in fortran). The indexes can be positive
195 if (corner.z < zmin) zmin = corner.z; << 114 * or negative. Negative sign means that the corresponding edge is
196 if (corner.z > zmax) zmax = corner.z; << 115 * invisible. The normal of the face should be directed to exterior
197 } << 116 * of the polyhedron.
198 << 117 *
199 if (IsOpen()) << 118 * @param Nnodes number of nodes
200 { << 119 * @param Nfaces number of faces
201 G4TwoVector vmin,vmax; << 120 * @param xyz nodes
202 G4GeomTools::DiskExtent(rmin,rmax, << 121 * @param faces_vec faces (quadrilaterals or triangles)
203 GetSinStartPhi(),G << 122 * @return status of the operation - is non-zero in case of problem
204 GetSinEndPhi(),Get << 123 */
205 vmin,vmax); << 124 const G4int numSide =
206 pMin.set(vmin.x(),vmin.y(),zmin); << 125 G4int(G4Polyhedron::GetNumberOfRotationSteps()
207 pMax.set(vmax.x(),vmax.y(),zmax); << 126 * (GetEndPhi() - GetStartPhi()) / twopi) + 1;
208 } << 127 G4int nNodes;
209 else << 128 G4int nFaces;
210 { << 129 typedef G4double double3[3];
211 pMin.set(-rmax,-rmax, zmin); << 130 double3* xyz;
212 pMax.set( rmax, rmax, zmax); << 131 typedef G4int int4[4];
213 } << 132 int4* faces_vec;
214 << 133 if (IsOpen())
215 // Check correctness of the bounding box << 134 {
216 // << 135 // Triangulate open ends. Simple ear-chopping algorithm...
217 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 136 // I'm not sure how robust this algorithm is (J.Allison).
218 { << 137 //
219 std::ostringstream message; << 138 std::vector<G4bool> chopped(GetNumRZCorner(), false);
220 message << "Bad bounding box (min >= max) << 139 std::vector<G4int*> triQuads;
221 << GetName() << " !" << 140 G4int remaining = GetNumRZCorner();
222 << "\npMin = " << pMin << 141 G4int iStarter = 0;
223 << "\npMax = " << pMax; << 142 while (remaining >= 3)
224 G4Exception("G4UGenericPolycone::BoundingL << 143 {
225 JustWarning, message); << 144 // Find unchopped corners...
226 StreamInfo(G4cout); << 145 //
227 } << 146 G4int A = -1, B = -1, C = -1;
228 } << 147 G4int iStepper = iStarter;
>> 148 do
>> 149 {
>> 150 if (A < 0) { A = iStepper; }
>> 151 else if (B < 0) { B = iStepper; }
>> 152 else if (C < 0) { C = iStepper; }
>> 153 do
>> 154 {
>> 155 if (++iStepper >= GetNumRZCorner()) { iStepper = 0; }
>> 156 }
>> 157 while (chopped[iStepper]);
>> 158 }
>> 159 while (C < 0 && iStepper != iStarter);
>> 160
>> 161 // Check triangle at B is pointing outward (an "ear").
>> 162 // Sign of z cross product determines...
>> 163 //
>> 164 G4double BAr = GetCorner(A).r - GetCorner(B).r;
>> 165 G4double BAz = GetCorner(A).z - GetCorner(B).z;
>> 166 G4double BCr = GetCorner(C).r - GetCorner(B).r;
>> 167 G4double BCz = GetCorner(C).z - GetCorner(B).z;
>> 168 if (BAr * BCz - BAz * BCr < kCarTolerance)
>> 169 {
>> 170 G4int* tq = new G4int[3];
>> 171 tq[0] = A + 1;
>> 172 tq[1] = B + 1;
>> 173 tq[2] = C + 1;
>> 174 triQuads.push_back(tq);
>> 175 chopped[B] = true;
>> 176 --remaining;
>> 177 }
>> 178 else
>> 179 {
>> 180 do
>> 181 {
>> 182 if (++iStarter >= GetNumRZCorner()) { iStarter = 0; }
>> 183 }
>> 184 while (chopped[iStarter]);
>> 185 }
>> 186 }
>> 187 // Transfer to faces...
>> 188 //
>> 189 nNodes = (numSide + 1) * GetNumRZCorner();
>> 190 nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size();
>> 191 faces_vec = new int4[nFaces];
>> 192 G4int iface = 0;
>> 193 G4int addition = GetNumRZCorner() * numSide;
>> 194 G4int d = GetNumRZCorner() - 1;
>> 195 for (G4int iEnd = 0; iEnd < 2; ++iEnd)
>> 196 {
>> 197 for (size_t i = 0; i < triQuads.size(); ++i)
>> 198 {
>> 199 // Negative for soft/auxiliary/normally invisible edges...
>> 200 //
>> 201 G4int a, b, c;
>> 202 if (iEnd == 0)
>> 203 {
>> 204 a = triQuads[i][0];
>> 205 b = triQuads[i][1];
>> 206 c = triQuads[i][2];
>> 207 }
>> 208 else
>> 209 {
>> 210 a = triQuads[i][0] + addition;
>> 211 b = triQuads[i][2] + addition;
>> 212 c = triQuads[i][1] + addition;
>> 213 }
>> 214 G4int ab = std::abs(b - a);
>> 215 G4int bc = std::abs(c - b);
>> 216 G4int ca = std::abs(a - c);
>> 217 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
>> 218 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
>> 219 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
>> 220 faces_vec[iface][3] = 0;
>> 221 ++iface;
>> 222 }
>> 223 }
229 224
230 ////////////////////////////////////////////// << 225 // Continue with sides...
231 // <<
232 // Calculate extent under transform and specif <<
233 226
234 G4bool << 227 xyz = new double3[nNodes];
235 G4UGenericPolycone::CalculateExtent(const EAxi << 228 const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
236 const G4Vo << 229 G4double phi = GetStartPhi();
237 const G4Af << 230 G4int ixyz = 0;
238 G4do << 231 for (G4int iSide = 0; iSide < numSide; ++iSide)
239 { << 232 {
240 G4ThreeVector bmin, bmax; << 233 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
241 G4bool exist; << 234 {
>> 235 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
>> 236 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
>> 237 xyz[ixyz][2] = GetCorner(iCorner).z;
>> 238 if (iSide == 0) // startPhi
>> 239 {
>> 240 if (iCorner < GetNumRZCorner() - 1)
>> 241 {
>> 242 faces_vec[iface][0] = ixyz + 1;
>> 243 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 244 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
>> 245 faces_vec[iface][3] = ixyz + 2;
>> 246 }
>> 247 else
>> 248 {
>> 249 faces_vec[iface][0] = ixyz + 1;
>> 250 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 251 faces_vec[iface][2] = ixyz + 2;
>> 252 faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
>> 253 }
>> 254 }
>> 255 else if (iSide == numSide - 1) // endPhi
>> 256 {
>> 257 if (iCorner < GetNumRZCorner() - 1)
>> 258 {
>> 259 faces_vec[iface][0] = ixyz + 1;
>> 260 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
>> 261 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
>> 262 faces_vec[iface][3] = -(ixyz + 2);
>> 263 }
>> 264 else
>> 265 {
>> 266 faces_vec[iface][0] = ixyz + 1;
>> 267 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
>> 268 faces_vec[iface][2] = ixyz + 2;
>> 269 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
>> 270 }
>> 271 }
>> 272 else
>> 273 {
>> 274 if (iCorner < GetNumRZCorner() - 1)
>> 275 {
>> 276 faces_vec[iface][0] = ixyz + 1;
>> 277 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 278 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
>> 279 faces_vec[iface][3] = -(ixyz + 2);
>> 280 }
>> 281 else
>> 282 {
>> 283 faces_vec[iface][0] = ixyz + 1;
>> 284 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 285 faces_vec[iface][2] = ixyz + 2;
>> 286 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
>> 287 }
>> 288 }
>> 289 ++iface;
>> 290 ++ixyz;
>> 291 }
>> 292 phi += dPhi;
>> 293 }
242 294
243 // Check bounding box (bbox) << 295 // Last corners...
244 // <<
245 BoundingLimits(bmin,bmax); <<
246 G4BoundingEnvelope bbox(bmin,bmax); <<
247 #ifdef G4BBOX_EXTENT <<
248 return bbox.CalculateExtent(pAxis,pVoxelLimi <<
249 #endif <<
250 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox <<
251 { <<
252 return exist = pMin < pMax; <<
253 } <<
254 <<
255 // To find the extent, RZ contour of the pol <<
256 // in triangles. The extent is calculated as <<
257 // all sub-polycones formed by rotation of t <<
258 // <<
259 G4TwoVectorList contourRZ; <<
260 G4TwoVectorList triangles; <<
261 G4double eminlim = pVoxelLimit.GetMinExtent( <<
262 G4double emaxlim = pVoxelLimit.GetMaxExtent( <<
263 <<
264 // get RZ contour, ensure anticlockwise orde <<
265 for (G4int i=0; i<GetNumRZCorner(); ++i) <<
266 { <<
267 G4PolyconeSideRZ corner = GetCorner(i); <<
268 contourRZ.emplace_back(corner.r,corner.z); <<
269 } <<
270 G4double area = G4GeomTools::PolygonArea(con <<
271 if (area < 0.) std::reverse(contourRZ.begin( <<
272 <<
273 // triangulate RZ countour <<
274 if (!G4GeomTools::TriangulatePolygon(contour <<
275 { <<
276 std::ostringstream message; <<
277 message << "Triangulation of RZ contour ha <<
278 << GetName() << " !" <<
279 << "\nExtent has been calculated u <<
280 G4Exception("G4UGenericPolycone::Calculate <<
281 "GeomMgt1002", JustWarning, me <<
282 return bbox.CalculateExtent(pAxis,pVoxelLi <<
283 } <<
284 <<
285 // set trigonometric values <<
286 const G4int NSTEPS = 24; // numbe <<
287 G4double astep = twopi/NSTEPS; // max a <<
288 <<
289 G4double sphi = GetStartPhi(); <<
290 G4double ephi = GetEndPhi(); <<
291 G4double dphi = IsOpen() ? ephi-sphi : two <<
292 G4int ksteps = (dphi <= astep) ? 1 : (G4i <<
293 G4double ang = dphi/ksteps; <<
294 <<
295 G4double sinHalf = std::sin(0.5*ang); <<
296 G4double cosHalf = std::cos(0.5*ang); <<
297 G4double sinStep = 2.*sinHalf*cosHalf; <<
298 G4double cosStep = 1. - 2.*sinHalf*sinHalf; <<
299 <<
300 G4double sinStart = GetSinStartPhi(); <<
301 G4double cosStart = GetCosStartPhi(); <<
302 G4double sinEnd = GetSinEndPhi(); <<
303 G4double cosEnd = GetCosEndPhi(); <<
304 <<
305 // define vectors and arrays <<
306 std::vector<const G4ThreeVectorList *> polyg <<
307 polygons.resize(ksteps+2); <<
308 G4ThreeVectorList pols[NSTEPS+2]; <<
309 for (G4int k=0; k<ksteps+2; ++k) pols[k].res <<
310 for (G4int k=0; k<ksteps+2; ++k) polygons[k] <<
311 G4double r0[6],z0[6]; // contour with origin <<
312 G4double r1[6]; // shifted radii of ex <<
313 <<
314 // main loop along triangles <<
315 pMin = kInfinity; <<
316 pMax =-kInfinity; <<
317 G4int ntria = triangles.size()/3; <<
318 for (G4int i=0; i<ntria; ++i) <<
319 { <<
320 G4int i3 = i*3; <<
321 for (G4int k=0; k<3; ++k) <<
322 { <<
323 G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3; <<
324 G4int k2 = k*2; <<
325 // set contour with original edges of tr <<
326 r0[k2+0] = triangles[e0].x(); z0[k2+0] = <<
327 r0[k2+1] = triangles[e1].x(); z0[k2+1] = <<
328 // set shifted radii <<
329 r1[k2+0] = r0[k2+0]; <<
330 r1[k2+1] = r0[k2+1]; <<
331 if (z0[k2+1] - z0[k2+0] <= 0) continue; <<
332 r1[k2+0] /= cosHalf; <<
333 r1[k2+1] /= cosHalf; <<
334 } <<
335 296
336 // rotate countour, set sequence of 6-side << 297 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
337 G4double sinCur = sinStart*cosHalf + cosSt << 298 {
338 G4double cosCur = cosStart*cosHalf - sinSt << 299 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
339 for (G4int j=0; j<6; ++j) << 300 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
340 { << 301 xyz[ixyz][2] = GetCorner(iCorner).z;
341 pols[0][j].set(r0[j]*cosStart,r0[j]*sinS << 302 ++ixyz;
>> 303 }
342 } 304 }
343 for (G4int k=1; k<ksteps+1; ++k) << 305 else // !phiIsOpen - i.e., a complete 360 degrees.
344 { 306 {
345 for (G4int j=0; j<6; ++j) << 307 nNodes = numSide * GetNumRZCorner();
>> 308 nFaces = numSide * GetNumRZCorner();;
>> 309 xyz = new double3[nNodes];
>> 310 faces_vec = new int4[nFaces];
>> 311 const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
>> 312 G4double phi = GetStartPhi();
>> 313 G4int ixyz = 0, iface = 0;
>> 314 for (G4int iSide = 0; iSide < numSide; ++iSide)
346 { 315 {
347 pols[k][j].set(r1[j]*cosCur,r1[j]*sinC << 316 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
>> 317 {
>> 318 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
>> 319 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
>> 320 xyz[ixyz][2] = GetCorner(iCorner).z;
>> 321
>> 322 if (iSide < numSide - 1)
>> 323 {
>> 324 if (iCorner < GetNumRZCorner() - 1)
>> 325 {
>> 326 faces_vec[iface][0] = ixyz + 1;
>> 327 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 328 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
>> 329 faces_vec[iface][3] = -(ixyz + 2);
>> 330 }
>> 331 else
>> 332 {
>> 333 faces_vec[iface][0] = ixyz + 1;
>> 334 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
>> 335 faces_vec[iface][2] = ixyz + 2;
>> 336 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
>> 337 }
>> 338 }
>> 339 else // Last side joins ends...
>> 340 {
>> 341 if (iCorner < GetNumRZCorner() - 1)
>> 342 {
>> 343 faces_vec[iface][0] = ixyz + 1;
>> 344 faces_vec[iface][1] = -(ixyz + GetNumRZCorner() - nFaces + 1);
>> 345 faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2;
>> 346 faces_vec[iface][3] = -(ixyz + 2);
>> 347 }
>> 348 else
>> 349 {
>> 350 faces_vec[iface][0] = ixyz + 1;
>> 351 faces_vec[iface][1] = -(ixyz - nFaces + GetNumRZCorner() + 1);
>> 352 faces_vec[iface][2] = ixyz - nFaces + 2;
>> 353 faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
>> 354 }
>> 355 }
>> 356 ++ixyz;
>> 357 ++iface;
>> 358 }
>> 359 phi += dPhi;
348 } 360 }
349 G4double sinTmp = sinCur; <<
350 sinCur = sinCur*cosStep + cosCur*sinStep <<
351 cosCur = cosCur*cosStep - sinTmp*sinStep <<
352 } 361 }
353 for (G4int j=0; j<6; ++j) << 362 G4Polyhedron* polyhedron = new G4Polyhedron;
>> 363 G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
>> 364 delete [] faces_vec;
>> 365 delete [] xyz;
>> 366 if (problem)
354 { 367 {
355 pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j] << 368 std::ostringstream message;
>> 369 message << "Problem creating G4Polyhedron for: " << GetName();
>> 370 G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002",
>> 371 JustWarning, message);
>> 372 delete polyhedron;
>> 373 return 0;
>> 374 }
>> 375 else
>> 376 {
>> 377 return polyhedron;
356 } 378 }
357 <<
358 // set sub-envelope and adjust extent <<
359 G4double emin,emax; <<
360 G4BoundingEnvelope benv(polygons); <<
361 if (!benv.CalculateExtent(pAxis,pVoxelLimi <<
362 if (emin < pMin) pMin = emin; <<
363 if (emax > pMax) pMax = emax; <<
364 if (eminlim > pMin && emaxlim < pMax) retu <<
365 } <<
366 return (pMin < pMax); <<
367 } <<
368 <<
369 ////////////////////////////////////////////// <<
370 // <<
371 // CreatePolyhedron <<
372 <<
373 G4Polyhedron* G4UGenericPolycone::CreatePolyhe <<
374 { <<
375 return new G4PolyhedronPcon(wrStart, wrDelta <<
376 } 379 }
377 380
378 #endif // G4GEOM_USE_USOLIDS <<
379 381