source: trunk/source/geometry/solids/specific/src/G4Polycone.cc@ 1347

Last change on this file since 1347 was 1337, checked in by garnier, 15 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 35.0 KB
RevLine 
[831]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//
[1337]27// $Id: G4Polycone.cc,v 1.44 2010/06/16 08:24:14 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
[831]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
53using namespace CLHEP;
54
55//
56// Constructor (GEANT3 style parameters)
57//
58G4Polycone::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//
122G4Polycone::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//
148void 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//
315G4Polycone::G4Polycone( __void__& a )
316 : G4VCSGfaceted(a), genericPcon(false), corners(0),
317 original_parameters(0), enclosingCylinder(0)
318{
319}
320
321
322//
323// Destructor
324//
325G4Polycone::~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//
337G4Polycone::G4Polycone( const G4Polycone &source )
338 : G4VCSGfaceted( source )
339{
340 CopyStuff( source );
341}
342
343
344//
345// Assignment operator
346//
347const 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//
367void 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//
409G4bool 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//
449EInside 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//
469G4double 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//
488G4double G4Polycone::DistanceToIn( const G4ThreeVector &p ) const
489{
490 return G4VCSGfaceted::DistanceToIn(p);
491}
492
493
494//
495// ComputeDimensions
496//
497void 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//
507G4GeometryType G4Polycone::GetEntityType() const
508{
509 return G4String("G4Polycone");
510}
511
512//
513// Stream object contents to an output stream
514//
515std::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//
566G4ThreeVector 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//
660G4ThreeVector 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//
719G4ThreeVector 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//
763G4ThreeVector 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//
782G4ThreeVector G4Polycone::GetPointOnSurface() const
[850]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;
[831]789
[850]790 phi = RandFlat::shoot(startPhi,endPhi);
791 cosphi = std::cos(phi);
792 sinphi = std::sin(phi);
[831]793
[850]794 rRand = RandFlat::shoot(original_parameters->Rmin[0],
795 original_parameters->Rmax[0]);
[831]796
[850]797 std::vector<G4double> areas; // (numPlanes+1);
798 std::vector<G4ThreeVector> points; // (numPlanes-1);
[831]799
[850]800 areas.push_back(pi*(sqr(original_parameters->Rmax[0])
801 -sqr(original_parameters->Rmin[0])));
[831]802
[850]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);
[831]818
[850]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);
[831]837
[850]838 if( (chose>=0.) && (chose<areas[0]) )
[831]839 {
[850]840 return G4ThreeVector(rRand*cosphi, rRand*sinphi,
841 original_parameters->Z_values[0]);
[831]842 }
843
[850]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]);
[831]861
[850]862 return G4ThreeVector(rRand*cosphi,rRand*sinphi,
863 original_parameters->Z_values[numPlanes-1]);
[831]864
865 }
[850]866 else // Generic Polycone
[831]867 {
[850]868 return GetPointOnSurfaceGeneric();
[831]869 }
870}
871
872//
873// CreatePolyhedron
874//
875G4Polyhedron* 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);
[1337]1153 delete [] faces_vec;
1154 delete [] xyz;
[831]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//
1175G4NURBS *G4Polycone::CreateNURBS() const
1176{
1177 return 0;
1178}
1179
1180
1181//
1182// G4PolyconeHistorical stuff
1183//
1184
1185G4PolyconeHistorical::G4PolyconeHistorical()
1186 : Z_values(0), Rmin(0), Rmax(0)
1187{
1188}
1189
1190G4PolyconeHistorical::~G4PolyconeHistorical()
1191{
1192 delete [] Z_values;
1193 delete [] Rmin;
1194 delete [] Rmax;
1195}
1196
1197G4PolyconeHistorical::
1198G4PolyconeHistorical( 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
1216G4PolyconeHistorical&
1217G4PolyconeHistorical::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}
Note: See TracBrowser for help on using the repository browser.