| 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 // G4GenericPolycone implementation 26 // G4GenericPolycone implementation
27 // 27 //
28 // Authors: T.Nikitina, G.Cosmo - CERN 28 // Authors: T.Nikitina, G.Cosmo - CERN
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4GenericPolycone.hh" 31 #include "G4GenericPolycone.hh"
32 32
33 #if !defined(G4GEOM_USE_UGENERICPOLYCONE) 33 #if !defined(G4GEOM_USE_UGENERICPOLYCONE)
34 34
35 #include "G4PolyconeSide.hh" 35 #include "G4PolyconeSide.hh"
36 #include "G4PolyPhiFace.hh" 36 #include "G4PolyPhiFace.hh"
37 37
38 #include "G4GeomTools.hh" 38 #include "G4GeomTools.hh"
39 #include "G4VoxelLimits.hh" 39 #include "G4VoxelLimits.hh"
40 #include "G4AffineTransform.hh" 40 #include "G4AffineTransform.hh"
41 #include "G4BoundingEnvelope.hh" 41 #include "G4BoundingEnvelope.hh"
42 42
43 #include "G4QuickRand.hh" << 43 #include "Randomize.hh"
44 44
45 #include "G4Polyhedron.hh" 45 #include "G4Polyhedron.hh"
46 #include "G4EnclosingCylinder.hh" 46 #include "G4EnclosingCylinder.hh"
47 #include "G4ReduciblePolygon.hh" 47 #include "G4ReduciblePolygon.hh"
48 #include "G4VPVParameterisation.hh" 48 #include "G4VPVParameterisation.hh"
49 49
50 namespace <<
51 { <<
52 G4Mutex surface_elementsMutex = G4MUTEX_INIT <<
53 } <<
54 <<
55 using namespace CLHEP; 50 using namespace CLHEP;
56 51
57 // Constructor (generic parameters) 52 // Constructor (generic parameters)
58 // 53 //
59 G4GenericPolycone::G4GenericPolycone( const G4 << 54 G4GenericPolycone::G4GenericPolycone( const G4String& name,
60 G4double phiStar 55 G4double phiStart,
61 G4double phiTota 56 G4double phiTotal,
62 G4int numRZ, 57 G4int numRZ,
63 const G4double r[], 58 const G4double r[],
64 const G4double z[] ) 59 const G4double z[] )
65 : G4VCSGfaceted( name ) 60 : G4VCSGfaceted( name )
66 { << 61 {
67 << 62
68 auto rz = new G4ReduciblePolygon( r, z, numR << 63 G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
69 << 64
70 Create( phiStart, phiTotal, rz ); 65 Create( phiStart, phiTotal, rz );
71 << 66
72 // Set original_parameters struct for consis 67 // Set original_parameters struct for consistency
73 // 68 //
74 //SetOriginalParameters(rz); 69 //SetOriginalParameters(rz);
75 70
76 delete rz; 71 delete rz;
77 } 72 }
78 73
79 // Create 74 // Create
80 // 75 //
81 // Generic create routine, called by each cons 76 // Generic create routine, called by each constructor after
82 // conversion of arguments 77 // conversion of arguments
83 // 78 //
84 void G4GenericPolycone::Create( G4double phiSt 79 void G4GenericPolycone::Create( G4double phiStart,
85 G4double phiTo << 80 G4double phiTotal,
86 G4ReduciblePol << 81 G4ReduciblePolygon *rz )
87 { 82 {
88 // 83 //
89 // Perform checks of rz values 84 // Perform checks of rz values
90 // 85 //
91 if (rz->Amin() < 0.0) 86 if (rz->Amin() < 0.0)
92 { 87 {
93 std::ostringstream message; 88 std::ostringstream message;
94 message << "Illegal input parameters - " < 89 message << "Illegal input parameters - " << GetName() << G4endl
95 << " All R values must be > 90 << " All R values must be >= 0 !";
96 G4Exception("G4GenericPolycone::Create()", 91 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
97 FatalErrorInArgument, message) 92 FatalErrorInArgument, message);
98 } 93 }
99 << 94
100 G4double rzArea = rz->Area(); 95 G4double rzArea = rz->Area();
101 if (rzArea < -kCarTolerance) 96 if (rzArea < -kCarTolerance)
102 { 97 {
103 rz->ReverseOrder(); 98 rz->ReverseOrder();
104 } 99 }
105 else if (rzArea < kCarTolerance) 100 else if (rzArea < kCarTolerance)
106 { 101 {
107 std::ostringstream message; 102 std::ostringstream message;
108 message << "Illegal input parameters - " < 103 message << "Illegal input parameters - " << GetName() << G4endl
109 << " R/Z cross section is z 104 << " R/Z cross section is zero or near zero: " << rzArea;
110 G4Exception("G4GenericPolycone::Create()", 105 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
111 FatalErrorInArgument, message) 106 FatalErrorInArgument, message);
112 } 107 }
113 << 108
114 if ( (!rz->RemoveDuplicateVertices( kCarTole 109 if ( (!rz->RemoveDuplicateVertices( kCarTolerance ))
115 || (!rz->RemoveRedundantVertices( kCarTole << 110 || (!rz->RemoveRedundantVertices( kCarTolerance )) )
116 { 111 {
117 std::ostringstream message; 112 std::ostringstream message;
118 message << "Illegal input parameters - " < 113 message << "Illegal input parameters - " << GetName() << G4endl
119 << " Too few unique R/Z val 114 << " Too few unique R/Z values !";
120 G4Exception("G4GenericPolycone::Create()", 115 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
121 FatalErrorInArgument, message) 116 FatalErrorInArgument, message);
122 } 117 }
123 118
124 if (rz->CrossesItself(1/kInfinity)) << 119 if (rz->CrossesItself(1/kInfinity))
125 { 120 {
126 std::ostringstream message; 121 std::ostringstream message;
127 message << "Illegal input parameters - " < 122 message << "Illegal input parameters - " << GetName() << G4endl
128 << " R/Z segments cross !"; 123 << " R/Z segments cross !";
129 G4Exception("G4GenericPolycone::Create()", 124 G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
130 FatalErrorInArgument, message) 125 FatalErrorInArgument, message);
131 } 126 }
132 127
133 numCorner = rz->NumVertices(); 128 numCorner = rz->NumVertices();
134 129
135 // 130 //
136 // Phi opening? Account for some possible ro 131 // Phi opening? Account for some possible roundoff, and interpret
137 // nonsense value as representing no phi ope 132 // nonsense value as representing no phi opening
138 // 133 //
139 if (phiTotal <= 0 || phiTotal > twopi-1E-10) 134 if (phiTotal <= 0 || phiTotal > twopi-1E-10)
140 { 135 {
141 phiIsOpen = false; 136 phiIsOpen = false;
142 startPhi = 0; 137 startPhi = 0;
143 endPhi = twopi; 138 endPhi = twopi;
144 } 139 }
145 else 140 else
146 { 141 {
147 phiIsOpen = true; 142 phiIsOpen = true;
148 << 143
149 // 144 //
150 // Convert phi into our convention 145 // Convert phi into our convention
151 // 146 //
152 startPhi = phiStart; 147 startPhi = phiStart;
153 while( startPhi < 0 ) // Loop checking, 148 while( startPhi < 0 ) // Loop checking, 13.08.2015, G.Cosmo
154 startPhi += twopi; 149 startPhi += twopi;
155 << 150
156 endPhi = phiStart+phiTotal; 151 endPhi = phiStart+phiTotal;
157 while( endPhi < startPhi ) // Loop chec 152 while( endPhi < startPhi ) // Loop checking, 13.08.2015, G.Cosmo
158 endPhi += twopi; 153 endPhi += twopi;
159 } 154 }
160 << 155
161 // 156 //
162 // Allocate corner array. << 157 // Allocate corner array.
163 // 158 //
164 corners = new G4PolyconeSideRZ[numCorner]; 159 corners = new G4PolyconeSideRZ[numCorner];
165 160
166 // 161 //
167 // Copy corners 162 // Copy corners
168 // 163 //
169 G4ReduciblePolygonIterator iterRZ(rz); 164 G4ReduciblePolygonIterator iterRZ(rz);
170 << 165
171 G4PolyconeSideRZ* next = corners; 166 G4PolyconeSideRZ* next = corners;
172 iterRZ.Begin(); 167 iterRZ.Begin();
173 do // Loop checking, 13.08.2015, G.Cosmo 168 do // Loop checking, 13.08.2015, G.Cosmo
174 { 169 {
175 next->r = iterRZ.GetA(); 170 next->r = iterRZ.GetA();
176 next->z = iterRZ.GetB(); 171 next->z = iterRZ.GetB();
177 } while( ++next, iterRZ.Next() ); 172 } while( ++next, iterRZ.Next() );
178 << 173
179 // 174 //
180 // Allocate face pointer array 175 // Allocate face pointer array
181 // 176 //
182 numFace = phiIsOpen ? numCorner+2 : numCorne 177 numFace = phiIsOpen ? numCorner+2 : numCorner;
183 faces = new G4VCSGface*[numFace]; 178 faces = new G4VCSGface*[numFace];
184 << 179
185 // 180 //
186 // Construct conical faces 181 // Construct conical faces
187 // 182 //
188 // But! Don't construct a face if both point 183 // But! Don't construct a face if both points are at zero radius!
189 // 184 //
190 G4PolyconeSideRZ *corner = corners, 185 G4PolyconeSideRZ *corner = corners,
191 *prev = corners + numCorner 186 *prev = corners + numCorner-1,
192 *nextNext; 187 *nextNext;
193 G4VCSGface** face = faces; 188 G4VCSGface** face = faces;
194 do // Loop checking, 13.08.2015, G.Cosmo 189 do // Loop checking, 13.08.2015, G.Cosmo
195 { 190 {
196 next = corner+1; 191 next = corner+1;
197 if (next >= corners+numCorner) next = corn 192 if (next >= corners+numCorner) next = corners;
198 nextNext = next+1; 193 nextNext = next+1;
199 if (nextNext >= corners+numCorner) nextNex 194 if (nextNext >= corners+numCorner) nextNext = corners;
200 << 195
201 if (corner->r < 1/kInfinity && next->r < 1 196 if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
202 << 197
203 // 198 //
204 // We must decide here if we can dare decl 199 // We must decide here if we can dare declare one of our faces
205 // as having a "valid" normal (i.e. allBeh 200 // as having a "valid" normal (i.e. allBehind = true). This
206 // is never possible if the face faces "in 201 // is never possible if the face faces "inward" in r.
207 // 202 //
208 G4bool allBehind; 203 G4bool allBehind;
209 if (corner->z > next->z) 204 if (corner->z > next->z)
210 { 205 {
211 allBehind = false; 206 allBehind = false;
212 } 207 }
213 else 208 else
214 { 209 {
215 // 210 //
216 // Otherwise, it is only true if the lin 211 // Otherwise, it is only true if the line passing
217 // through the two points of the segment 212 // through the two points of the segment do not
218 // split the r/z cross section 213 // split the r/z cross section
219 // 214 //
220 allBehind = !rz->BisectedBy( corner->r, 215 allBehind = !rz->BisectedBy( corner->r, corner->z,
221 next->r, next->z, kCarToleran 216 next->r, next->z, kCarTolerance );
222 } 217 }
223 << 218
224 *face++ = new G4PolyconeSide( prev, corner 219 *face++ = new G4PolyconeSide( prev, corner, next, nextNext,
225 startPhi, endPhi-startPhi, phi 220 startPhi, endPhi-startPhi, phiIsOpen, allBehind );
226 } while( prev=corner, corner=next, corner > 221 } while( prev=corner, corner=next, corner > corners );
227 << 222
228 if (phiIsOpen) 223 if (phiIsOpen)
229 { 224 {
230 // 225 //
231 // Construct phi open edges 226 // Construct phi open edges
232 // 227 //
233 *face++ = new G4PolyPhiFace( rz, startPhi, 228 *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi );
234 *face++ = new G4PolyPhiFace( rz, endPhi, 229 *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi );
235 } 230 }
236 << 231
237 // 232 //
238 // We might have dropped a face or two: reca 233 // We might have dropped a face or two: recalculate numFace
239 // 234 //
240 numFace = (G4int)(face-faces); << 235 numFace = face-faces;
241 << 236
242 // 237 //
243 // Make enclosingCylinder 238 // Make enclosingCylinder
244 // 239 //
245 enclosingCylinder = 240 enclosingCylinder =
246 new G4EnclosingCylinder( rz, phiIsOpen, ph 241 new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
247 } 242 }
248 243
249 // Fake default constructor - sets only member 244 // Fake default constructor - sets only member data and allocates memory
250 // for usage restri 245 // for usage restricted to object persistency.
251 // 246 //
252 G4GenericPolycone::G4GenericPolycone( __void__ 247 G4GenericPolycone::G4GenericPolycone( __void__& a )
253 : G4VCSGfaceted(a), startPhi(0.), endPhi(0.) 248 : G4VCSGfaceted(a), startPhi(0.), endPhi(0.), numCorner(0)
254 { 249 {
255 } 250 }
256 251
257 // Destructor 252 // Destructor
258 // 253 //
259 G4GenericPolycone::~G4GenericPolycone() 254 G4GenericPolycone::~G4GenericPolycone()
260 { 255 {
261 delete [] corners; 256 delete [] corners;
262 delete enclosingCylinder; 257 delete enclosingCylinder;
263 delete fElements; <<
264 delete fpPolyhedron; <<
265 corners = nullptr; <<
266 enclosingCylinder = nullptr; <<
267 fElements = nullptr; <<
268 fpPolyhedron = nullptr; <<
269 } 258 }
270 259
271 // Copy constructor 260 // Copy constructor
272 // 261 //
273 G4GenericPolycone::G4GenericPolycone( const G4 262 G4GenericPolycone::G4GenericPolycone( const G4GenericPolycone& source )
274 : G4VCSGfaceted( source ) 263 : G4VCSGfaceted( source )
275 { 264 {
276 CopyStuff( source ); 265 CopyStuff( source );
277 } 266 }
278 267
279 // Assignment operator 268 // Assignment operator
280 // 269 //
281 G4GenericPolycone& 270 G4GenericPolycone&
282 G4GenericPolycone::operator=( const G4GenericP 271 G4GenericPolycone::operator=( const G4GenericPolycone& source )
283 { 272 {
284 if (this == &source) return *this; 273 if (this == &source) return *this;
285 << 274
286 G4VCSGfaceted::operator=( source ); 275 G4VCSGfaceted::operator=( source );
287 << 276
288 delete [] corners; 277 delete [] corners;
289 // if (original_parameters) delete original_ 278 // if (original_parameters) delete original_parameters;
290 << 279
291 delete enclosingCylinder; 280 delete enclosingCylinder;
292 << 281
293 CopyStuff( source ); 282 CopyStuff( source );
294 << 283
295 return *this; 284 return *this;
296 } 285 }
297 286
298 // CopyStuff 287 // CopyStuff
299 // 288 //
300 void G4GenericPolycone::CopyStuff( const G4Gen 289 void G4GenericPolycone::CopyStuff( const G4GenericPolycone& source )
301 { 290 {
302 // 291 //
303 // Simple stuff 292 // Simple stuff
304 // 293 //
305 startPhi = source.startPhi; 294 startPhi = source.startPhi;
306 endPhi = source.endPhi; 295 endPhi = source.endPhi;
307 phiIsOpen = source.phiIsOpen; 296 phiIsOpen = source.phiIsOpen;
308 numCorner = source.numCorner; 297 numCorner = source.numCorner;
309 298
310 // 299 //
311 // The corner array 300 // The corner array
312 // 301 //
313 corners = new G4PolyconeSideRZ[numCorner]; 302 corners = new G4PolyconeSideRZ[numCorner];
314 << 303
315 G4PolyconeSideRZ *corn = corners, 304 G4PolyconeSideRZ *corn = corners,
316 *sourceCorn = source.corners; 305 *sourceCorn = source.corners;
317 do // Loop checking, 13.08.2015, G.Cosmo 306 do // Loop checking, 13.08.2015, G.Cosmo
318 { 307 {
319 *corn = *sourceCorn; 308 *corn = *sourceCorn;
320 } while( ++sourceCorn, ++corn < corners+numC 309 } while( ++sourceCorn, ++corn < corners+numCorner );
321 << 310
322 // 311 //
323 // Enclosing cylinder 312 // Enclosing cylinder
324 // 313 //
325 enclosingCylinder = new G4EnclosingCylinder( 314 enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder );
326 315
327 // <<
328 // Surface elements <<
329 // <<
330 delete fElements; <<
331 fElements = nullptr; <<
332 <<
333 // Polyhedron <<
334 // <<
335 fRebuildPolyhedron = false; 316 fRebuildPolyhedron = false;
336 delete fpPolyhedron; <<
337 fpPolyhedron = nullptr; 317 fpPolyhedron = nullptr;
338 } 318 }
339 319
340 // Reset 320 // Reset
341 // 321 //
342 G4bool G4GenericPolycone::Reset() 322 G4bool G4GenericPolycone::Reset()
343 { 323 {
344 std::ostringstream message; 324 std::ostringstream message;
345 message << "Solid " << GetName() << " built 325 message << "Solid " << GetName() << " built using generic construct."
346 << G4endl << "Not applicable to the 326 << G4endl << "Not applicable to the generic construct !";
347 G4Exception("G4GenericPolycone::Reset()", "G 327 G4Exception("G4GenericPolycone::Reset()", "GeomSolids1001",
348 JustWarning, message, "Parameter 328 JustWarning, message, "Parameters NOT resetted.");
349 return true; 329 return true;
350 } 330 }
351 331
352 // Inside 332 // Inside
353 // 333 //
354 // This is an override of G4VCSGfaceted::Insid 334 // This is an override of G4VCSGfaceted::Inside, created in order
355 // to speed things up by first checking with G 335 // to speed things up by first checking with G4EnclosingCylinder.
356 // 336 //
357 EInside G4GenericPolycone::Inside( const G4Thr 337 EInside G4GenericPolycone::Inside( const G4ThreeVector& p ) const
358 { 338 {
359 // 339 //
360 // Quick test 340 // Quick test
361 // 341 //
362 if (enclosingCylinder->MustBeOutside(p)) ret 342 if (enclosingCylinder->MustBeOutside(p)) return kOutside;
363 343
364 // 344 //
365 // Long answer 345 // Long answer
366 // 346 //
367 return G4VCSGfaceted::Inside(p); 347 return G4VCSGfaceted::Inside(p);
368 } 348 }
369 349
370 // DistanceToIn 350 // DistanceToIn
371 // 351 //
372 // This is an override of G4VCSGfaceted::Insid 352 // This is an override of G4VCSGfaceted::Inside, created in order
373 // to speed things up by first checking with G 353 // to speed things up by first checking with G4EnclosingCylinder.
374 // 354 //
375 G4double G4GenericPolycone::DistanceToIn( cons 355 G4double G4GenericPolycone::DistanceToIn( const G4ThreeVector& p,
376 cons 356 const G4ThreeVector& v ) const
377 { 357 {
378 // 358 //
379 // Quick test 359 // Quick test
380 // 360 //
381 if (enclosingCylinder->ShouldMiss(p,v)) 361 if (enclosingCylinder->ShouldMiss(p,v))
382 return kInfinity; 362 return kInfinity;
383 << 363
384 // 364 //
385 // Long answer 365 // Long answer
386 // 366 //
387 return G4VCSGfaceted::DistanceToIn( p, v ); 367 return G4VCSGfaceted::DistanceToIn( p, v );
388 } 368 }
389 369
390 // DistanceToIn 370 // DistanceToIn
391 // 371 //
392 G4double G4GenericPolycone::DistanceToIn( cons 372 G4double G4GenericPolycone::DistanceToIn( const G4ThreeVector& p ) const
393 { 373 {
394 return G4VCSGfaceted::DistanceToIn(p); 374 return G4VCSGfaceted::DistanceToIn(p);
395 } 375 }
396 376
397 // BoundingLimits 377 // BoundingLimits
398 // 378 //
399 // Get bounding box 379 // Get bounding box
400 // 380 //
401 void 381 void
402 G4GenericPolycone::BoundingLimits(G4ThreeVecto 382 G4GenericPolycone::BoundingLimits(G4ThreeVector& pMin,
403 G4ThreeVecto 383 G4ThreeVector& pMax) const
404 { 384 {
405 G4double rmin = kInfinity, rmax = -kInfinity 385 G4double rmin = kInfinity, rmax = -kInfinity;
406 G4double zmin = kInfinity, zmax = -kInfinity 386 G4double zmin = kInfinity, zmax = -kInfinity;
407 387
408 for (G4int i=0; i<GetNumRZCorner(); ++i) 388 for (G4int i=0; i<GetNumRZCorner(); ++i)
409 { 389 {
410 G4PolyconeSideRZ corner = GetCorner(i); 390 G4PolyconeSideRZ corner = GetCorner(i);
411 if (corner.r < rmin) rmin = corner.r; 391 if (corner.r < rmin) rmin = corner.r;
412 if (corner.r > rmax) rmax = corner.r; 392 if (corner.r > rmax) rmax = corner.r;
413 if (corner.z < zmin) zmin = corner.z; 393 if (corner.z < zmin) zmin = corner.z;
414 if (corner.z > zmax) zmax = corner.z; 394 if (corner.z > zmax) zmax = corner.z;
415 } 395 }
416 396
417 if (IsOpen()) 397 if (IsOpen())
418 { 398 {
419 G4TwoVector vmin,vmax; 399 G4TwoVector vmin,vmax;
420 G4GeomTools::DiskExtent(rmin,rmax, 400 G4GeomTools::DiskExtent(rmin,rmax,
421 GetSinStartPhi(),G 401 GetSinStartPhi(),GetCosStartPhi(),
422 GetSinEndPhi(),Get 402 GetSinEndPhi(),GetCosEndPhi(),
423 vmin,vmax); 403 vmin,vmax);
424 pMin.set(vmin.x(),vmin.y(),zmin); 404 pMin.set(vmin.x(),vmin.y(),zmin);
425 pMax.set(vmax.x(),vmax.y(),zmax); 405 pMax.set(vmax.x(),vmax.y(),zmax);
426 } 406 }
427 else 407 else
428 { 408 {
429 pMin.set(-rmax,-rmax, zmin); 409 pMin.set(-rmax,-rmax, zmin);
430 pMax.set( rmax, rmax, zmax); 410 pMax.set( rmax, rmax, zmax);
431 } 411 }
432 412
433 // Check correctness of the bounding box 413 // Check correctness of the bounding box
434 // 414 //
435 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 415 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
436 { 416 {
437 std::ostringstream message; 417 std::ostringstream message;
438 message << "Bad bounding box (min >= max) 418 message << "Bad bounding box (min >= max) for solid: "
439 << GetName() << " !" 419 << GetName() << " !"
440 << "\npMin = " << pMin 420 << "\npMin = " << pMin
441 << "\npMax = " << pMax; 421 << "\npMax = " << pMax;
442 G4Exception("GenericG4Polycone::BoundingLi 422 G4Exception("GenericG4Polycone::BoundingLimits()", "GeomMgt0001",
443 JustWarning, message); 423 JustWarning, message);
444 DumpInfo(); 424 DumpInfo();
445 } 425 }
446 } 426 }
447 427
448 // CalculateExtent 428 // CalculateExtent
449 // 429 //
450 // Calculate extent under transform and specif 430 // Calculate extent under transform and specified limit
451 // 431 //
452 G4bool 432 G4bool
453 G4GenericPolycone::CalculateExtent(const EAxis 433 G4GenericPolycone::CalculateExtent(const EAxis pAxis,
454 const G4Vox 434 const G4VoxelLimits& pVoxelLimit,
455 const G4Aff 435 const G4AffineTransform& pTransform,
456 G4dou 436 G4double& pMin, G4double& pMax) const
457 { 437 {
458 G4ThreeVector bmin, bmax; 438 G4ThreeVector bmin, bmax;
459 G4bool exist; 439 G4bool exist;
460 440
461 // Check bounding box (bbox) 441 // Check bounding box (bbox)
462 // 442 //
463 BoundingLimits(bmin,bmax); 443 BoundingLimits(bmin,bmax);
464 G4BoundingEnvelope bbox(bmin,bmax); 444 G4BoundingEnvelope bbox(bmin,bmax);
465 #ifdef G4BBOX_EXTENT 445 #ifdef G4BBOX_EXTENT
466 return bbox.CalculateExtent(pAxis,pVoxelLimi 446 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
467 #endif 447 #endif
468 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 448 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
469 { 449 {
470 return exist = pMin < pMax; << 450 return exist = (pMin < pMax) ? true : false;
471 } 451 }
472 452
473 // To find the extent, RZ contour of the pol 453 // To find the extent, RZ contour of the polycone is subdivided
474 // in triangles. The extent is calculated as 454 // in triangles. The extent is calculated as cumulative extent of
475 // all sub-polycones formed by rotation of t 455 // all sub-polycones formed by rotation of triangles around Z
476 // 456 //
477 G4TwoVectorList contourRZ; 457 G4TwoVectorList contourRZ;
478 G4TwoVectorList triangles; 458 G4TwoVectorList triangles;
479 G4double eminlim = pVoxelLimit.GetMinExtent( 459 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
480 G4double emaxlim = pVoxelLimit.GetMaxExtent( 460 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
481 461
482 // get RZ contour, ensure anticlockwise orde 462 // get RZ contour, ensure anticlockwise order of corners
483 for (G4int i=0; i<GetNumRZCorner(); ++i) 463 for (G4int i=0; i<GetNumRZCorner(); ++i)
484 { 464 {
485 G4PolyconeSideRZ corner = GetCorner(i); 465 G4PolyconeSideRZ corner = GetCorner(i);
486 contourRZ.emplace_back(corner.r,corner.z); << 466 contourRZ.push_back(G4TwoVector(corner.r,corner.z));
487 } 467 }
488 G4double area = G4GeomTools::PolygonArea(con 468 G4double area = G4GeomTools::PolygonArea(contourRZ);
489 if (area < 0.) std::reverse(contourRZ.begin( 469 if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
490 470
491 // triangulate RZ countour 471 // triangulate RZ countour
492 if (!G4GeomTools::TriangulatePolygon(contour 472 if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
493 { 473 {
494 std::ostringstream message; 474 std::ostringstream message;
495 message << "Triangulation of RZ contour ha 475 message << "Triangulation of RZ contour has failed for solid: "
496 << GetName() << " !" 476 << GetName() << " !"
497 << "\nExtent has been calculated u 477 << "\nExtent has been calculated using boundary box";
498 G4Exception("G4GenericPolycone::CalculateE 478 G4Exception("G4GenericPolycone::CalculateExtent()",
499 "GeomMgt1002", JustWarning, me 479 "GeomMgt1002", JustWarning, message);
500 return bbox.CalculateExtent(pAxis,pVoxelLi 480 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
501 } 481 }
502 482
503 // set trigonometric values 483 // set trigonometric values
504 const G4int NSTEPS = 24; // numbe 484 const G4int NSTEPS = 24; // number of steps for whole circle
505 G4double astep = twopi/NSTEPS; // max a 485 G4double astep = twopi/NSTEPS; // max angle for one step
506 486
507 G4double sphi = GetStartPhi(); 487 G4double sphi = GetStartPhi();
508 G4double ephi = GetEndPhi(); 488 G4double ephi = GetEndPhi();
509 G4double dphi = IsOpen() ? ephi-sphi : two 489 G4double dphi = IsOpen() ? ephi-sphi : twopi;
510 G4int ksteps = (dphi <= astep) ? 1 : (G4i 490 G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
511 G4double ang = dphi/ksteps; 491 G4double ang = dphi/ksteps;
512 492
513 G4double sinHalf = std::sin(0.5*ang); 493 G4double sinHalf = std::sin(0.5*ang);
514 G4double cosHalf = std::cos(0.5*ang); 494 G4double cosHalf = std::cos(0.5*ang);
515 G4double sinStep = 2.*sinHalf*cosHalf; 495 G4double sinStep = 2.*sinHalf*cosHalf;
516 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 496 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
517 497
518 G4double sinStart = GetSinStartPhi(); 498 G4double sinStart = GetSinStartPhi();
519 G4double cosStart = GetCosStartPhi(); 499 G4double cosStart = GetCosStartPhi();
520 G4double sinEnd = GetSinEndPhi(); 500 G4double sinEnd = GetSinEndPhi();
521 G4double cosEnd = GetCosEndPhi(); 501 G4double cosEnd = GetCosEndPhi();
522 502
523 // define vectors and arrays 503 // define vectors and arrays
524 std::vector<const G4ThreeVectorList *> polyg 504 std::vector<const G4ThreeVectorList *> polygons;
525 polygons.resize(ksteps+2); 505 polygons.resize(ksteps+2);
526 G4ThreeVectorList pols[NSTEPS+2]; 506 G4ThreeVectorList pols[NSTEPS+2];
527 for (G4int k=0; k<ksteps+2; ++k) pols[k].res 507 for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(6);
528 for (G4int k=0; k<ksteps+2; ++k) polygons[k] 508 for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
529 G4double r0[6],z0[6]; // contour with origin 509 G4double r0[6],z0[6]; // contour with original edges of triangle
530 G4double r1[6]; // shifted radii of ex 510 G4double r1[6]; // shifted radii of external edges of triangle
531 511
532 // main loop along triangles 512 // main loop along triangles
533 pMin = kInfinity; 513 pMin = kInfinity;
534 pMax =-kInfinity; 514 pMax =-kInfinity;
535 G4int ntria = (G4int)triangles.size()/3; << 515 G4int ntria = triangles.size()/3;
536 for (G4int i=0; i<ntria; ++i) 516 for (G4int i=0; i<ntria; ++i)
537 { 517 {
538 G4int i3 = i*3; 518 G4int i3 = i*3;
539 for (G4int k=0; k<3; ++k) 519 for (G4int k=0; k<3; ++k)
540 { 520 {
541 G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3; 521 G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3;
542 G4int k2 = k*2; 522 G4int k2 = k*2;
543 // set contour with original edges of tr 523 // set contour with original edges of triangle
544 r0[k2+0] = triangles[e0].x(); z0[k2+0] = 524 r0[k2+0] = triangles[e0].x(); z0[k2+0] = triangles[e0].y();
545 r0[k2+1] = triangles[e1].x(); z0[k2+1] = 525 r0[k2+1] = triangles[e1].x(); z0[k2+1] = triangles[e1].y();
546 // set shifted radii 526 // set shifted radii
547 r1[k2+0] = r0[k2+0]; 527 r1[k2+0] = r0[k2+0];
548 r1[k2+1] = r0[k2+1]; 528 r1[k2+1] = r0[k2+1];
549 if (z0[k2+1] - z0[k2+0] <= 0) continue; 529 if (z0[k2+1] - z0[k2+0] <= 0) continue;
550 r1[k2+0] /= cosHalf; 530 r1[k2+0] /= cosHalf;
551 r1[k2+1] /= cosHalf; 531 r1[k2+1] /= cosHalf;
552 } 532 }
553 533
554 // rotate countour, set sequence of 6-side 534 // rotate countour, set sequence of 6-sided polygons
555 G4double sinCur = sinStart*cosHalf + cosSt 535 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
556 G4double cosCur = cosStart*cosHalf - sinSt 536 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
557 for (G4int j=0; j<6; ++j) 537 for (G4int j=0; j<6; ++j)
558 { 538 {
559 pols[0][j].set(r0[j]*cosStart,r0[j]*sinS 539 pols[0][j].set(r0[j]*cosStart,r0[j]*sinStart,z0[j]);
560 } 540 }
561 for (G4int k=1; k<ksteps+1; ++k) 541 for (G4int k=1; k<ksteps+1; ++k)
562 { 542 {
563 for (G4int j=0; j<6; ++j) 543 for (G4int j=0; j<6; ++j)
564 { 544 {
565 pols[k][j].set(r1[j]*cosCur,r1[j]*sinC 545 pols[k][j].set(r1[j]*cosCur,r1[j]*sinCur,z0[j]);
566 } 546 }
567 G4double sinTmp = sinCur; 547 G4double sinTmp = sinCur;
568 sinCur = sinCur*cosStep + cosCur*sinStep 548 sinCur = sinCur*cosStep + cosCur*sinStep;
569 cosCur = cosCur*cosStep - sinTmp*sinStep 549 cosCur = cosCur*cosStep - sinTmp*sinStep;
570 } 550 }
571 for (G4int j=0; j<6; ++j) 551 for (G4int j=0; j<6; ++j)
572 { 552 {
573 pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j] 553 pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j]*sinEnd,z0[j]);
574 } 554 }
575 555
576 // set sub-envelope and adjust extent 556 // set sub-envelope and adjust extent
577 G4double emin,emax; 557 G4double emin,emax;
578 G4BoundingEnvelope benv(polygons); 558 G4BoundingEnvelope benv(polygons);
579 if (!benv.CalculateExtent(pAxis,pVoxelLimi 559 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
580 if (emin < pMin) pMin = emin; 560 if (emin < pMin) pMin = emin;
581 if (emax > pMax) pMax = emax; 561 if (emax > pMax) pMax = emax;
582 if (eminlim > pMin && emaxlim < pMax) retu 562 if (eminlim > pMin && emaxlim < pMax) return true; // max possible extent
583 } 563 }
584 return (pMin < pMax); 564 return (pMin < pMax);
585 } 565 }
586 566
587 // GetEntityType 567 // GetEntityType
588 // 568 //
589 G4GeometryType G4GenericPolycone::GetEntityTy 569 G4GeometryType G4GenericPolycone::GetEntityType() const
590 { 570 {
591 return {"G4GenericPolycone"}; << 571 return G4String("G4GenericPolycone");
592 } 572 }
593 573
594 // Clone 574 // Clone
595 // 575 //
596 // Make a clone of the object 576 // Make a clone of the object
597 // 577 //
598 G4VSolid* G4GenericPolycone::Clone() const 578 G4VSolid* G4GenericPolycone::Clone() const
599 { 579 {
600 return new G4GenericPolycone(*this); 580 return new G4GenericPolycone(*this);
601 } 581 }
602 582
603 // StreamInfo 583 // StreamInfo
604 // 584 //
605 // Stream object contents to an output stream 585 // Stream object contents to an output stream
606 // 586 //
607 std::ostream& G4GenericPolycone::StreamInfo( s 587 std::ostream& G4GenericPolycone::StreamInfo( std::ostream& os ) const
608 { 588 {
609 G4long oldprc = os.precision(16); << 589 G4int oldprc = os.precision(16);
610 os << "------------------------------------- 590 os << "-----------------------------------------------------------\n"
611 << " *** Dump for solid - " << GetName 591 << " *** Dump for solid - " << GetName() << " ***\n"
612 << " ================================= 592 << " ===================================================\n"
613 << " Solid type: G4GenericPolycone\n" 593 << " Solid type: G4GenericPolycone\n"
614 << " Parameters: \n" 594 << " Parameters: \n"
615 << " starting phi angle : " << startPh 595 << " starting phi angle : " << startPhi/degree << " degrees \n"
616 << " ending phi angle : " << endPhi/ 596 << " ending phi angle : " << endPhi/degree << " degrees \n";
617 G4int i=0; 597 G4int i=0;
618 << 598
619 os << " number of RZ points: " << numCorn 599 os << " number of RZ points: " << numCorner << "\n"
620 << " RZ values (corners): \n 600 << " RZ values (corners): \n";
621 for (i=0; i<numCorner; i++) 601 for (i=0; i<numCorner; i++)
622 { 602 {
623 os << " " 603 os << " "
624 << corners[i].r << ", " << corners[i 604 << corners[i].r << ", " << corners[i].z << "\n";
625 } 605 }
626 os << "------------------------------------- 606 os << "-----------------------------------------------------------\n";
627 os.precision(oldprc); 607 os.precision(oldprc);
628 608
629 return os; 609 return os;
630 } 610 }
631 611
632 ////////////////////////////////////////////// << 612 // GetPointOnSurface
633 // 613 //
634 // Return volume << 614 G4ThreeVector G4GenericPolycone::GetPointOnSurface() const
635 <<
636 G4double G4GenericPolycone::GetCubicVolume() <<
637 { 615 {
638 if (fCubicVolume == 0.) << 616 return GetPointOnSurfaceGeneric();
639 { << 617
640 G4double total = 0.; <<
641 G4int nrz = GetNumRZCorner(); <<
642 G4PolyconeSideRZ a = GetCorner(nrz - 1); <<
643 for (G4int i=0; i<nrz; ++i) <<
644 { <<
645 G4PolyconeSideRZ b = GetCorner(i); <<
646 total += (b.r*b.r + b.r*a.r + a.r*a.r)*( <<
647 a = b; <<
648 } <<
649 fCubicVolume = std::abs(total)*(GetEndPhi( <<
650 } <<
651 return fCubicVolume; <<
652 } 618 }
653 619
654 ////////////////////////////////////////////// << 620 // CreatePolyhedron
655 // 621 //
656 // Return surface area << 622 G4Polyhedron* G4GenericPolycone::CreatePolyhedron() const
657 << 623 {
658 G4double G4GenericPolycone::GetSurfaceArea() << 624 // The following code prepares for:
659 { << 625 // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
660 if (fSurfaceArea == 0.) << 626 // const double xyz[][3],
661 { << 627 // const int faces_vec[][4])
662 // phi cut area << 628 // Here is an extract from the header file HepPolyhedron.h:
663 G4int nrz = GetNumRZCorner(); << 629 /**
664 G4double scut = 0.; << 630 * Creates user defined polyhedron.
665 if (IsOpen()) << 631 * This function allows to the user to define arbitrary polyhedron.
>> 632 * The faces of the polyhedron should be either triangles or planar
>> 633 * quadrilateral. Nodes of a face are defined by indexes pointing to
>> 634 * the elements in the xyz array. Numeration of the elements in the
>> 635 * array starts from 1 (like in fortran). The indexes can be positive
>> 636 * or negative. Negative sign means that the corresponding edge is
>> 637 * invisible. The normal of the face should be directed to exterior
>> 638 * of the polyhedron.
>> 639 *
>> 640 * @param Nnodes number of nodes
>> 641 * @param Nfaces number of faces
>> 642 * @param xyz nodes
>> 643 * @param faces_vec faces (quadrilaterals or triangles)
>> 644 * @return status of the operation - is non-zero in case of problem
>> 645 */
>> 646 const G4int numSide =
>> 647 G4int(G4Polyhedron::GetNumberOfRotationSteps()
>> 648 * (endPhi - startPhi) / twopi) + 1;
>> 649 G4int nNodes;
>> 650 G4int nFaces;
>> 651 typedef G4double double3[3];
>> 652 double3* xyz;
>> 653 typedef G4int int4[4];
>> 654 int4* faces_vec;
>> 655 if (phiIsOpen)
666 { 656 {
667 G4PolyconeSideRZ a = GetCorner(nrz - 1); << 657 // Triangulate open ends. Simple ear-chopping algorithm...
668 for (G4int i=0; i<nrz; ++i) << 658 // I'm not sure how robust this algorithm is (J.Allison).
>> 659 //
>> 660 std::vector<G4bool> chopped(numCorner, false);
>> 661 std::vector<G4int*> triQuads;
>> 662 G4int remaining = numCorner;
>> 663 G4int iStarter = 0;
>> 664 while (remaining >= 3) // Loop checking, 13.08.2015, G.Cosmo
669 { 665 {
670 G4PolyconeSideRZ b = GetCorner(i); << 666 // Find unchopped corners...
671 scut += a.r*b.z - a.z*b.r; << 667 //
672 a = b; << 668 G4int A = -1, B = -1, C = -1;
>> 669 G4int iStepper = iStarter;
>> 670 do // Loop checking, 13.08.2015, G.Cosmo
>> 671 {
>> 672 if (A < 0) { A = iStepper; }
>> 673 else if (B < 0) { B = iStepper; }
>> 674 else if (C < 0) { C = iStepper; }
>> 675 do // Loop checking, 13.08.2015, G.Cosmo
>> 676 {
>> 677 if (++iStepper >= numCorner) { iStepper = 0; }
>> 678 }
>> 679 while (chopped[iStepper]);
>> 680 }
>> 681 while (C < 0 && iStepper != iStarter);
>> 682
>> 683 // Check triangle at B is pointing outward (an "ear").
>> 684 // Sign of z cross product determines...
>> 685 //
>> 686 G4double BAr = corners[A].r - corners[B].r;
>> 687 G4double BAz = corners[A].z - corners[B].z;
>> 688 G4double BCr = corners[C].r - corners[B].r;
>> 689 G4double BCz = corners[C].z - corners[B].z;
>> 690 if (BAr * BCz - BAz * BCr < kCarTolerance)
>> 691 {
>> 692 G4int* tq = new G4int[3];
>> 693 tq[0] = A + 1;
>> 694 tq[1] = B + 1;
>> 695 tq[2] = C + 1;
>> 696 triQuads.push_back(tq);
>> 697 chopped[B] = true;
>> 698 --remaining;
>> 699 }
>> 700 else
>> 701 {
>> 702 do // Loop checking, 13.08.2015, G.Cosmo
>> 703 {
>> 704 if (++iStarter >= numCorner) { iStarter = 0; }
>> 705 }
>> 706 while (chopped[iStarter]);
>> 707 }
>> 708 }
>> 709 // Transfer to faces...
>> 710 //
>> 711 nNodes = (numSide + 1) * numCorner;
>> 712 nFaces = numSide * numCorner + 2 * triQuads.size();
>> 713 faces_vec = new int4[nFaces];
>> 714 G4int iface = 0;
>> 715 G4int addition = numCorner * numSide;
>> 716 G4int d = numCorner - 1;
>> 717 for (G4int iEnd = 0; iEnd < 2; ++iEnd)
>> 718 {
>> 719 for (size_t i = 0; i < triQuads.size(); ++i)
>> 720 {
>> 721 // Negative for soft/auxiliary/normally invisible edges...
>> 722 //
>> 723 G4int a, b, c;
>> 724 if (iEnd == 0)
>> 725 {
>> 726 a = triQuads[i][0];
>> 727 b = triQuads[i][1];
>> 728 c = triQuads[i][2];
>> 729 }
>> 730 else
>> 731 {
>> 732 a = triQuads[i][0] + addition;
>> 733 b = triQuads[i][2] + addition;
>> 734 c = triQuads[i][1] + addition;
>> 735 }
>> 736 G4int ab = std::abs(b - a);
>> 737 G4int bc = std::abs(c - b);
>> 738 G4int ca = std::abs(a - c);
>> 739 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
>> 740 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
>> 741 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
>> 742 faces_vec[iface][3] = 0;
>> 743 ++iface;
>> 744 }
673 } 745 }
674 scut = std::abs(scut); <<
675 } <<
676 // lateral surface area <<
677 G4double slat = 0; <<
678 G4PolyconeSideRZ a = GetCorner(nrz - 1); <<
679 for (G4int i=0; i<nrz; ++i) <<
680 { <<
681 G4PolyconeSideRZ b = GetCorner(i); <<
682 G4double h = std::sqrt((b.r - a.r)*(b.r <<
683 slat += (b.r + a.r)*h; <<
684 a = b; <<
685 } <<
686 slat *= (GetEndPhi() - GetStartPhi())/2.; <<
687 fSurfaceArea = scut + slat; <<
688 } <<
689 return fSurfaceArea; <<
690 } <<
691 746
692 ////////////////////////////////////////////// << 747 // Continue with sides...
693 // <<
694 // Set vector of surface elements, auxiliary m <<
695 // random points on surface <<
696 748
697 void G4GenericPolycone::SetSurfaceElements() c << 749 xyz = new double3[nNodes];
698 { << 750 const G4double dPhi = (endPhi - startPhi) / numSide;
699 fElements = new std::vector<G4GenericPolycon << 751 G4double phi = startPhi;
700 G4double sarea = 0.; << 752 G4int ixyz = 0;
701 G4int nrz = GetNumRZCorner(); << 753 for (G4int iSide = 0; iSide < numSide; ++iSide)
>> 754 {
>> 755 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
>> 756 {
>> 757 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
>> 758 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
>> 759 xyz[ixyz][2] = corners[iCorner].z;
>> 760 if (iSide == 0) // startPhi
>> 761 {
>> 762 if (iCorner < numCorner - 1)
>> 763 {
>> 764 faces_vec[iface][0] = ixyz + 1;
>> 765 faces_vec[iface][1] = -(ixyz + numCorner + 1);
>> 766 faces_vec[iface][2] = ixyz + numCorner + 2;
>> 767 faces_vec[iface][3] = ixyz + 2;
>> 768 }
>> 769 else
>> 770 {
>> 771 faces_vec[iface][0] = ixyz + 1;
>> 772 faces_vec[iface][1] = -(ixyz + numCorner + 1);
>> 773 faces_vec[iface][2] = ixyz + 2;
>> 774 faces_vec[iface][3] = ixyz - numCorner + 2;
>> 775 }
>> 776 }
>> 777 else if (iSide == numSide - 1) // endPhi
>> 778 {
>> 779 if (iCorner < numCorner - 1)
>> 780 {
>> 781 faces_vec[iface][0] = ixyz + 1;
>> 782 faces_vec[iface][1] = ixyz + numCorner + 1;
>> 783 faces_vec[iface][2] = ixyz + numCorner + 2;
>> 784 faces_vec[iface][3] = -(ixyz + 2);
>> 785 }
>> 786 else
>> 787 {
>> 788 faces_vec[iface][0] = ixyz + 1;
>> 789 faces_vec[iface][1] = ixyz + numCorner + 1;
>> 790 faces_vec[iface][2] = ixyz + 2;
>> 791 faces_vec[iface][3] = -(ixyz - numCorner + 2);
>> 792 }
>> 793 }
>> 794 else
>> 795 {
>> 796 if (iCorner < numCorner - 1)
>> 797 {
>> 798 faces_vec[iface][0] = ixyz + 1;
>> 799 faces_vec[iface][1] = -(ixyz + numCorner + 1);
>> 800 faces_vec[iface][2] = ixyz + numCorner + 2;
>> 801 faces_vec[iface][3] = -(ixyz + 2);
>> 802 }
>> 803 else
>> 804 {
>> 805 faces_vec[iface][0] = ixyz + 1;
>> 806 faces_vec[iface][1] = -(ixyz + numCorner + 1);
>> 807 faces_vec[iface][2] = ixyz + 2;
>> 808 faces_vec[iface][3] = -(ixyz - numCorner + 2);
>> 809 }
>> 810 }
>> 811 ++iface;
>> 812 ++ixyz;
>> 813 }
>> 814 phi += dPhi;
>> 815 }
702 816
703 // set lateral surface elements << 817 // Last corners...
704 G4double dphi = GetEndPhi() - GetStartPhi(); <<
705 G4int ia = nrz - 1; <<
706 for (G4int ib=0; ib<nrz; ++ib) <<
707 { <<
708 G4PolyconeSideRZ a = GetCorner(ia); <<
709 G4PolyconeSideRZ b = GetCorner(ib); <<
710 G4GenericPolycone::surface_element selem; <<
711 selem.i0 = ia; <<
712 selem.i1 = ib; <<
713 selem.i2 = -1; <<
714 ia = ib; <<
715 if (a.r == 0. && b.r == 0.) continue; <<
716 G4double h = std::sqrt((b.r - a.r)*(b.r - <<
717 sarea += 0.5*dphi*(b.r + a.r)*h; <<
718 selem.area = sarea; <<
719 fElements->push_back(selem); <<
720 } <<
721 818
722 // set elements for phi cuts << 819 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
723 if (IsOpen()) << 820 {
724 { << 821 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
725 G4TwoVectorList contourRZ; << 822 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
726 std::vector<G4int> triangles; << 823 xyz[ixyz][2] = corners[iCorner].z;
727 for (G4int i=0; i<nrz; ++i) << 824 ++ixyz;
728 { << 825 }
729 G4PolyconeSideRZ corner = GetCorner(i); <<
730 contourRZ.emplace_back(corner.r, corner. <<
731 } 826 }
732 G4GeomTools::TriangulatePolygon(contourRZ, << 827 else // !phiIsOpen - i.e., a complete 360 degrees.
733 auto ntria = (G4int)triangles.size(); <<
734 for (G4int i=0; i<ntria; i+=3) <<
735 { 828 {
736 G4GenericPolycone::surface_element selem << 829 nNodes = numSide * numCorner;
737 selem.i0 = triangles[i]; << 830 nFaces = numSide * numCorner;;
738 selem.i1 = triangles[i+1]; << 831 xyz = new double3[nNodes];
739 selem.i2 = triangles[i+2]; << 832 faces_vec = new int4[nFaces];
740 G4PolyconeSideRZ a = GetCorner(selem.i0) << 833 const G4double dPhi = (endPhi - startPhi) / numSide;
741 G4PolyconeSideRZ b = GetCorner(selem.i1) << 834 G4double phi = startPhi;
742 G4PolyconeSideRZ c = GetCorner(selem.i2) << 835 G4int ixyz = 0, iface = 0;
743 G4double stria = << 836 for (G4int iSide = 0; iSide < numSide; ++iSide)
744 std::abs(G4GeomTools::TriangleArea(a.r << 837 {
745 sarea += stria; << 838 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
746 selem.area = sarea; << 839 {
747 fElements->push_back(selem); // start ph << 840 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
748 sarea += stria; << 841 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
749 selem.area = sarea; << 842 xyz[ixyz][2] = corners[iCorner].z;
750 selem.i0 += nrz; << 843
751 fElements->push_back(selem); // end phi << 844 if (iSide < numSide - 1)
752 } << 845 {
753 } << 846 if (iCorner < numCorner - 1)
754 } << 847 {
755 << 848 faces_vec[iface][0] = ixyz + 1;
756 ////////////////////////////////////////////// << 849 faces_vec[iface][1] = -(ixyz + numCorner + 1);
757 // << 850 faces_vec[iface][2] = ixyz + numCorner + 2;
758 // Generate random point on surface << 851 faces_vec[iface][3] = -(ixyz + 2);
759 << 852 }
760 G4ThreeVector G4GenericPolycone::GetPointOnSur << 853 else
761 { << 854 {
762 // Set surface elements << 855 faces_vec[iface][0] = ixyz + 1;
763 if (fElements == nullptr) << 856 faces_vec[iface][1] = -(ixyz + numCorner + 1);
764 { << 857 faces_vec[iface][2] = ixyz + 2;
765 G4AutoLock l(&surface_elementsMutex); << 858 faces_vec[iface][3] = -(ixyz - numCorner + 2);
766 SetSurfaceElements(); << 859 }
767 l.unlock(); << 860 }
768 } << 861 else // Last side joins ends...
769 << 862 {
770 // Select surface element << 863 if (iCorner < numCorner - 1)
771 G4GenericPolycone::surface_element selem; << 864 {
772 selem = fElements->back(); << 865 faces_vec[iface][0] = ixyz + 1;
773 G4double select = selem.area*G4QuickRand(); << 866 faces_vec[iface][1] = -(ixyz + numCorner - nFaces + 1);
774 auto it = std::lower_bound(fElements->begin( << 867 faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
775 [](const G4Generi << 868 faces_vec[iface][3] = -(ixyz + 2);
776 -> G4bool { retur << 869 }
777 << 870 else
778 // Generate random point << 871 {
779 G4double r = 0, z = 0, phi = 0; << 872 faces_vec[iface][0] = ixyz + 1;
780 G4double u = G4QuickRand(); << 873 faces_vec[iface][1] = -(ixyz - nFaces + numCorner + 1);
781 G4double v = G4QuickRand(); << 874 faces_vec[iface][2] = ixyz - nFaces + 2;
782 G4int i0 = (*it).i0; << 875 faces_vec[iface][3] = -(ixyz - numCorner + 2);
783 G4int i1 = (*it).i1; << 876 }
784 G4int i2 = (*it).i2; << 877 }
785 if (i2 < 0) // lateral surface << 878 ++ixyz;
786 { << 879 ++iface;
787 G4PolyconeSideRZ p0 = GetCorner(i0); << 880 }
788 G4PolyconeSideRZ p1 = GetCorner(i1); << 881 phi += dPhi;
789 if (p1.r < p0.r) << 882 }
790 { <<
791 p0 = GetCorner(i1); <<
792 p1 = GetCorner(i0); <<
793 } 883 }
794 if (p1.r - p0.r < kCarTolerance) // cylind << 884 G4Polyhedron* polyhedron = new G4Polyhedron;
795 { << 885 G4int prob = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
796 r = (p1.r - p0.r)*u + p0.r; << 886 delete [] faces_vec;
797 z = (p1.z - p0.z)*u + p0.z; << 887 delete [] xyz;
>> 888 if (prob)
>> 889 {
>> 890 std::ostringstream message;
>> 891 message << "Problem creating G4Polyhedron for: " << GetName();
>> 892 G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002",
>> 893 JustWarning, message);
>> 894 delete polyhedron;
>> 895 return nullptr;
798 } 896 }
799 else // conical surface << 897 else
800 { 898 {
801 r = std::sqrt(p1.r*p1.r*u + p0.r*p0.r*(1 << 899 return polyhedron;
802 z = p0.z + (p1.z - p0.z)*(r - p0.r)/(p1. <<
803 } 900 }
804 phi = (GetEndPhi() - GetStartPhi())*v + Ge <<
805 } <<
806 else // phi cut <<
807 { <<
808 G4int nrz = GetNumRZCorner(); <<
809 phi = (i0 < nrz) ? GetStartPhi() : GetEndP <<
810 if (i0 >= nrz) { i0 -= nrz; } <<
811 G4PolyconeSideRZ p0 = GetCorner(i0); <<
812 G4PolyconeSideRZ p1 = GetCorner(i1); <<
813 G4PolyconeSideRZ p2 = GetCorner(i2); <<
814 if (u + v > 1.) { u = 1. - u; v = 1. - v; <<
815 r = (p1.r - p0.r)*u + (p2.r - p0.r)*v + p <<
816 z = (p1.z - p0.z)*u + (p2.z - p0.z)*v + p <<
817 } <<
818 return {r*std::cos(phi), r*std::sin(phi), z} <<
819 } <<
820 <<
821 ////////////////////////////////////////////// <<
822 // <<
823 // CreatePolyhedron <<
824 <<
825 G4Polyhedron* G4GenericPolycone::CreatePolyhed <<
826 { <<
827 std::vector<G4TwoVector> rz(numCorner); <<
828 for (G4int i = 0; i < numCorner; ++i) <<
829 rz[i].set(corners[i].r, corners[i].z); <<
830 return new G4PolyhedronPcon(startPhi, endPhi <<
831 } 901 }
832 902
833 #endif 903 #endif
834 904