1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: G4Polycone.cc,v 1.44 2010/06/16 08:24:14 gcosmo Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-04-beta-01 $ |
---|
29 | // |
---|
30 | // |
---|
31 | // -------------------------------------------------------------------- |
---|
32 | // GEANT 4 class source file |
---|
33 | // |
---|
34 | // |
---|
35 | // G4Polycone.cc |
---|
36 | // |
---|
37 | // Implementation of a CSG polycone |
---|
38 | // |
---|
39 | // -------------------------------------------------------------------- |
---|
40 | |
---|
41 | #include "G4Polycone.hh" |
---|
42 | |
---|
43 | #include "G4PolyconeSide.hh" |
---|
44 | #include "G4PolyPhiFace.hh" |
---|
45 | |
---|
46 | #include "Randomize.hh" |
---|
47 | |
---|
48 | #include "G4Polyhedron.hh" |
---|
49 | #include "G4EnclosingCylinder.hh" |
---|
50 | #include "G4ReduciblePolygon.hh" |
---|
51 | #include "G4VPVParameterisation.hh" |
---|
52 | |
---|
53 | using namespace CLHEP; |
---|
54 | |
---|
55 | // |
---|
56 | // Constructor (GEANT3 style parameters) |
---|
57 | // |
---|
58 | G4Polycone::G4Polycone( const G4String& name, |
---|
59 | G4double phiStart, |
---|
60 | G4double phiTotal, |
---|
61 | G4int numZPlanes, |
---|
62 | const G4double zPlane[], |
---|
63 | const G4double rInner[], |
---|
64 | const G4double rOuter[] ) |
---|
65 | : G4VCSGfaceted( name ), genericPcon(false) |
---|
66 | { |
---|
67 | // |
---|
68 | // Some historical ugliness |
---|
69 | // |
---|
70 | original_parameters = new G4PolyconeHistorical(); |
---|
71 | |
---|
72 | original_parameters->Start_angle = phiStart; |
---|
73 | original_parameters->Opening_angle = phiTotal; |
---|
74 | original_parameters->Num_z_planes = numZPlanes; |
---|
75 | original_parameters->Z_values = new G4double[numZPlanes]; |
---|
76 | original_parameters->Rmin = new G4double[numZPlanes]; |
---|
77 | original_parameters->Rmax = new G4double[numZPlanes]; |
---|
78 | |
---|
79 | G4int i; |
---|
80 | for (i=0; i<numZPlanes; i++) |
---|
81 | { |
---|
82 | if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] )) |
---|
83 | { |
---|
84 | if( (rInner[i] > rOuter[i+1]) |
---|
85 | ||(rInner[i+1] > rOuter[i]) ) |
---|
86 | { |
---|
87 | DumpInfo(); |
---|
88 | G4cerr << "ERROR - G4Polycone::G4Polycone()" |
---|
89 | << G4endl |
---|
90 | << " Segments are not contiguous !" << G4endl |
---|
91 | << " rMin[" << i << "] = " << rInner[i] |
---|
92 | << " -- rMax[" << i+1 << "] = " << rOuter[i+1] << G4endl |
---|
93 | << " rMin[" << i+1 << "] = " << rInner[i+1] |
---|
94 | << " -- rMax[" << i << "] = " << rOuter[i] << G4endl; |
---|
95 | G4Exception("G4Polycone::G4Polycone()", "InvalidSetup", FatalException, |
---|
96 | "Cannot create a Polycone with no contiguous segments."); |
---|
97 | } |
---|
98 | } |
---|
99 | original_parameters->Z_values[i] = zPlane[i]; |
---|
100 | original_parameters->Rmin[i] = rInner[i]; |
---|
101 | original_parameters->Rmax[i] = rOuter[i]; |
---|
102 | } |
---|
103 | |
---|
104 | // |
---|
105 | // Build RZ polygon using special PCON/PGON GEANT3 constructor |
---|
106 | // |
---|
107 | G4ReduciblePolygon *rz = |
---|
108 | new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes ); |
---|
109 | |
---|
110 | // |
---|
111 | // Do the real work |
---|
112 | // |
---|
113 | Create( phiStart, phiTotal, rz ); |
---|
114 | |
---|
115 | delete rz; |
---|
116 | } |
---|
117 | |
---|
118 | |
---|
119 | // |
---|
120 | // Constructor (generic parameters) |
---|
121 | // |
---|
122 | G4Polycone::G4Polycone( const G4String& name, |
---|
123 | G4double phiStart, |
---|
124 | G4double phiTotal, |
---|
125 | G4int numRZ, |
---|
126 | const G4double r[], |
---|
127 | const G4double z[] ) |
---|
128 | : G4VCSGfaceted( name ), genericPcon(true) |
---|
129 | { |
---|
130 | G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ ); |
---|
131 | |
---|
132 | Create( phiStart, phiTotal, rz ); |
---|
133 | |
---|
134 | // Set original_parameters struct for consistency |
---|
135 | // |
---|
136 | SetOriginalParameters(); |
---|
137 | |
---|
138 | delete rz; |
---|
139 | } |
---|
140 | |
---|
141 | |
---|
142 | // |
---|
143 | // Create |
---|
144 | // |
---|
145 | // Generic create routine, called by each constructor after |
---|
146 | // conversion of arguments |
---|
147 | // |
---|
148 | void G4Polycone::Create( G4double phiStart, |
---|
149 | G4double phiTotal, |
---|
150 | G4ReduciblePolygon *rz ) |
---|
151 | { |
---|
152 | // |
---|
153 | // Perform checks of rz values |
---|
154 | // |
---|
155 | if (rz->Amin() < 0.0) |
---|
156 | { |
---|
157 | G4cerr << "ERROR - G4Polycone::Create(): " << GetName() << G4endl |
---|
158 | << " All R values must be >= 0 !" |
---|
159 | << G4endl; |
---|
160 | G4Exception("G4Polycone::Create()", "InvalidSetup", FatalException, |
---|
161 | "Illegal input parameters."); |
---|
162 | } |
---|
163 | |
---|
164 | G4double rzArea = rz->Area(); |
---|
165 | if (rzArea < -kCarTolerance) |
---|
166 | rz->ReverseOrder(); |
---|
167 | |
---|
168 | else if (rzArea < -kCarTolerance) |
---|
169 | { |
---|
170 | G4cerr << "ERROR - G4Polycone::Create(): " << GetName() << G4endl |
---|
171 | << " R/Z cross section is zero or near zero: " |
---|
172 | << rzArea << G4endl; |
---|
173 | G4Exception("G4Polycone::Create()", "InvalidSetup", FatalException, |
---|
174 | "Illegal input parameters."); |
---|
175 | } |
---|
176 | |
---|
177 | if ( (!rz->RemoveDuplicateVertices( kCarTolerance )) |
---|
178 | || (!rz->RemoveRedundantVertices( kCarTolerance )) ) |
---|
179 | { |
---|
180 | G4cerr << "ERROR - G4Polycone::Create(): " << GetName() << G4endl |
---|
181 | << " Too few unique R/Z values !" |
---|
182 | << G4endl; |
---|
183 | G4Exception("G4Polycone::Create()", "InvalidSetup", FatalException, |
---|
184 | "Illegal input parameters."); |
---|
185 | } |
---|
186 | |
---|
187 | if (rz->CrossesItself(1/kInfinity)) |
---|
188 | { |
---|
189 | G4cerr << "ERROR - G4Polycone::Create(): " << GetName() << G4endl |
---|
190 | << " R/Z segments cross !" |
---|
191 | << G4endl; |
---|
192 | G4Exception("G4Polycone::Create()", "InvalidSetup", FatalException, |
---|
193 | "Illegal input parameters."); |
---|
194 | } |
---|
195 | |
---|
196 | numCorner = rz->NumVertices(); |
---|
197 | |
---|
198 | // |
---|
199 | // Phi opening? Account for some possible roundoff, and interpret |
---|
200 | // nonsense value as representing no phi opening |
---|
201 | // |
---|
202 | if (phiTotal <= 0 || phiTotal > twopi-1E-10) |
---|
203 | { |
---|
204 | phiIsOpen = false; |
---|
205 | startPhi = 0; |
---|
206 | endPhi = twopi; |
---|
207 | } |
---|
208 | else |
---|
209 | { |
---|
210 | phiIsOpen = true; |
---|
211 | |
---|
212 | // |
---|
213 | // Convert phi into our convention |
---|
214 | // |
---|
215 | startPhi = phiStart; |
---|
216 | while( startPhi < 0 ) startPhi += twopi; |
---|
217 | |
---|
218 | endPhi = phiStart+phiTotal; |
---|
219 | while( endPhi < startPhi ) endPhi += twopi; |
---|
220 | } |
---|
221 | |
---|
222 | // |
---|
223 | // Allocate corner array. |
---|
224 | // |
---|
225 | corners = new G4PolyconeSideRZ[numCorner]; |
---|
226 | |
---|
227 | // |
---|
228 | // Copy corners |
---|
229 | // |
---|
230 | G4ReduciblePolygonIterator iterRZ(rz); |
---|
231 | |
---|
232 | G4PolyconeSideRZ *next = corners; |
---|
233 | iterRZ.Begin(); |
---|
234 | do |
---|
235 | { |
---|
236 | next->r = iterRZ.GetA(); |
---|
237 | next->z = iterRZ.GetB(); |
---|
238 | } while( ++next, iterRZ.Next() ); |
---|
239 | |
---|
240 | // |
---|
241 | // Allocate face pointer array |
---|
242 | // |
---|
243 | numFace = phiIsOpen ? numCorner+2 : numCorner; |
---|
244 | faces = new G4VCSGface*[numFace]; |
---|
245 | |
---|
246 | // |
---|
247 | // Construct conical faces |
---|
248 | // |
---|
249 | // But! Don't construct a face if both points are at zero radius! |
---|
250 | // |
---|
251 | G4PolyconeSideRZ *corner = corners, |
---|
252 | *prev = corners + numCorner-1, |
---|
253 | *nextNext; |
---|
254 | G4VCSGface **face = faces; |
---|
255 | do |
---|
256 | { |
---|
257 | next = corner+1; |
---|
258 | if (next >= corners+numCorner) next = corners; |
---|
259 | nextNext = next+1; |
---|
260 | if (nextNext >= corners+numCorner) nextNext = corners; |
---|
261 | |
---|
262 | if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue; |
---|
263 | |
---|
264 | // |
---|
265 | // We must decide here if we can dare declare one of our faces |
---|
266 | // as having a "valid" normal (i.e. allBehind = true). This |
---|
267 | // is never possible if the face faces "inward" in r. |
---|
268 | // |
---|
269 | G4bool allBehind; |
---|
270 | if (corner->z > next->z) |
---|
271 | { |
---|
272 | allBehind = false; |
---|
273 | } |
---|
274 | else |
---|
275 | { |
---|
276 | // |
---|
277 | // Otherwise, it is only true if the line passing |
---|
278 | // through the two points of the segment do not |
---|
279 | // split the r/z cross section |
---|
280 | // |
---|
281 | allBehind = !rz->BisectedBy( corner->r, corner->z, |
---|
282 | next->r, next->z, kCarTolerance ); |
---|
283 | } |
---|
284 | |
---|
285 | *face++ = new G4PolyconeSide( prev, corner, next, nextNext, |
---|
286 | startPhi, endPhi-startPhi, phiIsOpen, allBehind ); |
---|
287 | } while( prev=corner, corner=next, corner > corners ); |
---|
288 | |
---|
289 | if (phiIsOpen) |
---|
290 | { |
---|
291 | // |
---|
292 | // Construct phi open edges |
---|
293 | // |
---|
294 | *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi ); |
---|
295 | *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi ); |
---|
296 | } |
---|
297 | |
---|
298 | // |
---|
299 | // We might have dropped a face or two: recalculate numFace |
---|
300 | // |
---|
301 | numFace = face-faces; |
---|
302 | |
---|
303 | // |
---|
304 | // Make enclosingCylinder |
---|
305 | // |
---|
306 | enclosingCylinder = |
---|
307 | new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal ); |
---|
308 | } |
---|
309 | |
---|
310 | |
---|
311 | // |
---|
312 | // Fake default constructor - sets only member data and allocates memory |
---|
313 | // for usage restricted to object persistency. |
---|
314 | // |
---|
315 | G4Polycone::G4Polycone( __void__& a ) |
---|
316 | : G4VCSGfaceted(a), genericPcon(false), corners(0), |
---|
317 | original_parameters(0), enclosingCylinder(0) |
---|
318 | { |
---|
319 | } |
---|
320 | |
---|
321 | |
---|
322 | // |
---|
323 | // Destructor |
---|
324 | // |
---|
325 | G4Polycone::~G4Polycone() |
---|
326 | { |
---|
327 | delete [] corners; |
---|
328 | |
---|
329 | if (original_parameters) delete original_parameters; |
---|
330 | if (enclosingCylinder) delete enclosingCylinder; |
---|
331 | } |
---|
332 | |
---|
333 | |
---|
334 | // |
---|
335 | // Copy constructor |
---|
336 | // |
---|
337 | G4Polycone::G4Polycone( const G4Polycone &source ) |
---|
338 | : G4VCSGfaceted( source ) |
---|
339 | { |
---|
340 | CopyStuff( source ); |
---|
341 | } |
---|
342 | |
---|
343 | |
---|
344 | // |
---|
345 | // Assignment operator |
---|
346 | // |
---|
347 | const G4Polycone &G4Polycone::operator=( const G4Polycone &source ) |
---|
348 | { |
---|
349 | if (this == &source) return *this; |
---|
350 | |
---|
351 | G4VCSGfaceted::operator=( source ); |
---|
352 | |
---|
353 | delete [] corners; |
---|
354 | if (original_parameters) delete original_parameters; |
---|
355 | |
---|
356 | delete enclosingCylinder; |
---|
357 | |
---|
358 | CopyStuff( source ); |
---|
359 | |
---|
360 | return *this; |
---|
361 | } |
---|
362 | |
---|
363 | |
---|
364 | // |
---|
365 | // CopyStuff |
---|
366 | // |
---|
367 | void G4Polycone::CopyStuff( const G4Polycone &source ) |
---|
368 | { |
---|
369 | // |
---|
370 | // Simple stuff |
---|
371 | // |
---|
372 | startPhi = source.startPhi; |
---|
373 | endPhi = source.endPhi; |
---|
374 | phiIsOpen = source.phiIsOpen; |
---|
375 | numCorner = source.numCorner; |
---|
376 | genericPcon= source.genericPcon; |
---|
377 | |
---|
378 | // |
---|
379 | // The corner array |
---|
380 | // |
---|
381 | corners = new G4PolyconeSideRZ[numCorner]; |
---|
382 | |
---|
383 | G4PolyconeSideRZ *corn = corners, |
---|
384 | *sourceCorn = source.corners; |
---|
385 | do |
---|
386 | { |
---|
387 | *corn = *sourceCorn; |
---|
388 | } while( ++sourceCorn, ++corn < corners+numCorner ); |
---|
389 | |
---|
390 | // |
---|
391 | // Original parameters |
---|
392 | // |
---|
393 | if (source.original_parameters) |
---|
394 | { |
---|
395 | original_parameters = |
---|
396 | new G4PolyconeHistorical( *source.original_parameters ); |
---|
397 | } |
---|
398 | |
---|
399 | // |
---|
400 | // Enclosing cylinder |
---|
401 | // |
---|
402 | enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder ); |
---|
403 | } |
---|
404 | |
---|
405 | |
---|
406 | // |
---|
407 | // Reset |
---|
408 | // |
---|
409 | G4bool G4Polycone::Reset() |
---|
410 | { |
---|
411 | if (genericPcon) |
---|
412 | { |
---|
413 | G4cerr << "Solid " << GetName() << " built using generic construct." |
---|
414 | << G4endl << "Not applicable to the generic construct !" << G4endl; |
---|
415 | G4Exception("G4Polycone::Reset()", "NotApplicableConstruct", |
---|
416 | JustWarning, "Parameters NOT resetted."); |
---|
417 | return 1; |
---|
418 | } |
---|
419 | |
---|
420 | // |
---|
421 | // Clear old setup |
---|
422 | // |
---|
423 | G4VCSGfaceted::DeleteStuff(); |
---|
424 | delete [] corners; |
---|
425 | delete enclosingCylinder; |
---|
426 | |
---|
427 | // |
---|
428 | // Rebuild polycone |
---|
429 | // |
---|
430 | G4ReduciblePolygon *rz = |
---|
431 | new G4ReduciblePolygon( original_parameters->Rmin, |
---|
432 | original_parameters->Rmax, |
---|
433 | original_parameters->Z_values, |
---|
434 | original_parameters->Num_z_planes ); |
---|
435 | Create( original_parameters->Start_angle, |
---|
436 | original_parameters->Opening_angle, rz ); |
---|
437 | delete rz; |
---|
438 | |
---|
439 | return 0; |
---|
440 | } |
---|
441 | |
---|
442 | |
---|
443 | // |
---|
444 | // Inside |
---|
445 | // |
---|
446 | // This is an override of G4VCSGfaceted::Inside, created in order |
---|
447 | // to speed things up by first checking with G4EnclosingCylinder. |
---|
448 | // |
---|
449 | EInside G4Polycone::Inside( const G4ThreeVector &p ) const |
---|
450 | { |
---|
451 | // |
---|
452 | // Quick test |
---|
453 | // |
---|
454 | if (enclosingCylinder->MustBeOutside(p)) return kOutside; |
---|
455 | |
---|
456 | // |
---|
457 | // Long answer |
---|
458 | // |
---|
459 | return G4VCSGfaceted::Inside(p); |
---|
460 | } |
---|
461 | |
---|
462 | |
---|
463 | // |
---|
464 | // DistanceToIn |
---|
465 | // |
---|
466 | // This is an override of G4VCSGfaceted::Inside, created in order |
---|
467 | // to speed things up by first checking with G4EnclosingCylinder. |
---|
468 | // |
---|
469 | G4double G4Polycone::DistanceToIn( const G4ThreeVector &p, |
---|
470 | const G4ThreeVector &v ) const |
---|
471 | { |
---|
472 | // |
---|
473 | // Quick test |
---|
474 | // |
---|
475 | if (enclosingCylinder->ShouldMiss(p,v)) |
---|
476 | return kInfinity; |
---|
477 | |
---|
478 | // |
---|
479 | // Long answer |
---|
480 | // |
---|
481 | return G4VCSGfaceted::DistanceToIn( p, v ); |
---|
482 | } |
---|
483 | |
---|
484 | |
---|
485 | // |
---|
486 | // DistanceToIn |
---|
487 | // |
---|
488 | G4double G4Polycone::DistanceToIn( const G4ThreeVector &p ) const |
---|
489 | { |
---|
490 | return G4VCSGfaceted::DistanceToIn(p); |
---|
491 | } |
---|
492 | |
---|
493 | |
---|
494 | // |
---|
495 | // ComputeDimensions |
---|
496 | // |
---|
497 | void G4Polycone::ComputeDimensions( G4VPVParameterisation* p, |
---|
498 | const G4int n, |
---|
499 | const G4VPhysicalVolume* pRep ) |
---|
500 | { |
---|
501 | p->ComputeDimensions(*this,n,pRep); |
---|
502 | } |
---|
503 | |
---|
504 | // |
---|
505 | // GetEntityType |
---|
506 | // |
---|
507 | G4GeometryType G4Polycone::GetEntityType() const |
---|
508 | { |
---|
509 | return G4String("G4Polycone"); |
---|
510 | } |
---|
511 | |
---|
512 | // |
---|
513 | // Stream object contents to an output stream |
---|
514 | // |
---|
515 | std::ostream& G4Polycone::StreamInfo( std::ostream& os ) const |
---|
516 | { |
---|
517 | os << "-----------------------------------------------------------\n" |
---|
518 | << " *** Dump for solid - " << GetName() << " ***\n" |
---|
519 | << " ===================================================\n" |
---|
520 | << " Solid type: G4Polycone\n" |
---|
521 | << " Parameters: \n" |
---|
522 | << " starting phi angle : " << startPhi/degree << " degrees \n" |
---|
523 | << " ending phi angle : " << endPhi/degree << " degrees \n"; |
---|
524 | G4int i=0; |
---|
525 | if (!genericPcon) |
---|
526 | { |
---|
527 | G4int numPlanes = original_parameters->Num_z_planes; |
---|
528 | os << " number of Z planes: " << numPlanes << "\n" |
---|
529 | << " Z values: \n"; |
---|
530 | for (i=0; i<numPlanes; i++) |
---|
531 | { |
---|
532 | os << " Z plane " << i << ": " |
---|
533 | << original_parameters->Z_values[i] << "\n"; |
---|
534 | } |
---|
535 | os << " Tangent distances to inner surface (Rmin): \n"; |
---|
536 | for (i=0; i<numPlanes; i++) |
---|
537 | { |
---|
538 | os << " Z plane " << i << ": " |
---|
539 | << original_parameters->Rmin[i] << "\n"; |
---|
540 | } |
---|
541 | os << " Tangent distances to outer surface (Rmax): \n"; |
---|
542 | for (i=0; i<numPlanes; i++) |
---|
543 | { |
---|
544 | os << " Z plane " << i << ": " |
---|
545 | << original_parameters->Rmax[i] << "\n"; |
---|
546 | } |
---|
547 | } |
---|
548 | os << " number of RZ points: " << numCorner << "\n" |
---|
549 | << " RZ values (corners): \n"; |
---|
550 | for (i=0; i<numCorner; i++) |
---|
551 | { |
---|
552 | os << " " |
---|
553 | << corners[i].r << ", " << corners[i].z << "\n"; |
---|
554 | } |
---|
555 | os << "-----------------------------------------------------------\n"; |
---|
556 | |
---|
557 | return os; |
---|
558 | } |
---|
559 | |
---|
560 | |
---|
561 | // |
---|
562 | // GetPointOnCone |
---|
563 | // |
---|
564 | // Auxiliary method for Get Point On Surface |
---|
565 | // |
---|
566 | G4ThreeVector G4Polycone::GetPointOnCone(G4double fRmin1, G4double fRmax1, |
---|
567 | G4double fRmin2, G4double fRmax2, |
---|
568 | G4double zOne, G4double zTwo, |
---|
569 | G4double& totArea) const |
---|
570 | { |
---|
571 | // declare working variables |
---|
572 | // |
---|
573 | G4double Aone, Atwo, Afive, phi, zRand, fDPhi, fSPhi, cosu, sinu; |
---|
574 | G4double rRand1, chose, rone, rtwo, qone, qtwo, |
---|
575 | fDz = std::fabs((zTwo-zOne)/2.); |
---|
576 | G4ThreeVector point, offset; |
---|
577 | offset = G4ThreeVector(0.,0.,0.5*(zTwo+zOne)); |
---|
578 | fSPhi = startPhi; fDPhi = endPhi - startPhi; |
---|
579 | rone = (fRmax1-fRmax2)/(2.*fDz); |
---|
580 | rtwo = (fRmin1-fRmin2)/(2.*fDz); |
---|
581 | if(fRmax1==fRmax2){qone=0.;} |
---|
582 | else{ |
---|
583 | qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2); |
---|
584 | } |
---|
585 | if(fRmin1==fRmin2){qtwo=0.;} |
---|
586 | else{ |
---|
587 | qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2); |
---|
588 | } |
---|
589 | Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*(sqr(fRmin1-fRmin2)+sqr(zTwo-zOne)); |
---|
590 | Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*(sqr(fRmax1-fRmax2)+sqr(zTwo-zOne)); |
---|
591 | Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2); |
---|
592 | totArea = Aone+Atwo+2.*Afive; |
---|
593 | |
---|
594 | phi = RandFlat::shoot(startPhi,endPhi); |
---|
595 | cosu = std::cos(phi); |
---|
596 | sinu = std::sin(phi); |
---|
597 | |
---|
598 | |
---|
599 | if( (startPhi == 0) && (endPhi == twopi) ) { Afive = 0; } |
---|
600 | chose = RandFlat::shoot(0.,Aone+Atwo+2.*Afive); |
---|
601 | if( (chose >= 0) && (chose < Aone) ) |
---|
602 | { |
---|
603 | if(fRmax1 != fRmax2) |
---|
604 | { |
---|
605 | zRand = RandFlat::shoot(-1.*fDz,fDz); |
---|
606 | point = G4ThreeVector (rone*cosu*(qone-zRand), |
---|
607 | rone*sinu*(qone-zRand), zRand); |
---|
608 | |
---|
609 | |
---|
610 | } |
---|
611 | else |
---|
612 | { |
---|
613 | point = G4ThreeVector(fRmax1*cosu, fRmax1*sinu, |
---|
614 | RandFlat::shoot(-1.*fDz,fDz)); |
---|
615 | |
---|
616 | } |
---|
617 | } |
---|
618 | else if(chose >= Aone && chose < Aone + Atwo) |
---|
619 | { |
---|
620 | if(fRmin1 != fRmin2) |
---|
621 | { |
---|
622 | zRand = RandFlat::shoot(-1.*fDz,fDz); |
---|
623 | point = G4ThreeVector (rtwo*cosu*(qtwo-zRand), |
---|
624 | rtwo*sinu*(qtwo-zRand), zRand); |
---|
625 | |
---|
626 | } |
---|
627 | else |
---|
628 | { |
---|
629 | point = G4ThreeVector(fRmin1*cosu, fRmin1*sinu, |
---|
630 | RandFlat::shoot(-1.*fDz,fDz)); |
---|
631 | |
---|
632 | } |
---|
633 | } |
---|
634 | else if( (chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive) ) |
---|
635 | { |
---|
636 | zRand = RandFlat::shoot(-1.*fDz,fDz); |
---|
637 | rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2), |
---|
638 | fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2)); |
---|
639 | point = G4ThreeVector (rRand1*std::cos(startPhi), |
---|
640 | rRand1*std::sin(startPhi), zRand); |
---|
641 | } |
---|
642 | else |
---|
643 | { |
---|
644 | zRand = RandFlat::shoot(-1.*fDz,fDz); |
---|
645 | rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2), |
---|
646 | fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2)); |
---|
647 | point = G4ThreeVector (rRand1*std::cos(endPhi), |
---|
648 | rRand1*std::sin(endPhi), zRand); |
---|
649 | |
---|
650 | } |
---|
651 | return point+offset; |
---|
652 | } |
---|
653 | |
---|
654 | |
---|
655 | // |
---|
656 | // GetPointOnTubs |
---|
657 | // |
---|
658 | // Auxiliary method for GetPoint On Surface |
---|
659 | // |
---|
660 | G4ThreeVector G4Polycone::GetPointOnTubs(G4double fRMin, G4double fRMax, |
---|
661 | G4double zOne, G4double zTwo, |
---|
662 | G4double& totArea) const |
---|
663 | { |
---|
664 | G4double xRand,yRand,zRand,phi,cosphi,sinphi,chose, |
---|
665 | aOne,aTwo,aFou,rRand,fDz,fSPhi,fDPhi; |
---|
666 | fDz = std::fabs(0.5*(zTwo-zOne)); |
---|
667 | fSPhi = startPhi; |
---|
668 | fDPhi = endPhi-startPhi; |
---|
669 | |
---|
670 | aOne = 2.*fDz*fDPhi*fRMax; |
---|
671 | aTwo = 2.*fDz*fDPhi*fRMin; |
---|
672 | aFou = 2.*fDz*(fRMax-fRMin); |
---|
673 | totArea = aOne+aTwo+2.*aFou; |
---|
674 | phi = RandFlat::shoot(startPhi,endPhi); |
---|
675 | cosphi = std::cos(phi); |
---|
676 | sinphi = std::sin(phi); |
---|
677 | rRand = RandFlat::shoot(fRMin,fRMax); |
---|
678 | |
---|
679 | if(startPhi == 0 && endPhi == twopi) |
---|
680 | aFou = 0; |
---|
681 | |
---|
682 | chose = RandFlat::shoot(0.,aOne+aTwo+2.*aFou); |
---|
683 | if( (chose >= 0) && (chose < aOne) ) |
---|
684 | { |
---|
685 | xRand = fRMax*cosphi; |
---|
686 | yRand = fRMax*sinphi; |
---|
687 | zRand = RandFlat::shoot(-1.*fDz,fDz); |
---|
688 | return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne)); |
---|
689 | } |
---|
690 | else if( (chose >= aOne) && (chose < aOne + aTwo) ) |
---|
691 | { |
---|
692 | xRand = fRMin*cosphi; |
---|
693 | yRand = fRMin*sinphi; |
---|
694 | zRand = RandFlat::shoot(-1.*fDz,fDz); |
---|
695 | return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne)); |
---|
696 | } |
---|
697 | else if( (chose >= aOne+aTwo) && (chose <aOne+aTwo+aFou) ) |
---|
698 | { |
---|
699 | xRand = rRand*std::cos(fSPhi+fDPhi); |
---|
700 | yRand = rRand*std::sin(fSPhi+fDPhi); |
---|
701 | zRand = RandFlat::shoot(-1.*fDz,fDz); |
---|
702 | return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne)); |
---|
703 | } |
---|
704 | |
---|
705 | // else |
---|
706 | |
---|
707 | xRand = rRand*std::cos(fSPhi+fDPhi); |
---|
708 | yRand = rRand*std::sin(fSPhi+fDPhi); |
---|
709 | zRand = RandFlat::shoot(-1.*fDz,fDz); |
---|
710 | return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne)); |
---|
711 | } |
---|
712 | |
---|
713 | |
---|
714 | // |
---|
715 | // GetPointOnRing |
---|
716 | // |
---|
717 | // Auxiliary method for GetPoint On Surface |
---|
718 | // |
---|
719 | G4ThreeVector G4Polycone::GetPointOnRing(G4double fRMin1, G4double fRMax1, |
---|
720 | G4double fRMin2,G4double fRMax2, |
---|
721 | G4double zOne) const |
---|
722 | { |
---|
723 | G4double xRand,yRand,phi,cosphi,sinphi,rRand1,rRand2,A1,Atot,rCh; |
---|
724 | |
---|
725 | phi = RandFlat::shoot(startPhi,endPhi); |
---|
726 | cosphi = std::cos(phi); |
---|
727 | sinphi = std::sin(phi); |
---|
728 | |
---|
729 | if(fRMin1==fRMin2) |
---|
730 | { |
---|
731 | rRand1 = fRMin1; A1=0.; |
---|
732 | } |
---|
733 | else |
---|
734 | { |
---|
735 | rRand1 = RandFlat::shoot(fRMin1,fRMin2); |
---|
736 | A1=std::abs(fRMin2*fRMin2-fRMin1*fRMin1); |
---|
737 | } |
---|
738 | if(fRMax1==fRMax2) |
---|
739 | { |
---|
740 | rRand2=fRMax1; Atot=A1; |
---|
741 | } |
---|
742 | else |
---|
743 | { |
---|
744 | rRand2 = RandFlat::shoot(fRMax1,fRMax2); |
---|
745 | Atot = A1+std::abs(fRMax2*fRMax2-fRMax1*fRMax1); |
---|
746 | } |
---|
747 | rCh = RandFlat::shoot(0.,Atot); |
---|
748 | |
---|
749 | if(rCh>A1) { rRand1=rRand2; } |
---|
750 | |
---|
751 | xRand = rRand1*cosphi; |
---|
752 | yRand = rRand1*sinphi; |
---|
753 | |
---|
754 | return G4ThreeVector(xRand, yRand, zOne); |
---|
755 | } |
---|
756 | |
---|
757 | |
---|
758 | // |
---|
759 | // GetPointOnCut |
---|
760 | // |
---|
761 | // Auxiliary method for Get Point On Surface |
---|
762 | // |
---|
763 | G4ThreeVector G4Polycone::GetPointOnCut(G4double fRMin1, G4double fRMax1, |
---|
764 | G4double fRMin2, G4double fRMax2, |
---|
765 | G4double zOne, G4double zTwo, |
---|
766 | G4double& totArea) const |
---|
767 | { if(zOne==zTwo) |
---|
768 | { |
---|
769 | return GetPointOnRing(fRMin1, fRMax1,fRMin2,fRMax2,zOne); |
---|
770 | } |
---|
771 | if( (fRMin1 == fRMin2) && (fRMax1 == fRMax2) ) |
---|
772 | { |
---|
773 | return GetPointOnTubs(fRMin1, fRMax1,zOne,zTwo,totArea); |
---|
774 | } |
---|
775 | return GetPointOnCone(fRMin1,fRMax1,fRMin2,fRMax2,zOne,zTwo,totArea); |
---|
776 | } |
---|
777 | |
---|
778 | |
---|
779 | // |
---|
780 | // GetPointOnSurface |
---|
781 | // |
---|
782 | G4ThreeVector G4Polycone::GetPointOnSurface() const |
---|
783 | { |
---|
784 | if (!genericPcon) // Polycone by faces |
---|
785 | { |
---|
786 | G4double Area=0,totArea=0,Achose1=0,Achose2=0,phi,cosphi,sinphi,rRand; |
---|
787 | G4int i=0; |
---|
788 | G4int numPlanes = original_parameters->Num_z_planes; |
---|
789 | |
---|
790 | phi = RandFlat::shoot(startPhi,endPhi); |
---|
791 | cosphi = std::cos(phi); |
---|
792 | sinphi = std::sin(phi); |
---|
793 | |
---|
794 | rRand = RandFlat::shoot(original_parameters->Rmin[0], |
---|
795 | original_parameters->Rmax[0]); |
---|
796 | |
---|
797 | std::vector<G4double> areas; // (numPlanes+1); |
---|
798 | std::vector<G4ThreeVector> points; // (numPlanes-1); |
---|
799 | |
---|
800 | areas.push_back(pi*(sqr(original_parameters->Rmax[0]) |
---|
801 | -sqr(original_parameters->Rmin[0]))); |
---|
802 | |
---|
803 | for(i=0; i<numPlanes-1; i++) |
---|
804 | { |
---|
805 | Area = (original_parameters->Rmin[i]+original_parameters->Rmin[i+1]) |
---|
806 | * std::sqrt(sqr(original_parameters->Rmin[i] |
---|
807 | -original_parameters->Rmin[i+1])+ |
---|
808 | sqr(original_parameters->Z_values[i+1] |
---|
809 | -original_parameters->Z_values[i])); |
---|
810 | |
---|
811 | Area += (original_parameters->Rmax[i]+original_parameters->Rmax[i+1]) |
---|
812 | * std::sqrt(sqr(original_parameters->Rmax[i] |
---|
813 | -original_parameters->Rmax[i+1])+ |
---|
814 | sqr(original_parameters->Z_values[i+1] |
---|
815 | -original_parameters->Z_values[i])); |
---|
816 | |
---|
817 | Area *= 0.5*(endPhi-startPhi); |
---|
818 | |
---|
819 | if(startPhi==0.&& endPhi == twopi) |
---|
820 | { |
---|
821 | Area += std::fabs(original_parameters->Z_values[i+1] |
---|
822 | -original_parameters->Z_values[i])* |
---|
823 | (original_parameters->Rmax[i] |
---|
824 | +original_parameters->Rmax[i+1] |
---|
825 | -original_parameters->Rmin[i] |
---|
826 | -original_parameters->Rmin[i+1]); |
---|
827 | } |
---|
828 | areas.push_back(Area); |
---|
829 | totArea += Area; |
---|
830 | } |
---|
831 | |
---|
832 | areas.push_back(pi*(sqr(original_parameters->Rmax[numPlanes-1])- |
---|
833 | sqr(original_parameters->Rmin[numPlanes-1]))); |
---|
834 | |
---|
835 | totArea += (areas[0]+areas[numPlanes]); |
---|
836 | G4double chose = RandFlat::shoot(0.,totArea); |
---|
837 | |
---|
838 | if( (chose>=0.) && (chose<areas[0]) ) |
---|
839 | { |
---|
840 | return G4ThreeVector(rRand*cosphi, rRand*sinphi, |
---|
841 | original_parameters->Z_values[0]); |
---|
842 | } |
---|
843 | |
---|
844 | for (i=0; i<numPlanes-1; i++) |
---|
845 | { |
---|
846 | Achose1 += areas[i]; |
---|
847 | Achose2 = (Achose1+areas[i+1]); |
---|
848 | if(chose>=Achose1 && chose<Achose2) |
---|
849 | { |
---|
850 | return GetPointOnCut(original_parameters->Rmin[i], |
---|
851 | original_parameters->Rmax[i], |
---|
852 | original_parameters->Rmin[i+1], |
---|
853 | original_parameters->Rmax[i+1], |
---|
854 | original_parameters->Z_values[i], |
---|
855 | original_parameters->Z_values[i+1], Area); |
---|
856 | } |
---|
857 | } |
---|
858 | |
---|
859 | rRand = RandFlat::shoot(original_parameters->Rmin[numPlanes-1], |
---|
860 | original_parameters->Rmax[numPlanes-1]); |
---|
861 | |
---|
862 | return G4ThreeVector(rRand*cosphi,rRand*sinphi, |
---|
863 | original_parameters->Z_values[numPlanes-1]); |
---|
864 | |
---|
865 | } |
---|
866 | else // Generic Polycone |
---|
867 | { |
---|
868 | return GetPointOnSurfaceGeneric(); |
---|
869 | } |
---|
870 | } |
---|
871 | |
---|
872 | // |
---|
873 | // CreatePolyhedron |
---|
874 | // |
---|
875 | G4Polyhedron* G4Polycone::CreatePolyhedron() const |
---|
876 | { |
---|
877 | // |
---|
878 | // This has to be fixed in visualization. Fake it for the moment. |
---|
879 | // |
---|
880 | if (!genericPcon) |
---|
881 | { |
---|
882 | return new G4PolyhedronPcon( original_parameters->Start_angle, |
---|
883 | original_parameters->Opening_angle, |
---|
884 | original_parameters->Num_z_planes, |
---|
885 | original_parameters->Z_values, |
---|
886 | original_parameters->Rmin, |
---|
887 | original_parameters->Rmax ); |
---|
888 | } |
---|
889 | else |
---|
890 | { |
---|
891 | // The following code prepares for: |
---|
892 | // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces, |
---|
893 | // const double xyz[][3], |
---|
894 | // const int faces_vec[][4]) |
---|
895 | // Here is an extract from the header file HepPolyhedron.h: |
---|
896 | /** |
---|
897 | * Creates user defined polyhedron. |
---|
898 | * This function allows to the user to define arbitrary polyhedron. |
---|
899 | * The faces of the polyhedron should be either triangles or planar |
---|
900 | * quadrilateral. Nodes of a face are defined by indexes pointing to |
---|
901 | * the elements in the xyz array. Numeration of the elements in the |
---|
902 | * array starts from 1 (like in fortran). The indexes can be positive |
---|
903 | * or negative. Negative sign means that the corresponding edge is |
---|
904 | * invisible. The normal of the face should be directed to exterior |
---|
905 | * of the polyhedron. |
---|
906 | * |
---|
907 | * @param Nnodes number of nodes |
---|
908 | * @param Nfaces number of faces |
---|
909 | * @param xyz nodes |
---|
910 | * @param faces_vec faces (quadrilaterals or triangles) |
---|
911 | * @return status of the operation - is non-zero in case of problem |
---|
912 | */ |
---|
913 | const G4int numSide = |
---|
914 | G4int(G4Polyhedron::GetNumberOfRotationSteps() |
---|
915 | * (endPhi - startPhi) / twopi) + 1; |
---|
916 | G4int nNodes; |
---|
917 | G4int nFaces; |
---|
918 | typedef G4double double3[3]; |
---|
919 | double3* xyz; |
---|
920 | typedef G4int int4[4]; |
---|
921 | int4* faces_vec; |
---|
922 | if (phiIsOpen) |
---|
923 | { |
---|
924 | // Triangulate open ends. Simple ear-chopping algorithm... |
---|
925 | // I'm not sure how robust this algorithm is (J.Allison). |
---|
926 | // |
---|
927 | std::vector<G4bool> chopped(numCorner, false); |
---|
928 | std::vector<G4int*> triQuads; |
---|
929 | G4int remaining = numCorner; |
---|
930 | G4int iStarter = 0; |
---|
931 | while (remaining >= 3) |
---|
932 | { |
---|
933 | // Find unchopped corners... |
---|
934 | // |
---|
935 | G4int A = -1, B = -1, C = -1; |
---|
936 | G4int iStepper = iStarter; |
---|
937 | do |
---|
938 | { |
---|
939 | if (A < 0) { A = iStepper; } |
---|
940 | else if (B < 0) { B = iStepper; } |
---|
941 | else if (C < 0) { C = iStepper; } |
---|
942 | do |
---|
943 | { |
---|
944 | if (++iStepper >= numCorner) { iStepper = 0; } |
---|
945 | } |
---|
946 | while (chopped[iStepper]); |
---|
947 | } |
---|
948 | while (C < 0 && iStepper != iStarter); |
---|
949 | |
---|
950 | // Check triangle at B is pointing outward (an "ear"). |
---|
951 | // Sign of z cross product determines... |
---|
952 | // |
---|
953 | G4double BAr = corners[A].r - corners[B].r; |
---|
954 | G4double BAz = corners[A].z - corners[B].z; |
---|
955 | G4double BCr = corners[C].r - corners[B].r; |
---|
956 | G4double BCz = corners[C].z - corners[B].z; |
---|
957 | if (BAr * BCz - BAz * BCr < kCarTolerance) |
---|
958 | { |
---|
959 | G4int* tq = new G4int[3]; |
---|
960 | tq[0] = A + 1; |
---|
961 | tq[1] = B + 1; |
---|
962 | tq[2] = C + 1; |
---|
963 | triQuads.push_back(tq); |
---|
964 | chopped[B] = true; |
---|
965 | --remaining; |
---|
966 | } |
---|
967 | else |
---|
968 | { |
---|
969 | do |
---|
970 | { |
---|
971 | if (++iStarter >= numCorner) { iStarter = 0; } |
---|
972 | } |
---|
973 | while (chopped[iStarter]); |
---|
974 | } |
---|
975 | } |
---|
976 | // Transfer to faces... |
---|
977 | // |
---|
978 | nNodes = (numSide + 1) * numCorner; |
---|
979 | nFaces = numSide * numCorner + 2 * triQuads.size(); |
---|
980 | faces_vec = new int4[nFaces]; |
---|
981 | G4int iface = 0; |
---|
982 | G4int addition = numCorner * numSide; |
---|
983 | G4int d = numCorner - 1; |
---|
984 | for (G4int iEnd = 0; iEnd < 2; ++iEnd) |
---|
985 | { |
---|
986 | for (size_t i = 0; i < triQuads.size(); ++i) |
---|
987 | { |
---|
988 | // Negative for soft/auxiliary/normally invisible edges... |
---|
989 | // |
---|
990 | G4int a, b, c; |
---|
991 | if (iEnd == 0) |
---|
992 | { |
---|
993 | a = triQuads[i][0]; |
---|
994 | b = triQuads[i][1]; |
---|
995 | c = triQuads[i][2]; |
---|
996 | } |
---|
997 | else |
---|
998 | { |
---|
999 | a = triQuads[i][0] + addition; |
---|
1000 | b = triQuads[i][2] + addition; |
---|
1001 | c = triQuads[i][1] + addition; |
---|
1002 | } |
---|
1003 | G4int ab = std::abs(b - a); |
---|
1004 | G4int bc = std::abs(c - b); |
---|
1005 | G4int ca = std::abs(a - c); |
---|
1006 | faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a; |
---|
1007 | faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b; |
---|
1008 | faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c; |
---|
1009 | faces_vec[iface][3] = 0; |
---|
1010 | ++iface; |
---|
1011 | } |
---|
1012 | } |
---|
1013 | |
---|
1014 | // Continue with sides... |
---|
1015 | |
---|
1016 | xyz = new double3[nNodes]; |
---|
1017 | const G4double dPhi = (endPhi - startPhi) / numSide; |
---|
1018 | G4double phi = startPhi; |
---|
1019 | G4int ixyz = 0; |
---|
1020 | for (G4int iSide = 0; iSide < numSide; ++iSide) |
---|
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 | if (iSide == 0) // startPhi |
---|
1028 | { |
---|
1029 | if (iCorner < numCorner - 1) |
---|
1030 | { |
---|
1031 | faces_vec[iface][0] = ixyz + 1; |
---|
1032 | faces_vec[iface][1] = -(ixyz + numCorner + 1); |
---|
1033 | faces_vec[iface][2] = ixyz + numCorner + 2; |
---|
1034 | faces_vec[iface][3] = ixyz + 2; |
---|
1035 | } |
---|
1036 | else |
---|
1037 | { |
---|
1038 | faces_vec[iface][0] = ixyz + 1; |
---|
1039 | faces_vec[iface][1] = -(ixyz + numCorner + 1); |
---|
1040 | faces_vec[iface][2] = ixyz + 2; |
---|
1041 | faces_vec[iface][3] = ixyz - numCorner + 2; |
---|
1042 | } |
---|
1043 | } |
---|
1044 | else if (iSide == numSide - 1) // endPhi |
---|
1045 | { |
---|
1046 | if (iCorner < numCorner - 1) |
---|
1047 | { |
---|
1048 | faces_vec[iface][0] = ixyz + 1; |
---|
1049 | faces_vec[iface][1] = ixyz + numCorner + 1; |
---|
1050 | faces_vec[iface][2] = ixyz + numCorner + 2; |
---|
1051 | faces_vec[iface][3] = -(ixyz + 2); |
---|
1052 | } |
---|
1053 | else |
---|
1054 | { |
---|
1055 | faces_vec[iface][0] = ixyz + 1; |
---|
1056 | faces_vec[iface][1] = ixyz + numCorner + 1; |
---|
1057 | faces_vec[iface][2] = ixyz + 2; |
---|
1058 | faces_vec[iface][3] = -(ixyz - numCorner + 2); |
---|
1059 | } |
---|
1060 | } |
---|
1061 | else |
---|
1062 | { |
---|
1063 | if (iCorner < numCorner - 1) |
---|
1064 | { |
---|
1065 | faces_vec[iface][0] = ixyz + 1; |
---|
1066 | faces_vec[iface][1] = -(ixyz + numCorner + 1); |
---|
1067 | faces_vec[iface][2] = ixyz + numCorner + 2; |
---|
1068 | faces_vec[iface][3] = -(ixyz + 2); |
---|
1069 | } |
---|
1070 | else |
---|
1071 | { |
---|
1072 | faces_vec[iface][0] = ixyz + 1; |
---|
1073 | faces_vec[iface][1] = -(ixyz + numCorner + 1); |
---|
1074 | faces_vec[iface][2] = ixyz + 2; |
---|
1075 | faces_vec[iface][3] = -(ixyz - numCorner + 2); |
---|
1076 | } |
---|
1077 | } |
---|
1078 | ++iface; |
---|
1079 | ++ixyz; |
---|
1080 | } |
---|
1081 | phi += dPhi; |
---|
1082 | } |
---|
1083 | |
---|
1084 | // Last corners... |
---|
1085 | |
---|
1086 | for (G4int iCorner = 0; iCorner < numCorner; ++iCorner) |
---|
1087 | { |
---|
1088 | xyz[ixyz][0] = corners[iCorner].r * std::cos(phi); |
---|
1089 | xyz[ixyz][1] = corners[iCorner].r * std::sin(phi); |
---|
1090 | xyz[ixyz][2] = corners[iCorner].z; |
---|
1091 | ++ixyz; |
---|
1092 | } |
---|
1093 | } |
---|
1094 | else // !phiIsOpen - i.e., a complete 360 degrees. |
---|
1095 | { |
---|
1096 | nNodes = numSide * numCorner; |
---|
1097 | nFaces = numSide * numCorner;; |
---|
1098 | xyz = new double3[nNodes]; |
---|
1099 | faces_vec = new int4[nFaces]; |
---|
1100 | const G4double dPhi = (endPhi - startPhi) / numSide; |
---|
1101 | G4double phi = startPhi; |
---|
1102 | G4int ixyz = 0, iface = 0; |
---|
1103 | for (G4int iSide = 0; iSide < numSide; ++iSide) |
---|
1104 | { |
---|
1105 | for (G4int iCorner = 0; iCorner < numCorner; ++iCorner) |
---|
1106 | { |
---|
1107 | xyz[ixyz][0] = corners[iCorner].r * std::cos(phi); |
---|
1108 | xyz[ixyz][1] = corners[iCorner].r * std::sin(phi); |
---|
1109 | xyz[ixyz][2] = corners[iCorner].z; |
---|
1110 | |
---|
1111 | if (iSide < numSide - 1) |
---|
1112 | { |
---|
1113 | if (iCorner < numCorner - 1) |
---|
1114 | { |
---|
1115 | faces_vec[iface][0] = ixyz + 1; |
---|
1116 | faces_vec[iface][1] = -(ixyz + numCorner + 1); |
---|
1117 | faces_vec[iface][2] = ixyz + numCorner + 2; |
---|
1118 | faces_vec[iface][3] = -(ixyz + 2); |
---|
1119 | } |
---|
1120 | else |
---|
1121 | { |
---|
1122 | faces_vec[iface][0] = ixyz + 1; |
---|
1123 | faces_vec[iface][1] = -(ixyz + numCorner + 1); |
---|
1124 | faces_vec[iface][2] = ixyz + 2; |
---|
1125 | faces_vec[iface][3] = -(ixyz - numCorner + 2); |
---|
1126 | } |
---|
1127 | } |
---|
1128 | else // Last side joins ends... |
---|
1129 | { |
---|
1130 | if (iCorner < numCorner - 1) |
---|
1131 | { |
---|
1132 | faces_vec[iface][0] = ixyz + 1; |
---|
1133 | faces_vec[iface][1] = -(ixyz + numCorner - nFaces + 1); |
---|
1134 | faces_vec[iface][2] = ixyz + numCorner - nFaces + 2; |
---|
1135 | faces_vec[iface][3] = -(ixyz + 2); |
---|
1136 | } |
---|
1137 | else |
---|
1138 | { |
---|
1139 | faces_vec[iface][0] = ixyz + 1; |
---|
1140 | faces_vec[iface][1] = -(ixyz - nFaces + numCorner + 1); |
---|
1141 | faces_vec[iface][2] = ixyz - nFaces + 2; |
---|
1142 | faces_vec[iface][3] = -(ixyz - numCorner + 2); |
---|
1143 | } |
---|
1144 | } |
---|
1145 | ++ixyz; |
---|
1146 | ++iface; |
---|
1147 | } |
---|
1148 | phi += dPhi; |
---|
1149 | } |
---|
1150 | } |
---|
1151 | G4Polyhedron* polyhedron = new G4Polyhedron; |
---|
1152 | G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec); |
---|
1153 | delete [] faces_vec; |
---|
1154 | delete [] xyz; |
---|
1155 | if (problem) |
---|
1156 | { |
---|
1157 | std::ostringstream oss; |
---|
1158 | oss << "Problem creating G4Polyhedron for: " << GetName(); |
---|
1159 | G4Exception("G4Polycone::CreatePolyhedron()", "BadPolyhedron", |
---|
1160 | JustWarning, oss.str().c_str()); |
---|
1161 | delete polyhedron; |
---|
1162 | return 0; |
---|
1163 | } |
---|
1164 | else |
---|
1165 | { |
---|
1166 | return polyhedron; |
---|
1167 | } |
---|
1168 | } |
---|
1169 | } |
---|
1170 | |
---|
1171 | |
---|
1172 | // |
---|
1173 | // CreateNURBS |
---|
1174 | // |
---|
1175 | G4NURBS *G4Polycone::CreateNURBS() const |
---|
1176 | { |
---|
1177 | return 0; |
---|
1178 | } |
---|
1179 | |
---|
1180 | |
---|
1181 | // |
---|
1182 | // G4PolyconeHistorical stuff |
---|
1183 | // |
---|
1184 | |
---|
1185 | G4PolyconeHistorical::G4PolyconeHistorical() |
---|
1186 | : Z_values(0), Rmin(0), Rmax(0) |
---|
1187 | { |
---|
1188 | } |
---|
1189 | |
---|
1190 | G4PolyconeHistorical::~G4PolyconeHistorical() |
---|
1191 | { |
---|
1192 | delete [] Z_values; |
---|
1193 | delete [] Rmin; |
---|
1194 | delete [] Rmax; |
---|
1195 | } |
---|
1196 | |
---|
1197 | G4PolyconeHistorical:: |
---|
1198 | G4PolyconeHistorical( const G4PolyconeHistorical &source ) |
---|
1199 | { |
---|
1200 | Start_angle = source.Start_angle; |
---|
1201 | Opening_angle = source.Opening_angle; |
---|
1202 | Num_z_planes = source.Num_z_planes; |
---|
1203 | |
---|
1204 | Z_values = new G4double[Num_z_planes]; |
---|
1205 | Rmin = new G4double[Num_z_planes]; |
---|
1206 | Rmax = new G4double[Num_z_planes]; |
---|
1207 | |
---|
1208 | for( G4int i = 0; i < Num_z_planes; i++) |
---|
1209 | { |
---|
1210 | Z_values[i] = source.Z_values[i]; |
---|
1211 | Rmin[i] = source.Rmin[i]; |
---|
1212 | Rmax[i] = source.Rmax[i]; |
---|
1213 | } |
---|
1214 | } |
---|
1215 | |
---|
1216 | G4PolyconeHistorical& |
---|
1217 | G4PolyconeHistorical::operator=( const G4PolyconeHistorical& right ) |
---|
1218 | { |
---|
1219 | if ( &right == this ) return *this; |
---|
1220 | |
---|
1221 | if (&right) |
---|
1222 | { |
---|
1223 | Start_angle = right.Start_angle; |
---|
1224 | Opening_angle = right.Opening_angle; |
---|
1225 | Num_z_planes = right.Num_z_planes; |
---|
1226 | |
---|
1227 | delete [] Z_values; |
---|
1228 | delete [] Rmin; |
---|
1229 | delete [] Rmax; |
---|
1230 | Z_values = new G4double[Num_z_planes]; |
---|
1231 | Rmin = new G4double[Num_z_planes]; |
---|
1232 | Rmax = new G4double[Num_z_planes]; |
---|
1233 | |
---|
1234 | for( G4int i = 0; i < Num_z_planes; i++) |
---|
1235 | { |
---|
1236 | Z_values[i] = right.Z_values[i]; |
---|
1237 | Rmin[i] = right.Rmin[i]; |
---|
1238 | Rmax[i] = right.Rmax[i]; |
---|
1239 | } |
---|
1240 | } |
---|
1241 | return *this; |
---|
1242 | } |
---|