source: trunk/source/geometry/solids/specific/src/G4Polyhedra.cc@ 833

Last change on this file since 833 was 831, checked in by garnier, 17 years ago

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