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