| 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 const 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 <<
327 // <<
328 // Surface elements <<
329 // <<
330 delete fElements; <<
331 fElements = nullptr; <<
332 <<
333 // Polyhedron <<
334 // <<
335 fRebuildPolyhedron = false; <<
336 delete fpPolyhedron; <<
337 fpPolyhedron = nullptr; <<
338 } 328 }
339 329
>> 330
>> 331 //
340 // Reset 332 // Reset
341 // 333 //
342 G4bool G4GenericPolycone::Reset() 334 G4bool G4GenericPolycone::Reset()
343 { 335 {
344 std::ostringstream message; << 336
345 message << "Solid " << GetName() << " built << 337 std::ostringstream message;
346 << G4endl << "Not applicable to the << 338 message << "Solid " << GetName() << " built using generic construct."
347 G4Exception("G4GenericPolycone::Reset()", "G << 339 << G4endl << "Not applicable to the generic construct !";
348 JustWarning, message, "Parameter << 340 G4Exception("G4GenericPolycone::Reset()", "GeomSolids1001",
349 return true; << 341 JustWarning, message, "Parameters NOT resetted.");
>> 342 return 1;
>> 343
350 } 344 }
351 345
>> 346
>> 347 //
352 // Inside 348 // Inside
353 // 349 //
354 // This is an override of G4VCSGfaceted::Insid 350 // This is an override of G4VCSGfaceted::Inside, created in order
355 // to speed things up by first checking with G 351 // to speed things up by first checking with G4EnclosingCylinder.
356 // 352 //
357 EInside G4GenericPolycone::Inside( const G4Thr << 353 EInside G4GenericPolycone::Inside( const G4ThreeVector &p ) const
358 { 354 {
359 // 355 //
360 // Quick test 356 // Quick test
361 // 357 //
362 if (enclosingCylinder->MustBeOutside(p)) ret 358 if (enclosingCylinder->MustBeOutside(p)) return kOutside;
363 359
364 // 360 //
365 // Long answer 361 // Long answer
366 // 362 //
367 return G4VCSGfaceted::Inside(p); 363 return G4VCSGfaceted::Inside(p);
368 } 364 }
369 365
>> 366
>> 367 //
370 // DistanceToIn 368 // DistanceToIn
371 // 369 //
372 // This is an override of G4VCSGfaceted::Insid 370 // This is an override of G4VCSGfaceted::Inside, created in order
373 // to speed things up by first checking with G 371 // to speed things up by first checking with G4EnclosingCylinder.
374 // 372 //
375 G4double G4GenericPolycone::DistanceToIn( cons << 373 G4double G4GenericPolycone::DistanceToIn( const G4ThreeVector &p,
376 cons << 374 const G4ThreeVector &v ) const
377 { 375 {
378 // 376 //
379 // Quick test 377 // Quick test
380 // 378 //
381 if (enclosingCylinder->ShouldMiss(p,v)) 379 if (enclosingCylinder->ShouldMiss(p,v))
382 return kInfinity; 380 return kInfinity;
383 << 381
384 // 382 //
385 // Long answer 383 // Long answer
386 // 384 //
387 return G4VCSGfaceted::DistanceToIn( p, v ); 385 return G4VCSGfaceted::DistanceToIn( p, v );
388 } 386 }
389 387
>> 388
>> 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
398 // 398 //
399 // Get bounding box << 399 // ComputeDimensions
400 // 400 //
401 void << 401 /*void G4GenericPolycone::ComputeDimensions( G4VPVParameterisation* p,
402 G4GenericPolycone::BoundingLimits(G4ThreeVecto << 402 const G4int n,
403 G4ThreeVecto << 403 const G4VPhysicalVolume* pRep )
404 { 404 {
405 G4double rmin = kInfinity, rmax = -kInfinity << 405 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 } 406 }
447 << 407 */
448 // CalculateExtent <<
449 // <<
450 // Calculate extent under transform and specif <<
451 // 408 //
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 409 // GetEntityType
588 // 410 //
589 G4GeometryType G4GenericPolycone::GetEntityTy 411 G4GeometryType G4GenericPolycone::GetEntityType() const
590 { 412 {
591 return {"G4GenericPolycone"}; << 413 return G4String("G4GenericPolycone");
592 } 414 }
593 415
594 // Clone << 416
595 // 417 //
596 // Make a clone of the object 418 // Make a clone of the object
597 // 419 //
598 G4VSolid* G4GenericPolycone::Clone() const 420 G4VSolid* G4GenericPolycone::Clone() const
599 { 421 {
600 return new G4GenericPolycone(*this); 422 return new G4GenericPolycone(*this);
601 } 423 }
602 424
603 // StreamInfo <<
604 // 425 //
605 // Stream object contents to an output stream 426 // Stream object contents to an output stream
606 // 427 //
607 std::ostream& G4GenericPolycone::StreamInfo( s 428 std::ostream& G4GenericPolycone::StreamInfo( std::ostream& os ) const
608 { 429 {
609 G4long oldprc = os.precision(16); << 430 G4int oldprc = os.precision(16);
610 os << "------------------------------------- 431 os << "-----------------------------------------------------------\n"
611 << " *** Dump for solid - " << GetName 432 << " *** Dump for solid - " << GetName() << " ***\n"
612 << " ================================= 433 << " ===================================================\n"
613 << " Solid type: G4GenericPolycone\n" 434 << " Solid type: G4GenericPolycone\n"
614 << " Parameters: \n" 435 << " Parameters: \n"
615 << " starting phi angle : " << startPh 436 << " starting phi angle : " << startPhi/degree << " degrees \n"
616 << " ending phi angle : " << endPhi/ 437 << " ending phi angle : " << endPhi/degree << " degrees \n";
617 G4int i=0; 438 G4int i=0;
618 << 439
619 os << " number of RZ points: " << numCorn 440 os << " number of RZ points: " << numCorner << "\n"
620 << " RZ values (corners): \n 441 << " RZ values (corners): \n";
621 for (i=0; i<numCorner; i++) 442 for (i=0; i<numCorner; i++)
622 { 443 {
623 os << " " 444 os << " "
624 << corners[i].r << ", " << corners[i 445 << corners[i].r << ", " << corners[i].z << "\n";
625 } 446 }
626 os << "------------------------------------- 447 os << "-----------------------------------------------------------\n";
627 os.precision(oldprc); 448 os.precision(oldprc);
628 449
629 return os; 450 return os;
630 } 451 }
631 452
632 ////////////////////////////////////////////// <<
633 // <<
634 // Return volume <<
635 453
636 G4double G4GenericPolycone::GetCubicVolume() << 454
>> 455 //
>> 456 // GetPointOnSurface
>> 457 //
>> 458 G4ThreeVector G4GenericPolycone::GetPointOnSurface() const
637 { 459 {
638 if (fCubicVolume == 0.) << 460 return GetPointOnSurfaceGeneric();
639 { << 461
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 } 462 }
653 463
654 ////////////////////////////////////////////// <<
655 // 464 //
656 // Return surface area << 465 // CreatePolyhedron
657 << 466 //
658 G4double G4GenericPolycone::GetSurfaceArea() << 467 G4Polyhedron* G4GenericPolycone::CreatePolyhedron() const
659 { << 468 {
660 if (fSurfaceArea == 0.) << 469 // The following code prepares for:
661 { << 470 // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
662 // phi cut area << 471 // const double xyz[][3],
663 G4int nrz = GetNumRZCorner(); << 472 // const int faces_vec[][4])
664 G4double scut = 0.; << 473 // Here is an extract from the header file HepPolyhedron.h:
665 if (IsOpen()) << 474 /**
>> 475 * Creates user defined polyhedron.
>> 476 * This function allows to the user to define arbitrary polyhedron.
>> 477 * The faces of the polyhedron should be either triangles or planar
>> 478 * quadrilateral. Nodes of a face are defined by indexes pointing to
>> 479 * the elements in the xyz array. Numeration of the elements in the
>> 480 * array starts from 1 (like in fortran). The indexes can be positive
>> 481 * or negative. Negative sign means that the corresponding edge is
>> 482 * invisible. The normal of the face should be directed to exterior
>> 483 * of the polyhedron.
>> 484 *
>> 485 * @param Nnodes number of nodes
>> 486 * @param Nfaces number of faces
>> 487 * @param xyz nodes
>> 488 * @param faces_vec faces (quadrilaterals or triangles)
>> 489 * @return status of the operation - is non-zero in case of problem
>> 490 */
>> 491 const G4int numSide =
>> 492 G4int(G4Polyhedron::GetNumberOfRotationSteps()
>> 493 * (endPhi - startPhi) / twopi) + 1;
>> 494 G4int nNodes;
>> 495 G4int nFaces;
>> 496 typedef G4double double3[3];
>> 497 double3* xyz;
>> 498 typedef G4int int4[4];
>> 499 int4* faces_vec;
>> 500 if (phiIsOpen)
666 { 501 {
667 G4PolyconeSideRZ a = GetCorner(nrz - 1); << 502 // Triangulate open ends. Simple ear-chopping algorithm...
668 for (G4int i=0; i<nrz; ++i) << 503 // I'm not sure how robust this algorithm is (J.Allison).
>> 504 //
>> 505 std::vector<G4bool> chopped(numCorner, false);
>> 506 std::vector<G4int*> triQuads;
>> 507 G4int remaining = numCorner;
>> 508 G4int iStarter = 0;
>> 509 while (remaining >= 3)
669 { 510 {
670 G4PolyconeSideRZ b = GetCorner(i); << 511 // Find unchopped corners...
671 scut += a.r*b.z - a.z*b.r; << 512 //
672 a = b; << 513 G4int A = -1, B = -1, C = -1;
>> 514 G4int iStepper = iStarter;
>> 515 do
>> 516 {
>> 517 if (A < 0) { A = iStepper; }
>> 518 else if (B < 0) { B = iStepper; }
>> 519 else if (C < 0) { C = iStepper; }
>> 520 do
>> 521 {
>> 522 if (++iStepper >= numCorner) { iStepper = 0; }
>> 523 }
>> 524 while (chopped[iStepper]);
>> 525 }
>> 526 while (C < 0 && iStepper != iStarter);
>> 527
>> 528 // Check triangle at B is pointing outward (an "ear").
>> 529 // Sign of z cross product determines...
>> 530 //
>> 531 G4double BAr = corners[A].r - corners[B].r;
>> 532 G4double BAz = corners[A].z - corners[B].z;
>> 533 G4double BCr = corners[C].r - corners[B].r;
>> 534 G4double BCz = corners[C].z - corners[B].z;
>> 535 if (BAr * BCz - BAz * BCr < kCarTolerance)
>> 536 {
>> 537 G4int* tq = new G4int[3];
>> 538 tq[0] = A + 1;
>> 539 tq[1] = B + 1;
>> 540 tq[2] = C + 1;
>> 541 triQuads.push_back(tq);
>> 542 chopped[B] = true;
>> 543 --remaining;
>> 544 }
>> 545 else
>> 546 {
>> 547 do
>> 548 {
>> 549 if (++iStarter >= numCorner) { iStarter = 0; }
>> 550 }
>> 551 while (chopped[iStarter]);
>> 552 }
>> 553 }
>> 554 // Transfer to faces...
>> 555 //
>> 556 nNodes = (numSide + 1) * numCorner;
>> 557 nFaces = numSide * numCorner + 2 * triQuads.size();
>> 558 faces_vec = new int4[nFaces];
>> 559 G4int iface = 0;
>> 560 G4int addition = numCorner * numSide;
>> 561 G4int d = numCorner - 1;
>> 562 for (G4int iEnd = 0; iEnd < 2; ++iEnd)
>> 563 {
>> 564 for (size_t i = 0; i < triQuads.size(); ++i)
>> 565 {
>> 566 // Negative for soft/auxiliary/normally invisible edges...
>> 567 //
>> 568 G4int a, b, c;
>> 569 if (iEnd == 0)
>> 570 {
>> 571 a = triQuads[i][0];
>> 572 b = triQuads[i][1];
>> 573 c = triQuads[i][2];
>> 574 }
>> 575 else
>> 576 {
>> 577 a = triQuads[i][0] + addition;
>> 578 b = triQuads[i][2] + addition;
>> 579 c = triQuads[i][1] + addition;
>> 580 }
>> 581 G4int ab = std::abs(b - a);
>> 582 G4int bc = std::abs(c - b);
>> 583 G4int ca = std::abs(a - c);
>> 584 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
>> 585 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
>> 586 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
>> 587 faces_vec[iface][3] = 0;
>> 588 ++iface;
>> 589 }
673 } 590 }
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 591
692 ////////////////////////////////////////////// << 592 // Continue with sides...
693 // <<
694 // Set vector of surface elements, auxiliary m <<
695 // random points on surface <<
696 593
697 void G4GenericPolycone::SetSurfaceElements() c << 594 xyz = new double3[nNodes];
698 { << 595 const G4double dPhi = (endPhi - startPhi) / numSide;
699 fElements = new std::vector<G4GenericPolycon << 596 G4double phi = startPhi;
700 G4double sarea = 0.; << 597 G4int ixyz = 0;
701 G4int nrz = GetNumRZCorner(); << 598 for (G4int iSide = 0; iSide < numSide; ++iSide)
>> 599 {
>> 600 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
>> 601 {
>> 602 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
>> 603 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
>> 604 xyz[ixyz][2] = corners[iCorner].z;
>> 605 if (iSide == 0) // startPhi
>> 606 {
>> 607 if (iCorner < numCorner - 1)
>> 608 {
>> 609 faces_vec[iface][0] = ixyz + 1;
>> 610 faces_vec[iface][1] = -(ixyz + numCorner + 1);
>> 611 faces_vec[iface][2] = ixyz + numCorner + 2;
>> 612 faces_vec[iface][3] = ixyz + 2;
>> 613 }
>> 614 else
>> 615 {
>> 616 faces_vec[iface][0] = ixyz + 1;
>> 617 faces_vec[iface][1] = -(ixyz + numCorner + 1);
>> 618 faces_vec[iface][2] = ixyz + 2;
>> 619 faces_vec[iface][3] = ixyz - numCorner + 2;
>> 620 }
>> 621 }
>> 622 else if (iSide == numSide - 1) // endPhi
>> 623 {
>> 624 if (iCorner < numCorner - 1)
>> 625 {
>> 626 faces_vec[iface][0] = ixyz + 1;
>> 627 faces_vec[iface][1] = ixyz + numCorner + 1;
>> 628 faces_vec[iface][2] = ixyz + numCorner + 2;
>> 629 faces_vec[iface][3] = -(ixyz + 2);
>> 630 }
>> 631 else
>> 632 {
>> 633 faces_vec[iface][0] = ixyz + 1;
>> 634 faces_vec[iface][1] = ixyz + numCorner + 1;
>> 635 faces_vec[iface][2] = ixyz + 2;
>> 636 faces_vec[iface][3] = -(ixyz - numCorner + 2);
>> 637 }
>> 638 }
>> 639 else
>> 640 {
>> 641 if (iCorner < numCorner - 1)
>> 642 {
>> 643 faces_vec[iface][0] = ixyz + 1;
>> 644 faces_vec[iface][1] = -(ixyz + numCorner + 1);
>> 645 faces_vec[iface][2] = ixyz + numCorner + 2;
>> 646 faces_vec[iface][3] = -(ixyz + 2);
>> 647 }
>> 648 else
>> 649 {
>> 650 faces_vec[iface][0] = ixyz + 1;
>> 651 faces_vec[iface][1] = -(ixyz + numCorner + 1);
>> 652 faces_vec[iface][2] = ixyz + 2;
>> 653 faces_vec[iface][3] = -(ixyz - numCorner + 2);
>> 654 }
>> 655 }
>> 656 ++iface;
>> 657 ++ixyz;
>> 658 }
>> 659 phi += dPhi;
>> 660 }
702 661
703 // set lateral surface elements << 662 // 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 663
722 // set elements for phi cuts << 664 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
723 if (IsOpen()) << 665 {
724 { << 666 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
725 G4TwoVectorList contourRZ; << 667 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
726 std::vector<G4int> triangles; << 668 xyz[ixyz][2] = corners[iCorner].z;
727 for (G4int i=0; i<nrz; ++i) << 669 ++ixyz;
728 { << 670 }
729 G4PolyconeSideRZ corner = GetCorner(i); <<
730 contourRZ.emplace_back(corner.r, corner. <<
731 } 671 }
732 G4GeomTools::TriangulatePolygon(contourRZ, << 672 else // !phiIsOpen - i.e., a complete 360 degrees.
733 auto ntria = (G4int)triangles.size(); <<
734 for (G4int i=0; i<ntria; i+=3) <<
735 { 673 {
736 G4GenericPolycone::surface_element selem << 674 nNodes = numSide * numCorner;
737 selem.i0 = triangles[i]; << 675 nFaces = numSide * numCorner;;
738 selem.i1 = triangles[i+1]; << 676 xyz = new double3[nNodes];
739 selem.i2 = triangles[i+2]; << 677 faces_vec = new int4[nFaces];
740 G4PolyconeSideRZ a = GetCorner(selem.i0) << 678 const G4double dPhi = (endPhi - startPhi) / numSide;
741 G4PolyconeSideRZ b = GetCorner(selem.i1) << 679 G4double phi = startPhi;
742 G4PolyconeSideRZ c = GetCorner(selem.i2) << 680 G4int ixyz = 0, iface = 0;
743 G4double stria = << 681 for (G4int iSide = 0; iSide < numSide; ++iSide)
744 std::abs(G4GeomTools::TriangleArea(a.r << 682 {
745 sarea += stria; << 683 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
746 selem.area = sarea; << 684 {
747 fElements->push_back(selem); // start ph << 685 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
748 sarea += stria; << 686 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
749 selem.area = sarea; << 687 xyz[ixyz][2] = corners[iCorner].z;
750 selem.i0 += nrz; << 688
751 fElements->push_back(selem); // end phi << 689 if (iSide < numSide - 1)
752 } << 690 {
753 } << 691 if (iCorner < numCorner - 1)
754 } << 692 {
755 << 693 faces_vec[iface][0] = ixyz + 1;
756 ////////////////////////////////////////////// << 694 faces_vec[iface][1] = -(ixyz + numCorner + 1);
757 // << 695 faces_vec[iface][2] = ixyz + numCorner + 2;
758 // Generate random point on surface << 696 faces_vec[iface][3] = -(ixyz + 2);
759 << 697 }
760 G4ThreeVector G4GenericPolycone::GetPointOnSur << 698 else
761 { << 699 {
762 // Set surface elements << 700 faces_vec[iface][0] = ixyz + 1;
763 if (fElements == nullptr) << 701 faces_vec[iface][1] = -(ixyz + numCorner + 1);
764 { << 702 faces_vec[iface][2] = ixyz + 2;
765 G4AutoLock l(&surface_elementsMutex); << 703 faces_vec[iface][3] = -(ixyz - numCorner + 2);
766 SetSurfaceElements(); << 704 }
767 l.unlock(); << 705 }
768 } << 706 else // Last side joins ends...
769 << 707 {
770 // Select surface element << 708 if (iCorner < numCorner - 1)
771 G4GenericPolycone::surface_element selem; << 709 {
772 selem = fElements->back(); << 710 faces_vec[iface][0] = ixyz + 1;
773 G4double select = selem.area*G4QuickRand(); << 711 faces_vec[iface][1] = -(ixyz + numCorner - nFaces + 1);
774 auto it = std::lower_bound(fElements->begin( << 712 faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
775 [](const G4Generi << 713 faces_vec[iface][3] = -(ixyz + 2);
776 -> G4bool { retur << 714 }
777 << 715 else
778 // Generate random point << 716 {
779 G4double r = 0, z = 0, phi = 0; << 717 faces_vec[iface][0] = ixyz + 1;
780 G4double u = G4QuickRand(); << 718 faces_vec[iface][1] = -(ixyz - nFaces + numCorner + 1);
781 G4double v = G4QuickRand(); << 719 faces_vec[iface][2] = ixyz - nFaces + 2;
782 G4int i0 = (*it).i0; << 720 faces_vec[iface][3] = -(ixyz - numCorner + 2);
783 G4int i1 = (*it).i1; << 721 }
784 G4int i2 = (*it).i2; << 722 }
785 if (i2 < 0) // lateral surface << 723 ++ixyz;
786 { << 724 ++iface;
787 G4PolyconeSideRZ p0 = GetCorner(i0); << 725 }
788 G4PolyconeSideRZ p1 = GetCorner(i1); << 726 phi += dPhi;
789 if (p1.r < p0.r) << 727 }
790 { <<
791 p0 = GetCorner(i1); <<
792 p1 = GetCorner(i0); <<
793 } 728 }
794 if (p1.r - p0.r < kCarTolerance) // cylind << 729 G4Polyhedron* polyhedron = new G4Polyhedron;
795 { << 730 G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
796 r = (p1.r - p0.r)*u + p0.r; << 731 delete [] faces_vec;
797 z = (p1.z - p0.z)*u + p0.z; << 732 delete [] xyz;
>> 733 if (problem)
>> 734 {
>> 735 std::ostringstream message;
>> 736 message << "Problem creating G4Polyhedron for: " << GetName();
>> 737 G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002",
>> 738 JustWarning, message);
>> 739 delete polyhedron;
>> 740 return 0;
798 } 741 }
799 else // conical surface << 742 else
800 { 743 {
801 r = std::sqrt(p1.r*p1.r*u + p0.r*p0.r*(1 << 744 return polyhedron;
802 z = p0.z + (p1.z - p0.z)*(r - p0.r)/(p1. <<
803 } 745 }
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 } 746 }
832 747
833 #endif 748 #endif
834 749