| 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 G4UPolyhedra wrapper clas 26 // Implementation of G4UPolyhedra wrapper class
27 // 27 //
28 // 31.10.13 G.Cosmo, CERN 28 // 31.10.13 G.Cosmo, CERN
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4Polyhedra.hh" 31 #include "G4Polyhedra.hh"
32 #include "G4UPolyhedra.hh" 32 #include "G4UPolyhedra.hh"
33 33
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G 34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
35 35
36 #include "G4GeomTools.hh" 36 #include "G4GeomTools.hh"
37 #include "G4GeometryTolerance.hh" 37 #include "G4GeometryTolerance.hh"
38 #include "G4AffineTransform.hh" 38 #include "G4AffineTransform.hh"
39 #include "G4VPVParameterisation.hh" 39 #include "G4VPVParameterisation.hh"
40 #include "G4BoundingEnvelope.hh" 40 #include "G4BoundingEnvelope.hh"
41 41
42 using namespace CLHEP; 42 using namespace CLHEP;
43 43
44 ////////////////////////////////////////////// 44 ////////////////////////////////////////////////////////////////////////
45 // 45 //
46 // Constructor (GEANT3 style parameters) 46 // Constructor (GEANT3 style parameters)
47 // 47 //
48 // GEANT3 PGON radii are specified in the dist 48 // GEANT3 PGON radii are specified in the distance to the norm of each face.
49 // << 49 //
50 G4UPolyhedra::G4UPolyhedra(const G4String& nam << 50 G4UPolyhedra::G4UPolyhedra(const G4String& name,
51 G4double phiS 51 G4double phiStart,
52 G4double phiT 52 G4double phiTotal,
53 G4int numSide << 53 G4int numSide,
54 G4int numZPla 54 G4int numZPlanes,
55 const G4double zPla 55 const G4double zPlane[],
56 const G4double rInn 56 const G4double rInner[],
57 const G4double rOut 57 const G4double rOuter[] )
58 : Base_t(name, phiStart, phiTotal, numSide, 58 : Base_t(name, phiStart, phiTotal, numSide,
59 numZPlanes, zPlane, rInner, rOuter) 59 numZPlanes, zPlane, rInner, rOuter)
60 { 60 {
61 fGenericPgon = false; 61 fGenericPgon = false;
62 SetOriginalParameters(); 62 SetOriginalParameters();
63 wrStart = phiStart; 63 wrStart = phiStart;
64 while (wrStart < 0) 64 while (wrStart < 0)
65 { 65 {
66 wrStart += twopi; 66 wrStart += twopi;
67 } 67 }
68 wrDelta = phiTotal; 68 wrDelta = phiTotal;
69 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL 69 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL_EPSILON))
70 { 70 {
71 wrDelta = twopi; 71 wrDelta = twopi;
72 } 72 }
73 wrNumSide = numSide; 73 wrNumSide = numSide;
74 G4double convertRad = 1./std::cos(0.5*wrDelt 74 G4double convertRad = 1./std::cos(0.5*wrDelta/wrNumSide);
75 rzcorners.resize(0); 75 rzcorners.resize(0);
76 for (G4int i=0; i<numZPlanes; ++i) 76 for (G4int i=0; i<numZPlanes; ++i)
77 { 77 {
78 G4double z = zPlane[i]; 78 G4double z = zPlane[i];
79 G4double r = rOuter[i]*convertRad; 79 G4double r = rOuter[i]*convertRad;
80 rzcorners.emplace_back(r,z); << 80 rzcorners.push_back(G4TwoVector(r,z));
81 } 81 }
82 for (G4int i=numZPlanes-1; i>=0; --i) 82 for (G4int i=numZPlanes-1; i>=0; --i)
83 { 83 {
84 G4double z = zPlane[i]; 84 G4double z = zPlane[i];
85 G4double r = rInner[i]*convertRad; 85 G4double r = rInner[i]*convertRad;
86 rzcorners.emplace_back(r,z); << 86 rzcorners.push_back(G4TwoVector(r,z));
87 } 87 }
88 std::vector<G4int> iout; 88 std::vector<G4int> iout;
89 G4GeomTools::RemoveRedundantVertices(rzcorne 89 G4GeomTools::RemoveRedundantVertices(rzcorners,iout,2*kCarTolerance);
90 } 90 }
91 91
92 92
93 ////////////////////////////////////////////// 93 ////////////////////////////////////////////////////////////////////////
94 // 94 //
95 // Constructor (generic parameters) 95 // Constructor (generic parameters)
96 // 96 //
97 G4UPolyhedra::G4UPolyhedra(const G4String& nam << 97 G4UPolyhedra::G4UPolyhedra(const G4String& name,
98 G4double phiS 98 G4double phiStart,
99 G4double phiT 99 G4double phiTotal,
100 G4int numS << 100 G4int numSide,
101 G4int numR 101 G4int numRZ,
102 const G4double r[], 102 const G4double r[],
103 const G4double z[] 103 const G4double z[] )
104 : Base_t(name, phiStart, phiTotal, numSide, 104 : Base_t(name, phiStart, phiTotal, numSide, numRZ, r, z)
105 { 105 {
106 fGenericPgon = true; 106 fGenericPgon = true;
107 SetOriginalParameters(); 107 SetOriginalParameters();
108 wrStart = phiStart; 108 wrStart = phiStart;
109 while (wrStart < 0.) 109 while (wrStart < 0.)
110 { 110 {
111 wrStart += twopi; 111 wrStart += twopi;
112 } 112 }
113 wrDelta = phiTotal; 113 wrDelta = phiTotal;
114 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL 114 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL_EPSILON))
115 { 115 {
116 wrDelta = twopi; 116 wrDelta = twopi;
117 } 117 }
118 wrNumSide = numSide; 118 wrNumSide = numSide;
119 rzcorners.resize(0); 119 rzcorners.resize(0);
120 for (G4int i=0; i<numRZ; ++i) 120 for (G4int i=0; i<numRZ; ++i)
121 { 121 {
122 rzcorners.emplace_back(r[i],z[i]); << 122 rzcorners.push_back(G4TwoVector(r[i],z[i]));
123 } 123 }
124 std::vector<G4int> iout; 124 std::vector<G4int> iout;
125 G4GeomTools::RemoveRedundantVertices(rzcorne 125 G4GeomTools::RemoveRedundantVertices(rzcorners,iout,2*kCarTolerance);
126 } 126 }
127 127
128 128
129 ////////////////////////////////////////////// 129 ////////////////////////////////////////////////////////////////////////
130 // 130 //
131 // Fake default constructor - sets only member 131 // Fake default constructor - sets only member data and allocates memory
132 // for usage restri 132 // for usage restricted to object persistency.
133 // 133 //
134 G4UPolyhedra::G4UPolyhedra( __void__& a ) 134 G4UPolyhedra::G4UPolyhedra( __void__& a )
135 : Base_t(a) 135 : Base_t(a)
136 { 136 {
137 } 137 }
138 138
139 139
140 ////////////////////////////////////////////// 140 ////////////////////////////////////////////////////////////////////////
141 // 141 //
142 // Destructor 142 // Destructor
143 // 143 //
144 G4UPolyhedra::~G4UPolyhedra() = default; << 144 G4UPolyhedra::~G4UPolyhedra()
>> 145 {
>> 146 }
145 147
146 148
147 ////////////////////////////////////////////// 149 ////////////////////////////////////////////////////////////////////////
148 // 150 //
149 // Copy constructor 151 // Copy constructor
150 // 152 //
151 G4UPolyhedra::G4UPolyhedra( const G4UPolyhedra 153 G4UPolyhedra::G4UPolyhedra( const G4UPolyhedra& source )
152 : Base_t( source ) 154 : Base_t( source )
153 { 155 {
154 fGenericPgon = source.fGenericPgon; 156 fGenericPgon = source.fGenericPgon;
155 fOriginalParameters = source.fOriginalParame 157 fOriginalParameters = source.fOriginalParameters;
156 wrStart = source.wrStart; 158 wrStart = source.wrStart;
157 wrDelta = source.wrDelta; 159 wrDelta = source.wrDelta;
158 wrNumSide = source.wrNumSide; 160 wrNumSide = source.wrNumSide;
159 rzcorners = source.rzcorners; 161 rzcorners = source.rzcorners;
160 } 162 }
161 163
162 164
163 ////////////////////////////////////////////// 165 ////////////////////////////////////////////////////////////////////////
164 // 166 //
165 // Assignment operator 167 // Assignment operator
166 // 168 //
167 G4UPolyhedra& G4UPolyhedra::operator=( const G 169 G4UPolyhedra& G4UPolyhedra::operator=( const G4UPolyhedra& source )
168 { 170 {
169 if (this == &source) return *this; 171 if (this == &source) return *this;
170 172
171 Base_t::operator=( source ); 173 Base_t::operator=( source );
172 fGenericPgon = source.fGenericPgon; 174 fGenericPgon = source.fGenericPgon;
173 fOriginalParameters = source.fOriginalParame 175 fOriginalParameters = source.fOriginalParameters;
174 wrStart = source.wrStart; 176 wrStart = source.wrStart;
175 wrDelta = source.wrDelta; 177 wrDelta = source.wrDelta;
176 wrNumSide = source.wrNumSide; 178 wrNumSide = source.wrNumSide;
177 rzcorners = source.rzcorners; 179 rzcorners = source.rzcorners;
178 180
179 return *this; 181 return *this;
180 } 182 }
181 183
182 184
183 ////////////////////////////////////////////// 185 ////////////////////////////////////////////////////////////////////////
184 // 186 //
185 // Accessors & modifiers 187 // Accessors & modifiers
186 // 188 //
187 G4int G4UPolyhedra::GetNumSide() const 189 G4int G4UPolyhedra::GetNumSide() const
188 { 190 {
189 return wrNumSide; 191 return wrNumSide;
190 } 192 }
191 G4double G4UPolyhedra::GetStartPhi() const 193 G4double G4UPolyhedra::GetStartPhi() const
192 { 194 {
193 return wrStart; 195 return wrStart;
194 } 196 }
195 G4double G4UPolyhedra::GetEndPhi() const 197 G4double G4UPolyhedra::GetEndPhi() const
196 { 198 {
197 return (wrStart + wrDelta); 199 return (wrStart + wrDelta);
198 } 200 }
199 G4double G4UPolyhedra::GetSinStartPhi() const 201 G4double G4UPolyhedra::GetSinStartPhi() const
200 { 202 {
201 G4double phi = GetStartPhi(); 203 G4double phi = GetStartPhi();
202 return std::sin(phi); 204 return std::sin(phi);
203 } 205 }
204 G4double G4UPolyhedra::GetCosStartPhi() const 206 G4double G4UPolyhedra::GetCosStartPhi() const
205 { 207 {
206 G4double phi = GetStartPhi(); 208 G4double phi = GetStartPhi();
207 return std::cos(phi); 209 return std::cos(phi);
208 } 210 }
209 G4double G4UPolyhedra::GetSinEndPhi() const 211 G4double G4UPolyhedra::GetSinEndPhi() const
210 { 212 {
211 G4double phi = GetEndPhi(); 213 G4double phi = GetEndPhi();
212 return std::sin(phi); 214 return std::sin(phi);
213 } 215 }
214 G4double G4UPolyhedra::GetCosEndPhi() const 216 G4double G4UPolyhedra::GetCosEndPhi() const
215 { 217 {
216 G4double phi = GetEndPhi(); 218 G4double phi = GetEndPhi();
217 return std::cos(phi); 219 return std::cos(phi);
218 } 220 }
219 G4bool G4UPolyhedra::IsOpen() const 221 G4bool G4UPolyhedra::IsOpen() const
220 { 222 {
221 return (wrDelta < twopi); 223 return (wrDelta < twopi);
222 } 224 }
223 G4bool G4UPolyhedra::IsGeneric() const 225 G4bool G4UPolyhedra::IsGeneric() const
224 { 226 {
225 return fGenericPgon; 227 return fGenericPgon;
226 } 228 }
227 G4int G4UPolyhedra::GetNumRZCorner() const 229 G4int G4UPolyhedra::GetNumRZCorner() const
228 { 230 {
229 return rzcorners.size(); 231 return rzcorners.size();
230 } 232 }
231 G4PolyhedraSideRZ G4UPolyhedra::GetCorner(G4in 233 G4PolyhedraSideRZ G4UPolyhedra::GetCorner(G4int index) const
232 { 234 {
233 G4TwoVector rz = rzcorners.at(index); 235 G4TwoVector rz = rzcorners.at(index);
234 G4PolyhedraSideRZ psiderz = { rz.x(), rz.y() 236 G4PolyhedraSideRZ psiderz = { rz.x(), rz.y() };
235 237
236 return psiderz; 238 return psiderz;
237 } 239 }
238 G4PolyhedraHistorical* G4UPolyhedra::GetOrigin 240 G4PolyhedraHistorical* G4UPolyhedra::GetOriginalParameters() const
239 { 241 {
240 return new G4PolyhedraHistorical(fOriginalPa 242 return new G4PolyhedraHistorical(fOriginalParameters);
241 } 243 }
242 void G4UPolyhedra::SetOriginalParameters() 244 void G4UPolyhedra::SetOriginalParameters()
243 { 245 {
244 G4double startPhi = GetPhiStart(); 246 G4double startPhi = GetPhiStart();
245 G4double deltaPhi = GetPhiDelta(); 247 G4double deltaPhi = GetPhiDelta();
246 G4int numPlanes = GetZSegmentCount() + 1; 248 G4int numPlanes = GetZSegmentCount() + 1;
247 G4int numSides = GetSideCount(); 249 G4int numSides = GetSideCount();
248 250
249 fOriginalParameters.Start_angle = startPhi 251 fOriginalParameters.Start_angle = startPhi;
250 fOriginalParameters.Opening_angle = deltaPhi 252 fOriginalParameters.Opening_angle = deltaPhi;
251 fOriginalParameters.Num_z_planes = numPlane 253 fOriginalParameters.Num_z_planes = numPlanes;
252 fOriginalParameters.numSide = numSides 254 fOriginalParameters.numSide = numSides;
253 255
254 delete [] fOriginalParameters.Z_values; 256 delete [] fOriginalParameters.Z_values;
255 delete [] fOriginalParameters.Rmin; 257 delete [] fOriginalParameters.Rmin;
256 delete [] fOriginalParameters.Rmax; 258 delete [] fOriginalParameters.Rmax;
257 fOriginalParameters.Z_values = new G4double[ 259 fOriginalParameters.Z_values = new G4double[numPlanes];
258 fOriginalParameters.Rmin = new G4double[numP 260 fOriginalParameters.Rmin = new G4double[numPlanes];
259 fOriginalParameters.Rmax = new G4double[numP 261 fOriginalParameters.Rmax = new G4double[numPlanes];
260 262
261 G4double convertRad = fGenericPgon 263 G4double convertRad = fGenericPgon
262 ? 1.0 : std::cos(0.5*del 264 ? 1.0 : std::cos(0.5*deltaPhi/numSides);
263 for (G4int i=0; i<numPlanes; ++i) 265 for (G4int i=0; i<numPlanes; ++i)
264 { 266 {
265 fOriginalParameters.Z_values[i] = GetZPlan 267 fOriginalParameters.Z_values[i] = GetZPlanes()[i];
266 fOriginalParameters.Rmax[i] = GetRMax( 268 fOriginalParameters.Rmax[i] = GetRMax()[i]/convertRad;
267 fOriginalParameters.Rmin[i] = GetRMin( 269 fOriginalParameters.Rmin[i] = GetRMin()[i]/convertRad;
268 } 270 }
269 } 271 }
270 void G4UPolyhedra::SetOriginalParameters(G4Pol 272 void G4UPolyhedra::SetOriginalParameters(G4PolyhedraHistorical* pars)
271 { 273 {
272 fOriginalParameters = *pars; 274 fOriginalParameters = *pars;
273 fRebuildPolyhedron = true; 275 fRebuildPolyhedron = true;
274 Reset(); 276 Reset();
275 } 277 }
276 278
277 G4bool G4UPolyhedra::Reset() 279 G4bool G4UPolyhedra::Reset()
278 { 280 {
279 if (fGenericPgon) 281 if (fGenericPgon)
280 { 282 {
281 std::ostringstream message; 283 std::ostringstream message;
282 message << "Solid " << GetName() << " buil 284 message << "Solid " << GetName() << " built using generic construct."
283 << G4endl << "Not applicable to th 285 << G4endl << "Not applicable to the generic construct !";
284 G4Exception("G4UPolyhedra::Reset()", "Geom 286 G4Exception("G4UPolyhedra::Reset()", "GeomSolids1001",
285 JustWarning, message, "Paramet 287 JustWarning, message, "Parameters NOT reset.");
286 return true; // error code set 288 return true; // error code set
287 } 289 }
288 290
289 // 291 //
290 // Rebuild polyhedra based on original param 292 // Rebuild polyhedra based on original parameters
291 // 293 //
292 wrStart = fOriginalParameters.Start_angle; 294 wrStart = fOriginalParameters.Start_angle;
293 while (wrStart < 0.) 295 while (wrStart < 0.)
294 { 296 {
295 wrStart += twopi; 297 wrStart += twopi;
296 } 298 }
297 wrDelta = fOriginalParameters.Opening_angle; 299 wrDelta = fOriginalParameters.Opening_angle;
298 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL 300 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL_EPSILON))
299 { 301 {
300 wrDelta = twopi; 302 wrDelta = twopi;
301 } 303 }
302 wrNumSide = fOriginalParameters.numSide; 304 wrNumSide = fOriginalParameters.numSide;
303 rzcorners.resize(0); 305 rzcorners.resize(0);
304 for (G4int i=0; i<fOriginalParameters.Num_z_ 306 for (G4int i=0; i<fOriginalParameters.Num_z_planes; ++i)
305 { 307 {
306 G4double z = fOriginalParameters.Z_values[ 308 G4double z = fOriginalParameters.Z_values[i];
307 G4double r = fOriginalParameters.Rmax[i]; 309 G4double r = fOriginalParameters.Rmax[i];
308 rzcorners.emplace_back(r,z); << 310 rzcorners.push_back(G4TwoVector(r,z));
309 } 311 }
310 for (G4int i=fOriginalParameters.Num_z_plane 312 for (G4int i=fOriginalParameters.Num_z_planes-1; i>=0; --i)
311 { 313 {
312 G4double z = fOriginalParameters.Z_values[ 314 G4double z = fOriginalParameters.Z_values[i];
313 G4double r = fOriginalParameters.Rmin[i]; 315 G4double r = fOriginalParameters.Rmin[i];
314 rzcorners.emplace_back(r,z); << 316 rzcorners.push_back(G4TwoVector(r,z));
315 } 317 }
316 std::vector<G4int> iout; 318 std::vector<G4int> iout;
317 G4GeomTools::RemoveRedundantVertices(rzcorne 319 G4GeomTools::RemoveRedundantVertices(rzcorners,iout,2*kCarTolerance);
318 320
319 return false; // error code unset 321 return false; // error code unset
320 } 322 }
321 323
322 324
323 ////////////////////////////////////////////// 325 ////////////////////////////////////////////////////////////////////////
324 // 326 //
325 // Dispatch to parameterisation for replicatio 327 // Dispatch to parameterisation for replication mechanism dimension
326 // computation & modification. 328 // computation & modification.
327 // 329 //
328 void G4UPolyhedra::ComputeDimensions(G4VPVPara 330 void G4UPolyhedra::ComputeDimensions(G4VPVParameterisation* p,
329 const G4i 331 const G4int n,
330 const G4V 332 const G4VPhysicalVolume* pRep)
331 { 333 {
332 p->ComputeDimensions(*(G4Polyhedra*)this,n,p 334 p->ComputeDimensions(*(G4Polyhedra*)this,n,pRep);
333 } 335 }
334 336
335 337
336 ////////////////////////////////////////////// 338 //////////////////////////////////////////////////////////////////////////
337 // 339 //
338 // Make a clone of the object 340 // Make a clone of the object
339 341
340 G4VSolid* G4UPolyhedra::Clone() const 342 G4VSolid* G4UPolyhedra::Clone() const
341 { 343 {
342 return new G4UPolyhedra(*this); 344 return new G4UPolyhedra(*this);
343 } 345 }
344 346
345 347
346 ////////////////////////////////////////////// 348 //////////////////////////////////////////////////////////////////////////
347 // 349 //
348 // Get bounding box 350 // Get bounding box
349 351
350 void G4UPolyhedra::BoundingLimits(G4ThreeVecto 352 void G4UPolyhedra::BoundingLimits(G4ThreeVector& pMin,
351 G4ThreeVecto 353 G4ThreeVector& pMax) const
352 { 354 {
353 static G4bool checkBBox = true; 355 static G4bool checkBBox = true;
354 static G4bool checkPhi = true; 356 static G4bool checkPhi = true;
355 357
356 G4double rmin = kInfinity, rmax = -kInfinity 358 G4double rmin = kInfinity, rmax = -kInfinity;
357 G4double zmin = kInfinity, zmax = -kInfinity 359 G4double zmin = kInfinity, zmax = -kInfinity;
358 for (G4int i=0; i<GetNumRZCorner(); ++i) 360 for (G4int i=0; i<GetNumRZCorner(); ++i)
359 { 361 {
360 G4PolyhedraSideRZ corner = GetCorner(i); 362 G4PolyhedraSideRZ corner = GetCorner(i);
361 if (corner.r < rmin) rmin = corner.r; 363 if (corner.r < rmin) rmin = corner.r;
362 if (corner.r > rmax) rmax = corner.r; 364 if (corner.r > rmax) rmax = corner.r;
363 if (corner.z < zmin) zmin = corner.z; 365 if (corner.z < zmin) zmin = corner.z;
364 if (corner.z > zmax) zmax = corner.z; 366 if (corner.z > zmax) zmax = corner.z;
365 } 367 }
366 368
367 G4double sphi = GetStartPhi(); 369 G4double sphi = GetStartPhi();
368 G4double ephi = GetEndPhi(); 370 G4double ephi = GetEndPhi();
369 G4double dphi = IsOpen() ? ephi-sphi : tw 371 G4double dphi = IsOpen() ? ephi-sphi : twopi;
370 G4int ksteps = GetNumSide(); 372 G4int ksteps = GetNumSide();
371 G4double astep = dphi/ksteps; 373 G4double astep = dphi/ksteps;
372 G4double sinStep = std::sin(astep); 374 G4double sinStep = std::sin(astep);
373 G4double cosStep = std::cos(astep); 375 G4double cosStep = std::cos(astep);
374 376
375 G4double sinCur = GetSinStartPhi(); 377 G4double sinCur = GetSinStartPhi();
376 G4double cosCur = GetCosStartPhi(); 378 G4double cosCur = GetCosStartPhi();
377 if (!IsOpen()) rmin = 0.; 379 if (!IsOpen()) rmin = 0.;
378 G4double xmin = rmin*cosCur, xmax = xmin; 380 G4double xmin = rmin*cosCur, xmax = xmin;
379 G4double ymin = rmin*sinCur, ymax = ymin; 381 G4double ymin = rmin*sinCur, ymax = ymin;
380 for (G4int k=0; k<ksteps+1; ++k) 382 for (G4int k=0; k<ksteps+1; ++k)
381 { 383 {
382 G4double x = rmax*cosCur; 384 G4double x = rmax*cosCur;
383 if (x < xmin) xmin = x; 385 if (x < xmin) xmin = x;
384 if (x > xmax) xmax = x; 386 if (x > xmax) xmax = x;
385 G4double y = rmax*sinCur; 387 G4double y = rmax*sinCur;
386 if (y < ymin) ymin = y; 388 if (y < ymin) ymin = y;
387 if (y > ymax) ymax = y; 389 if (y > ymax) ymax = y;
388 if (rmin > 0.) 390 if (rmin > 0.)
389 { 391 {
390 G4double xx = rmin*cosCur; 392 G4double xx = rmin*cosCur;
391 if (xx < xmin) xmin = xx; 393 if (xx < xmin) xmin = xx;
392 if (xx > xmax) xmax = xx; 394 if (xx > xmax) xmax = xx;
393 G4double yy = rmin*sinCur; 395 G4double yy = rmin*sinCur;
394 if (yy < ymin) ymin = yy; 396 if (yy < ymin) ymin = yy;
395 if (yy > ymax) ymax = yy; 397 if (yy > ymax) ymax = yy;
396 } 398 }
397 G4double sinTmp = sinCur; 399 G4double sinTmp = sinCur;
398 sinCur = sinCur*cosStep + cosCur*sinStep; 400 sinCur = sinCur*cosStep + cosCur*sinStep;
399 cosCur = cosCur*cosStep - sinTmp*sinStep; 401 cosCur = cosCur*cosStep - sinTmp*sinStep;
400 } 402 }
401 pMin.set(xmin,ymin,zmin); 403 pMin.set(xmin,ymin,zmin);
402 pMax.set(xmax,ymax,zmax); 404 pMax.set(xmax,ymax,zmax);
403 405
404 // Check correctness of the bounding box 406 // Check correctness of the bounding box
405 // 407 //
406 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 408 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
407 { 409 {
408 std::ostringstream message; 410 std::ostringstream message;
409 message << "Bad bounding box (min >= max) 411 message << "Bad bounding box (min >= max) for solid: "
410 << GetName() << " !" 412 << GetName() << " !"
411 << "\npMin = " << pMin 413 << "\npMin = " << pMin
412 << "\npMax = " << pMax; 414 << "\npMax = " << pMax;
413 G4Exception("G4UPolyhedra::BoundingLimits( 415 G4Exception("G4UPolyhedra::BoundingLimits()", "GeomMgt0001",
414 JustWarning, message); 416 JustWarning, message);
415 StreamInfo(G4cout); 417 StreamInfo(G4cout);
416 } 418 }
417 419
418 // Check consistency of bounding boxes 420 // Check consistency of bounding boxes
419 // 421 //
420 if (checkBBox) 422 if (checkBBox)
421 { 423 {
422 U3Vector vmin, vmax; 424 U3Vector vmin, vmax;
423 Extent(vmin,vmax); 425 Extent(vmin,vmax);
424 if (std::abs(pMin.x()-vmin.x()) > kCarTole 426 if (std::abs(pMin.x()-vmin.x()) > kCarTolerance ||
425 std::abs(pMin.y()-vmin.y()) > kCarTole 427 std::abs(pMin.y()-vmin.y()) > kCarTolerance ||
426 std::abs(pMin.z()-vmin.z()) > kCarTole 428 std::abs(pMin.z()-vmin.z()) > kCarTolerance ||
427 std::abs(pMax.x()-vmax.x()) > kCarTole 429 std::abs(pMax.x()-vmax.x()) > kCarTolerance ||
428 std::abs(pMax.y()-vmax.y()) > kCarTole 430 std::abs(pMax.y()-vmax.y()) > kCarTolerance ||
429 std::abs(pMax.z()-vmax.z()) > kCarTole 431 std::abs(pMax.z()-vmax.z()) > kCarTolerance)
430 { 432 {
431 std::ostringstream message; 433 std::ostringstream message;
432 message << "Inconsistency in bounding bo 434 message << "Inconsistency in bounding boxes for solid: "
433 << GetName() << " !" 435 << GetName() << " !"
434 << "\nBBox min: wrapper = " << p 436 << "\nBBox min: wrapper = " << pMin << " solid = " << vmin
435 << "\nBBox max: wrapper = " << p 437 << "\nBBox max: wrapper = " << pMax << " solid = " << vmax;
436 G4Exception("G4UPolyhedra::BoundingLimit 438 G4Exception("G4UPolyhedra::BoundingLimits()", "GeomMgt0001",
437 JustWarning, message); 439 JustWarning, message);
438 checkBBox = false; 440 checkBBox = false;
439 } 441 }
440 } 442 }
441 443
442 // Check consistency of angles 444 // Check consistency of angles
443 // 445 //
444 if (checkPhi) 446 if (checkPhi)
445 { 447 {
446 if (GetStartPhi() != GetPhiStart() || << 448 if (GetStartPhi() != GetPhiStart() ||
447 GetEndPhi() != GetPhiEnd() || << 449 GetEndPhi() != GetPhiEnd() ||
448 GetNumSide() != GetSideCount() || << 450 GetNumSide() != GetSideCount() ||
449 IsOpen() != (Base_t::GetPhiDelta( 451 IsOpen() != (Base_t::GetPhiDelta() < twopi))
450 { 452 {
451 std::ostringstream message; 453 std::ostringstream message;
452 message << "Inconsistency in Phi angles 454 message << "Inconsistency in Phi angles or # of sides for solid: "
453 << GetName() << " !" 455 << GetName() << " !"
454 << "\nPhi start : wrapper = " < 456 << "\nPhi start : wrapper = " << GetStartPhi()
455 << " solid = " << GetPhiStar 457 << " solid = " << GetPhiStart()
456 << "\nPhi end : wrapper = " < 458 << "\nPhi end : wrapper = " << GetEndPhi()
457 << " solid = " << GetPhiEnd( 459 << " solid = " << GetPhiEnd()
458 << "\nPhi # sides: wrapper = " < 460 << "\nPhi # sides: wrapper = " << GetNumSide()
459 << " solid = " << GetSideCou 461 << " solid = " << GetSideCount()
460 << "\nPhi is open: wrapper = " < 462 << "\nPhi is open: wrapper = " << (IsOpen() ? "true" : "false")
461 << " solid = " 463 << " solid = "
462 << ((Base_t::GetPhiDelta() < two 464 << ((Base_t::GetPhiDelta() < twopi) ? "true" : "false");
463 G4Exception("G4UPolyhedra::BoundingLimit 465 G4Exception("G4UPolyhedra::BoundingLimits()", "GeomMgt0001",
464 JustWarning, message); 466 JustWarning, message);
465 checkPhi = false; 467 checkPhi = false;
466 } 468 }
467 } 469 }
468 } 470 }
469 471
470 ////////////////////////////////////////////// 472 //////////////////////////////////////////////////////////////////////////
471 // 473 //
472 // Calculate extent under transform and specif 474 // Calculate extent under transform and specified limit
473 475
474 G4bool 476 G4bool
475 G4UPolyhedra::CalculateExtent(const EAxis pAxi 477 G4UPolyhedra::CalculateExtent(const EAxis pAxis,
476 const G4VoxelLim 478 const G4VoxelLimits& pVoxelLimit,
477 const G4AffineTr 479 const G4AffineTransform& pTransform,
478 G4double& 480 G4double& pMin, G4double& pMax) const
479 { 481 {
480 G4ThreeVector bmin, bmax; 482 G4ThreeVector bmin, bmax;
481 G4bool exist; 483 G4bool exist;
482 484
483 // Check bounding box (bbox) 485 // Check bounding box (bbox)
484 // 486 //
485 BoundingLimits(bmin,bmax); 487 BoundingLimits(bmin,bmax);
486 G4BoundingEnvelope bbox(bmin,bmax); 488 G4BoundingEnvelope bbox(bmin,bmax);
487 #ifdef G4BBOX_EXTENT 489 #ifdef G4BBOX_EXTENT
488 if (true) return bbox.CalculateExtent(pAxis, 490 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
489 #endif 491 #endif
490 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 492 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
491 { 493 {
492 return exist = pMin < pMax; << 494 return exist = (pMin < pMax) ? true : false;
493 } 495 }
494 496
495 // To find the extent, RZ contour of the pol 497 // To find the extent, RZ contour of the polycone is subdivided
496 // in triangles. The extent is calculated as 498 // in triangles. The extent is calculated as cumulative extent of
497 // all sub-polycones formed by rotation of t 499 // all sub-polycones formed by rotation of triangles around Z
498 // 500 //
499 G4TwoVectorList contourRZ; 501 G4TwoVectorList contourRZ;
500 G4TwoVectorList triangles; 502 G4TwoVectorList triangles;
501 std::vector<G4int> iout; 503 std::vector<G4int> iout;
502 G4double eminlim = pVoxelLimit.GetMinExtent( 504 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
503 G4double emaxlim = pVoxelLimit.GetMaxExtent( 505 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
504 506
505 // get RZ contour, ensure anticlockwise orde 507 // get RZ contour, ensure anticlockwise order of corners
506 for (G4int i=0; i<GetNumRZCorner(); ++i) 508 for (G4int i=0; i<GetNumRZCorner(); ++i)
507 { 509 {
508 G4PolyhedraSideRZ corner = GetCorner(i); 510 G4PolyhedraSideRZ corner = GetCorner(i);
509 contourRZ.emplace_back(corner.r,corner.z); << 511 contourRZ.push_back(G4TwoVector(corner.r,corner.z));
510 } 512 }
511 G4GeomTools::RemoveRedundantVertices(contour 513 G4GeomTools::RemoveRedundantVertices(contourRZ,iout,2*kCarTolerance);
512 G4double area = G4GeomTools::PolygonArea(con 514 G4double area = G4GeomTools::PolygonArea(contourRZ);
513 if (area < 0.) std::reverse(contourRZ.begin( 515 if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
514 516
515 // triangulate RZ countour 517 // triangulate RZ countour
516 if (!G4GeomTools::TriangulatePolygon(contour 518 if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
517 { 519 {
518 std::ostringstream message; 520 std::ostringstream message;
519 message << "Triangulation of RZ contour ha 521 message << "Triangulation of RZ contour has failed for solid: "
520 << GetName() << " !" 522 << GetName() << " !"
521 << "\nExtent has been calculated u 523 << "\nExtent has been calculated using boundary box";
522 G4Exception("G4UPolyhedra::CalculateExtent 524 G4Exception("G4UPolyhedra::CalculateExtent()",
523 "GeomMgt1002",JustWarning,mess 525 "GeomMgt1002",JustWarning,message);
524 return bbox.CalculateExtent(pAxis,pVoxelLi 526 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
525 } 527 }
526 528
527 // set trigonometric values 529 // set trigonometric values
528 G4double sphi = GetStartPhi(); 530 G4double sphi = GetStartPhi();
529 G4double ephi = GetEndPhi(); 531 G4double ephi = GetEndPhi();
530 G4double dphi = IsOpen() ? ephi-sphi : t 532 G4double dphi = IsOpen() ? ephi-sphi : twopi;
531 G4int ksteps = GetNumSide(); 533 G4int ksteps = GetNumSide();
532 G4double astep = dphi/ksteps; 534 G4double astep = dphi/ksteps;
533 G4double sinStep = std::sin(astep); 535 G4double sinStep = std::sin(astep);
534 G4double cosStep = std::cos(astep); 536 G4double cosStep = std::cos(astep);
535 G4double sinStart = GetSinStartPhi(); 537 G4double sinStart = GetSinStartPhi();
536 G4double cosStart = GetCosStartPhi(); 538 G4double cosStart = GetCosStartPhi();
537 539
538 // allocate vector lists 540 // allocate vector lists
539 std::vector<const G4ThreeVectorList *> polyg 541 std::vector<const G4ThreeVectorList *> polygons;
540 polygons.resize(ksteps+1); 542 polygons.resize(ksteps+1);
541 for (G4int k=0; k<ksteps+1; ++k) 543 for (G4int k=0; k<ksteps+1; ++k)
542 { 544 {
543 polygons[k] = new G4ThreeVectorList(3); 545 polygons[k] = new G4ThreeVectorList(3);
544 } 546 }
545 547
546 // main loop along triangles 548 // main loop along triangles
547 pMin = kInfinity; 549 pMin = kInfinity;
548 pMax = -kInfinity; 550 pMax = -kInfinity;
549 G4int ntria = triangles.size()/3; 551 G4int ntria = triangles.size()/3;
550 for (G4int i=0; i<ntria; ++i) 552 for (G4int i=0; i<ntria; ++i)
551 { 553 {
552 G4double sinCur = sinStart; 554 G4double sinCur = sinStart;
553 G4double cosCur = cosStart; 555 G4double cosCur = cosStart;
554 G4int i3 = i*3; 556 G4int i3 = i*3;
555 for (G4int k=0; k<ksteps+1; ++k) // rotate 557 for (G4int k=0; k<ksteps+1; ++k) // rotate triangle
556 { 558 {
557 auto ptr = const_cast<G4ThreeVectorList* << 559 G4ThreeVectorList* ptr = const_cast<G4ThreeVectorList*>(polygons[k]);
558 auto iter = ptr->begin(); << 560 G4ThreeVectorList::iterator iter = ptr->begin();
559 iter->set(triangles[i3+0].x()*cosCur, 561 iter->set(triangles[i3+0].x()*cosCur,
560 triangles[i3+0].x()*sinCur, 562 triangles[i3+0].x()*sinCur,
561 triangles[i3+0].y()); 563 triangles[i3+0].y());
562 iter++; 564 iter++;
563 iter->set(triangles[i3+1].x()*cosCur, 565 iter->set(triangles[i3+1].x()*cosCur,
564 triangles[i3+1].x()*sinCur, 566 triangles[i3+1].x()*sinCur,
565 triangles[i3+1].y()); 567 triangles[i3+1].y());
566 iter++; 568 iter++;
567 iter->set(triangles[i3+2].x()*cosCur, 569 iter->set(triangles[i3+2].x()*cosCur,
568 triangles[i3+2].x()*sinCur, 570 triangles[i3+2].x()*sinCur,
569 triangles[i3+2].y()); 571 triangles[i3+2].y());
570 572
571 G4double sinTmp = sinCur; 573 G4double sinTmp = sinCur;
572 sinCur = sinCur*cosStep + cosCur*sinStep 574 sinCur = sinCur*cosStep + cosCur*sinStep;
573 cosCur = cosCur*cosStep - sinTmp*sinStep 575 cosCur = cosCur*cosStep - sinTmp*sinStep;
574 } 576 }
575 577
576 // set sub-envelope and adjust extent 578 // set sub-envelope and adjust extent
577 G4double emin,emax; 579 G4double emin,emax;
578 G4BoundingEnvelope benv(polygons); 580 G4BoundingEnvelope benv(polygons);
579 if (!benv.CalculateExtent(pAxis,pVoxelLimi 581 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
580 if (emin < pMin) pMin = emin; 582 if (emin < pMin) pMin = emin;
581 if (emax > pMax) pMax = emax; 583 if (emax > pMax) pMax = emax;
582 if (eminlim > pMin && emaxlim < pMax) brea 584 if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
583 } 585 }
584 // free memory 586 // free memory
585 for (G4int k=0; k<ksteps+1; ++k) { delete po << 587 for (G4int k=0; k<ksteps+1; ++k) { delete polygons[k]; polygons[k]=0;}
586 return (pMin < pMax); 588 return (pMin < pMax);
587 } 589 }
588 590
589 591
590 ////////////////////////////////////////////// 592 ////////////////////////////////////////////////////////////////////////
591 // 593 //
592 // CreatePolyhedron 594 // CreatePolyhedron
593 // 595 //
594 G4Polyhedron* G4UPolyhedra::CreatePolyhedron() 596 G4Polyhedron* G4UPolyhedra::CreatePolyhedron() const
595 { 597 {
596 return new G4PolyhedronPgon(wrStart, wrDelta << 598 if (!IsGeneric())
>> 599 {
>> 600 return new G4PolyhedronPgon( fOriginalParameters.Start_angle,
>> 601 fOriginalParameters.Opening_angle,
>> 602 fOriginalParameters.numSide,
>> 603 fOriginalParameters.Num_z_planes,
>> 604 fOriginalParameters.Z_values,
>> 605 fOriginalParameters.Rmin,
>> 606 fOriginalParameters.Rmax);
>> 607 }
>> 608 else
>> 609 {
>> 610 // The following code prepares for:
>> 611 // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
>> 612 // const double xyz[][3],
>> 613 // const int faces_vec[][4])
>> 614 // Here is an extract from the header file HepPolyhedron.h:
>> 615 /**
>> 616 * Creates user defined polyhedron.
>> 617 * This function allows the user to define arbitrary polyhedron.
>> 618 * The faces of the polyhedron should be either triangles or planar
>> 619 * quadrilateral. Nodes of a face are defined by indexes pointing to
>> 620 * the elements in the xyz array. Numeration of the elements in the
>> 621 * array starts from 1 (like in fortran). The indexes can be positive
>> 622 * or negative. Negative sign means that the corresponding edge is
>> 623 * invisible. The normal of the face should be directed to exterior
>> 624 * of the polyhedron.
>> 625 *
>> 626 * @param Nnodes number of nodes
>> 627 * @param Nfaces number of faces
>> 628 * @param xyz nodes
>> 629 * @param faces_vec faces (quadrilaterals or triangles)
>> 630 * @return status of the operation - is non-zero in case of problem
>> 631 */
>> 632 G4int nNodes;
>> 633 G4int nFaces;
>> 634 typedef G4double double3[3];
>> 635 double3* xyz;
>> 636 typedef G4int int4[4];
>> 637 int4* faces_vec;
>> 638 if (IsOpen())
>> 639 {
>> 640 // Triangulate open ends. Simple ear-chopping algorithm...
>> 641 // I'm not sure how robust this algorithm is (J.Allison).
>> 642 //
>> 643 std::vector<G4bool> chopped(GetNumRZCorner(), false);
>> 644 std::vector<G4int*> triQuads;
>> 645 G4int remaining = GetNumRZCorner();
>> 646 G4int iStarter = 0;
>> 647 while (remaining >= 3) // Loop checking, 13.08.2015, G.Cosmo
>> 648 {
>> 649 // Find unchopped corners...
>> 650 //
>> 651 G4int A = -1, B = -1, C = -1;
>> 652 G4int iStepper = iStarter;
>> 653 do // Loop checking, 13.08.2015, G.Cosmo
>> 654 {
>> 655 if (A < 0) { A = iStepper; }
>> 656 else if (B < 0) { B = iStepper; }
>> 657 else if (C < 0) { C = iStepper; }
>> 658 do // Loop checking, 13.08.2015, G.Cosmo
>> 659 {
>> 660 if (++iStepper >= GetNumRZCorner()) iStepper = 0;
>> 661 }
>> 662 while (chopped[iStepper]);
>> 663 }
>> 664 while (C < 0 && iStepper != iStarter);
>> 665
>> 666 // Check triangle at B is pointing outward (an "ear").
>> 667 // Sign of z cross product determines...
>> 668
>> 669 G4double BAr = GetCorner(A).r - GetCorner(B).r;
>> 670 G4double BAz = GetCorner(A).z - GetCorner(B).z;
>> 671 G4double BCr = GetCorner(C).r - GetCorner(B).r;
>> 672 G4double BCz = GetCorner(C).z - GetCorner(B).z;
>> 673 if (BAr * BCz - BAz * BCr < kCarTolerance)
>> 674 {
>> 675 G4int* tq = new G4int[3];
>> 676 tq[0] = A + 1;
>> 677 tq[1] = B + 1;
>> 678 tq[2] = C + 1;
>> 679 triQuads.push_back(tq);
>> 680 chopped[B] = true;
>> 681 --remaining;
>> 682 }
>> 683 else
>> 684 {
>> 685 do // Loop checking, 13.08.2015, G.Cosmo
>> 686 {
>> 687 if (++iStarter >= GetNumRZCorner()) { iStarter = 0; }
>> 688 }
>> 689 while (chopped[iStarter]);
>> 690 }
>> 691 }
>> 692
>> 693 // Transfer to faces...
>> 694 G4int numSide=GetNumSide();
>> 695 nNodes = (numSide + 1) * GetNumRZCorner();
>> 696 nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size();
>> 697 faces_vec = new int4[nFaces];
>> 698 G4int iface = 0;
>> 699 G4int addition = GetNumRZCorner() * numSide;
>> 700 G4int d = GetNumRZCorner() - 1;
>> 701 for (G4int iEnd = 0; iEnd < 2; ++iEnd)
>> 702 {
>> 703 for (size_t i = 0; i < triQuads.size(); ++i)
>> 704 {
>> 705 // Negative for soft/auxiliary/normally invisible edges...
>> 706 //
>> 707 G4int a, b, c;
>> 708 if (iEnd == 0)
>> 709 {
>> 710 a = triQuads[i][0];
>> 711 b = triQuads[i][1];
>> 712 c = triQuads[i][2];
>> 713 }
>> 714 else
>> 715 {
>> 716 a = triQuads[i][0] + addition;
>> 717 b = triQuads[i][2] + addition;
>> 718 c = triQuads[i][1] + addition;
>> 719 }
>> 720 G4int ab = std::abs(b - a);
>> 721 G4int bc = std::abs(c - b);
>> 722 G4int ca = std::abs(a - c);
>> 723 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
>> 724 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
>> 725 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
>> 726 faces_vec[iface][3] = 0;
>> 727 ++iface;
>> 728 }
>> 729 }
>> 730
>> 731 // Continue with sides...
>> 732
>> 733 xyz = new double3[nNodes];
>> 734 const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
>> 735 G4double phi = GetStartPhi();
>> 736 G4int ixyz = 0;
>> 737 for (G4int iSide = 0; iSide < numSide; ++iSide)
>> 738 {
>> 739 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
>> 740 {
>> 741 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
>> 742 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
>> 743 xyz[ixyz][2] = GetCorner(iCorner).z;
>> 744 if (iCorner < GetNumRZCorner() - 1)
>> 745 {
>> 746 faces_vec[iface][0] = ixyz + 1;
>> 747 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
>> 748 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
>> 749 faces_vec[iface][3] = ixyz + 2;
>> 750 }
>> 751 else
>> 752 {
>> 753 faces_vec[iface][0] = ixyz + 1;
>> 754 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
>> 755 faces_vec[iface][2] = ixyz + 2;
>> 756 faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
>> 757 }
>> 758 ++iface;
>> 759 ++ixyz;
>> 760 }
>> 761 phi += dPhi;
>> 762 }
>> 763
>> 764 // Last GetCorner...
>> 765
>> 766 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
>> 767 {
>> 768 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
>> 769 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
>> 770 xyz[ixyz][2] = GetCorner(iCorner).z;
>> 771 ++ixyz;
>> 772 }
>> 773 }
>> 774 else // !phiIsOpen - i.e., a complete 360 degrees.
>> 775 {
>> 776 nNodes = GetNumSide() * GetNumRZCorner();
>> 777 nFaces = GetNumSide() * GetNumRZCorner();;
>> 778 xyz = new double3[nNodes];
>> 779 faces_vec = new int4[nFaces];
>> 780 // const G4double dPhi = (endPhi - startPhi) / numSide;
>> 781 const G4double dPhi = twopi / GetNumSide();
>> 782 // !phiIsOpen endPhi-startPhi = 360 degrees.
>> 783 G4double phi = GetStartPhi();
>> 784 G4int ixyz = 0, iface = 0;
>> 785 for (G4int iSide = 0; iSide < GetNumSide(); ++iSide)
>> 786 {
>> 787 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
>> 788 {
>> 789 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
>> 790 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
>> 791 xyz[ixyz][2] = GetCorner(iCorner).z;
>> 792 if (iSide < GetNumSide() - 1)
>> 793 {
>> 794 if (iCorner < GetNumRZCorner() - 1)
>> 795 {
>> 796 faces_vec[iface][0] = ixyz + 1;
>> 797 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
>> 798 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
>> 799 faces_vec[iface][3] = ixyz + 2;
>> 800 }
>> 801 else
>> 802 {
>> 803 faces_vec[iface][0] = ixyz + 1;
>> 804 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
>> 805 faces_vec[iface][2] = ixyz + 2;
>> 806 faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
>> 807 }
>> 808 }
>> 809 else // Last side joins ends...
>> 810 {
>> 811 if (iCorner < GetNumRZCorner() - 1)
>> 812 {
>> 813 faces_vec[iface][0] = ixyz + 1;
>> 814 faces_vec[iface][1] = ixyz + GetNumRZCorner() - nFaces + 1;
>> 815 faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2;
>> 816 faces_vec[iface][3] = ixyz + 2;
>> 817 }
>> 818 else
>> 819 {
>> 820 faces_vec[iface][0] = ixyz + 1;
>> 821 faces_vec[iface][1] = ixyz - nFaces + GetNumRZCorner() + 1;
>> 822 faces_vec[iface][2] = ixyz - nFaces + 2;
>> 823 faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
>> 824 }
>> 825 }
>> 826 ++ixyz;
>> 827 ++iface;
>> 828 }
>> 829 phi += dPhi;
>> 830 }
>> 831 }
>> 832 G4Polyhedron* polyhedron = new G4Polyhedron;
>> 833 G4int prob = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
>> 834 delete [] faces_vec;
>> 835 delete [] xyz;
>> 836 if (prob)
>> 837 {
>> 838 std::ostringstream message;
>> 839 message << "Problem creating G4Polyhedron for: " << GetName();
>> 840 G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
>> 841 JustWarning, message);
>> 842 delete polyhedron;
>> 843 return nullptr;
>> 844 }
>> 845 else
>> 846 {
>> 847 return polyhedron;
>> 848 }
>> 849 }
597 } 850 }
598 851
599 #endif // G4GEOM_USE_USOLIDS 852 #endif // G4GEOM_USE_USOLIDS
600 853