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