Geant4  10.00.p01
G4GenericPolycone.cc
Go to the documentation of this file.
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: G4GenericPolycone.cc 72131 2013-07-11 07:12:24Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4GenericPolycone.cc
35 //
36 // Implementation of a CSG polycone
37 //
38 // --------------------------------------------------------------------
39 
40 #include "G4GenericPolycone.hh"
41 
42 #if !defined(G4GEOM_USE_UGENERICPOLYCONE)
43 
44 #include "G4PolyconeSide.hh"
45 #include "G4PolyPhiFace.hh"
46 
47 #include "Randomize.hh"
48 
49 #include "G4Polyhedron.hh"
50 #include "G4EnclosingCylinder.hh"
51 #include "G4ReduciblePolygon.hh"
52 #include "G4VPVParameterisation.hh"
53 
54 using namespace CLHEP;
55 
56 //
57 // Constructor (generic parameters)
58 //
60  G4double phiStart,
61  G4double phiTotal,
62  G4int numRZ,
63  const G4double r[],
64  const G4double z[] )
65  : G4VCSGfaceted( name )
66 {
67 
68  G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
69 
70  Create( phiStart, phiTotal, rz );
71 
72  // Set original_parameters struct for consistency
73  //
74  //SetOriginalParameters(rz);
75 
76  delete rz;
77 }
78 
79 //
80 // Create
81 //
82 // Generic create routine, called by each constructor after
83 // conversion of arguments
84 //
86  G4double phiTotal,
87  G4ReduciblePolygon *rz )
88 {
89  //
90  // Perform checks of rz values
91  //
92  if (rz->Amin() < 0.0)
93  {
94  std::ostringstream message;
95  message << "Illegal input parameters - " << GetName() << G4endl
96  << " All R values must be >= 0 !";
97  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
98  FatalErrorInArgument, message);
99  }
100 
101  G4double rzArea = rz->Area();
102  if (rzArea < -kCarTolerance)
103  rz->ReverseOrder();
104 
105  else if (rzArea < -kCarTolerance)
106  {
107  std::ostringstream message;
108  message << "Illegal input parameters - " << GetName() << G4endl
109  << " R/Z cross section is zero or near zero: " << rzArea;
110  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
111  FatalErrorInArgument, message);
112  }
113 
115  || (!rz->RemoveRedundantVertices( kCarTolerance )) )
116  {
117  std::ostringstream message;
118  message << "Illegal input parameters - " << GetName() << G4endl
119  << " Too few unique R/Z values !";
120  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
121  FatalErrorInArgument, message);
122  }
123 
124  if (rz->CrossesItself(1/kInfinity))
125  {
126  std::ostringstream message;
127  message << "Illegal input parameters - " << GetName() << G4endl
128  << " R/Z segments cross !";
129  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
130  FatalErrorInArgument, message);
131  }
132 
133  numCorner = rz->NumVertices();
134 
135  //
136  // Phi opening? Account for some possible roundoff, and interpret
137  // nonsense value as representing no phi opening
138  //
139  if (phiTotal <= 0 || phiTotal > twopi-1E-10)
140  {
141  phiIsOpen = false;
142  startPhi = 0;
143  endPhi = twopi;
144  }
145  else
146  {
147  phiIsOpen = true;
148 
149  //
150  // Convert phi into our convention
151  //
152  startPhi = phiStart;
153  while( startPhi < 0 ) startPhi += twopi;
154 
155  endPhi = phiStart+phiTotal;
156  while( endPhi < startPhi ) endPhi += twopi;
157  }
158 
159  //
160  // Allocate corner array.
161  //
163 
164  //
165  // Copy corners
166  //
167  G4ReduciblePolygonIterator iterRZ(rz);
168 
169  G4PolyconeSideRZ *next = corners;
170  iterRZ.Begin();
171  do
172  {
173  next->r = iterRZ.GetA();
174  next->z = iterRZ.GetB();
175  } while( ++next, iterRZ.Next() );
176 
177  //
178  // Allocate face pointer array
179  //
181  faces = new G4VCSGface*[numFace];
182 
183  //
184  // Construct conical faces
185  //
186  // But! Don't construct a face if both points are at zero radius!
187  //
188  G4PolyconeSideRZ *corner = corners,
189  *prev = corners + numCorner-1,
190  *nextNext;
191  G4VCSGface **face = faces;
192  do
193  {
194  next = corner+1;
195  if (next >= corners+numCorner) next = corners;
196  nextNext = next+1;
197  if (nextNext >= corners+numCorner) nextNext = corners;
198 
199  if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
200 
201  //
202  // We must decide here if we can dare declare one of our faces
203  // as having a "valid" normal (i.e. allBehind = true). This
204  // is never possible if the face faces "inward" in r.
205  //
206  G4bool allBehind;
207  if (corner->z > next->z)
208  {
209  allBehind = false;
210  }
211  else
212  {
213  //
214  // Otherwise, it is only true if the line passing
215  // through the two points of the segment do not
216  // split the r/z cross section
217  //
218  allBehind = !rz->BisectedBy( corner->r, corner->z,
219  next->r, next->z, kCarTolerance );
220  }
221 
222  *face++ = new G4PolyconeSide( prev, corner, next, nextNext,
223  startPhi, endPhi-startPhi, phiIsOpen, allBehind );
224  } while( prev=corner, corner=next, corner > corners );
225 
226  if (phiIsOpen)
227  {
228  //
229  // Construct phi open edges
230  //
231  *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi );
232  *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi );
233  }
234 
235  //
236  // We might have dropped a face or two: recalculate numFace
237  //
238  numFace = face-faces;
239 
240  //
241  // Make enclosingCylinder
242  //
244  new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
245 }
246 
247 
248 //
249 // Fake default constructor - sets only member data and allocates memory
250 // for usage restricted to object persistency.
251 //
253  : G4VCSGfaceted(a), startPhi(0.), endPhi(0.), phiIsOpen(false),
254  numCorner(0), corners(0), enclosingCylinder(0)
255 {
256 }
257 
258 
259 //
260 // Destructor
261 //
263 {
264  delete [] corners;
265  delete enclosingCylinder;
266 }
267 
268 
269 //
270 // Copy constructor
271 //
273  : G4VCSGfaceted( source )
274 {
275  CopyStuff( source );
276 }
277 
278 
279 //
280 // Assignment operator
281 //
283 {
284  if (this == &source) return *this;
285 
286  G4VCSGfaceted::operator=( source );
287 
288  delete [] corners;
289  // if (original_parameters) delete original_parameters;
290 
291  delete enclosingCylinder;
292 
293  CopyStuff( source );
294 
295  return *this;
296 }
297 
298 
299 //
300 // CopyStuff
301 //
303 {
304  //
305  // Simple stuff
306  //
307  startPhi = source.startPhi;
308  endPhi = source.endPhi;
309  phiIsOpen = source.phiIsOpen;
310  numCorner = source.numCorner;
311 
312  //
313  // The corner array
314  //
316 
317  G4PolyconeSideRZ *corn = corners,
318  *sourceCorn = source.corners;
319  do
320  {
321  *corn = *sourceCorn;
322  } while( ++sourceCorn, ++corn < corners+numCorner );
323 
324  //
325  // Enclosing cylinder
326  //
328 }
329 
330 
331 //
332 // Reset
333 //
335 {
336 
337  std::ostringstream message;
338  message << "Solid " << GetName() << " built using generic construct."
339  << G4endl << "Not applicable to the generic construct !";
340  G4Exception("G4GenericPolycone::Reset()", "GeomSolids1001",
341  JustWarning, message, "Parameters NOT resetted.");
342  return 1;
343 
344 }
345 
346 
347 //
348 // Inside
349 //
350 // This is an override of G4VCSGfaceted::Inside, created in order
351 // to speed things up by first checking with G4EnclosingCylinder.
352 //
354 {
355  //
356  // Quick test
357  //
358  if (enclosingCylinder->MustBeOutside(p)) return kOutside;
359 
360  //
361  // Long answer
362  //
363  return G4VCSGfaceted::Inside(p);
364 }
365 
366 
367 //
368 // DistanceToIn
369 //
370 // This is an override of G4VCSGfaceted::Inside, created in order
371 // to speed things up by first checking with G4EnclosingCylinder.
372 //
374  const G4ThreeVector &v ) const
375 {
376  //
377  // Quick test
378  //
379  if (enclosingCylinder->ShouldMiss(p,v))
380  return kInfinity;
381 
382  //
383  // Long answer
384  //
385  return G4VCSGfaceted::DistanceToIn( p, v );
386 }
387 
388 
389 //
390 // DistanceToIn
391 //
393 {
394  return G4VCSGfaceted::DistanceToIn(p);
395 }
396 
397 
398 //
399 // ComputeDimensions
400 //
401 /*void G4GenericPolycone::ComputeDimensions( G4VPVParameterisation* p,
402  const G4int n,
403  const G4VPhysicalVolume* pRep )
404 {
405  p->ComputeDimensions(*this,n,pRep);
406 }
407 */
408 //
409 // GetEntityType
410 //
412 {
413  return G4String("G4GenericPolycone");
414 }
415 
416 
417 //
418 // Make a clone of the object
419 //
421 {
422  return new G4GenericPolycone(*this);
423 }
424 
425 //
426 // Stream object contents to an output stream
427 //
428 std::ostream& G4GenericPolycone::StreamInfo( std::ostream& os ) const
429 {
430  G4int oldprc = os.precision(16);
431  os << "-----------------------------------------------------------\n"
432  << " *** Dump for solid - " << GetName() << " ***\n"
433  << " ===================================================\n"
434  << " Solid type: G4GenericPolycone\n"
435  << " Parameters: \n"
436  << " starting phi angle : " << startPhi/degree << " degrees \n"
437  << " ending phi angle : " << endPhi/degree << " degrees \n";
438  G4int i=0;
439 
440  os << " number of RZ points: " << numCorner << "\n"
441  << " RZ values (corners): \n";
442  for (i=0; i<numCorner; i++)
443  {
444  os << " "
445  << corners[i].r << ", " << corners[i].z << "\n";
446  }
447  os << "-----------------------------------------------------------\n";
448  os.precision(oldprc);
449 
450  return os;
451 }
452 
453 
454 
455 //
456 // GetPointOnSurface
457 //
459 {
460  return GetPointOnSurfaceGeneric();
461 
462 }
463 
464 //
465 // CreatePolyhedron
466 //
468 {
469  // The following code prepares for:
470  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
471  // const double xyz[][3],
472  // const int faces_vec[][4])
473  // Here is an extract from the header file HepPolyhedron.h:
491  const G4int numSide =
492  G4int(G4Polyhedron::GetNumberOfRotationSteps()
493  * (endPhi - startPhi) / twopi) + 1;
494  G4int nNodes;
495  G4int nFaces;
496  typedef G4double double3[3];
497  double3* xyz;
498  typedef G4int int4[4];
499  int4* faces_vec;
500  if (phiIsOpen)
501  {
502  // Triangulate open ends. Simple ear-chopping algorithm...
503  // I'm not sure how robust this algorithm is (J.Allison).
504  //
505  std::vector<G4bool> chopped(numCorner, false);
506  std::vector<G4int*> triQuads;
507  G4int remaining = numCorner;
508  G4int iStarter = 0;
509  while (remaining >= 3)
510  {
511  // Find unchopped corners...
512  //
513  G4int A = -1, B = -1, C = -1;
514  G4int iStepper = iStarter;
515  do
516  {
517  if (A < 0) { A = iStepper; }
518  else if (B < 0) { B = iStepper; }
519  else if (C < 0) { C = iStepper; }
520  do
521  {
522  if (++iStepper >= numCorner) { iStepper = 0; }
523  }
524  while (chopped[iStepper]);
525  }
526  while (C < 0 && iStepper != iStarter);
527 
528  // Check triangle at B is pointing outward (an "ear").
529  // Sign of z cross product determines...
530  //
531  G4double BAr = corners[A].r - corners[B].r;
532  G4double BAz = corners[A].z - corners[B].z;
533  G4double BCr = corners[C].r - corners[B].r;
534  G4double BCz = corners[C].z - corners[B].z;
535  if (BAr * BCz - BAz * BCr < kCarTolerance)
536  {
537  G4int* tq = new G4int[3];
538  tq[0] = A + 1;
539  tq[1] = B + 1;
540  tq[2] = C + 1;
541  triQuads.push_back(tq);
542  chopped[B] = true;
543  --remaining;
544  }
545  else
546  {
547  do
548  {
549  if (++iStarter >= numCorner) { iStarter = 0; }
550  }
551  while (chopped[iStarter]);
552  }
553  }
554  // Transfer to faces...
555  //
556  nNodes = (numSide + 1) * numCorner;
557  nFaces = numSide * numCorner + 2 * triQuads.size();
558  faces_vec = new int4[nFaces];
559  G4int iface = 0;
560  G4int addition = numCorner * numSide;
561  G4int d = numCorner - 1;
562  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
563  {
564  for (size_t i = 0; i < triQuads.size(); ++i)
565  {
566  // Negative for soft/auxiliary/normally invisible edges...
567  //
568  G4int a, b, c;
569  if (iEnd == 0)
570  {
571  a = triQuads[i][0];
572  b = triQuads[i][1];
573  c = triQuads[i][2];
574  }
575  else
576  {
577  a = triQuads[i][0] + addition;
578  b = triQuads[i][2] + addition;
579  c = triQuads[i][1] + addition;
580  }
581  G4int ab = std::abs(b - a);
582  G4int bc = std::abs(c - b);
583  G4int ca = std::abs(a - c);
584  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
585  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
586  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
587  faces_vec[iface][3] = 0;
588  ++iface;
589  }
590  }
591 
592  // Continue with sides...
593 
594  xyz = new double3[nNodes];
595  const G4double dPhi = (endPhi - startPhi) / numSide;
596  G4double phi = startPhi;
597  G4int ixyz = 0;
598  for (G4int iSide = 0; iSide < numSide; ++iSide)
599  {
600  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
601  {
602  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
603  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
604  xyz[ixyz][2] = corners[iCorner].z;
605  if (iSide == 0) // startPhi
606  {
607  if (iCorner < numCorner - 1)
608  {
609  faces_vec[iface][0] = ixyz + 1;
610  faces_vec[iface][1] = -(ixyz + numCorner + 1);
611  faces_vec[iface][2] = ixyz + numCorner + 2;
612  faces_vec[iface][3] = ixyz + 2;
613  }
614  else
615  {
616  faces_vec[iface][0] = ixyz + 1;
617  faces_vec[iface][1] = -(ixyz + numCorner + 1);
618  faces_vec[iface][2] = ixyz + 2;
619  faces_vec[iface][3] = ixyz - numCorner + 2;
620  }
621  }
622  else if (iSide == numSide - 1) // endPhi
623  {
624  if (iCorner < numCorner - 1)
625  {
626  faces_vec[iface][0] = ixyz + 1;
627  faces_vec[iface][1] = ixyz + numCorner + 1;
628  faces_vec[iface][2] = ixyz + numCorner + 2;
629  faces_vec[iface][3] = -(ixyz + 2);
630  }
631  else
632  {
633  faces_vec[iface][0] = ixyz + 1;
634  faces_vec[iface][1] = ixyz + numCorner + 1;
635  faces_vec[iface][2] = ixyz + 2;
636  faces_vec[iface][3] = -(ixyz - numCorner + 2);
637  }
638  }
639  else
640  {
641  if (iCorner < numCorner - 1)
642  {
643  faces_vec[iface][0] = ixyz + 1;
644  faces_vec[iface][1] = -(ixyz + numCorner + 1);
645  faces_vec[iface][2] = ixyz + numCorner + 2;
646  faces_vec[iface][3] = -(ixyz + 2);
647  }
648  else
649  {
650  faces_vec[iface][0] = ixyz + 1;
651  faces_vec[iface][1] = -(ixyz + numCorner + 1);
652  faces_vec[iface][2] = ixyz + 2;
653  faces_vec[iface][3] = -(ixyz - numCorner + 2);
654  }
655  }
656  ++iface;
657  ++ixyz;
658  }
659  phi += dPhi;
660  }
661 
662  // Last corners...
663 
664  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
665  {
666  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
667  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
668  xyz[ixyz][2] = corners[iCorner].z;
669  ++ixyz;
670  }
671  }
672  else // !phiIsOpen - i.e., a complete 360 degrees.
673  {
674  nNodes = numSide * numCorner;
675  nFaces = numSide * numCorner;;
676  xyz = new double3[nNodes];
677  faces_vec = new int4[nFaces];
678  const G4double dPhi = (endPhi - startPhi) / numSide;
679  G4double phi = startPhi;
680  G4int ixyz = 0, iface = 0;
681  for (G4int iSide = 0; iSide < numSide; ++iSide)
682  {
683  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
684  {
685  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
686  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
687  xyz[ixyz][2] = corners[iCorner].z;
688 
689  if (iSide < numSide - 1)
690  {
691  if (iCorner < numCorner - 1)
692  {
693  faces_vec[iface][0] = ixyz + 1;
694  faces_vec[iface][1] = -(ixyz + numCorner + 1);
695  faces_vec[iface][2] = ixyz + numCorner + 2;
696  faces_vec[iface][3] = -(ixyz + 2);
697  }
698  else
699  {
700  faces_vec[iface][0] = ixyz + 1;
701  faces_vec[iface][1] = -(ixyz + numCorner + 1);
702  faces_vec[iface][2] = ixyz + 2;
703  faces_vec[iface][3] = -(ixyz - numCorner + 2);
704  }
705  }
706  else // Last side joins ends...
707  {
708  if (iCorner < numCorner - 1)
709  {
710  faces_vec[iface][0] = ixyz + 1;
711  faces_vec[iface][1] = -(ixyz + numCorner - nFaces + 1);
712  faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
713  faces_vec[iface][3] = -(ixyz + 2);
714  }
715  else
716  {
717  faces_vec[iface][0] = ixyz + 1;
718  faces_vec[iface][1] = -(ixyz - nFaces + numCorner + 1);
719  faces_vec[iface][2] = ixyz - nFaces + 2;
720  faces_vec[iface][3] = -(ixyz - numCorner + 2);
721  }
722  }
723  ++ixyz;
724  ++iface;
725  }
726  phi += dPhi;
727  }
728  }
729  G4Polyhedron* polyhedron = new G4Polyhedron;
730  G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
731  delete [] faces_vec;
732  delete [] xyz;
733  if (problem)
734  {
735  std::ostringstream message;
736  message << "Problem creating G4Polyhedron for: " << GetName();
737  G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002",
738  JustWarning, message);
739  delete polyhedron;
740  return 0;
741  }
742  else
743  {
744  return polyhedron;
745  }
746 }
747 
748 #endif
G4String GetName() const
G4bool CrossesItself(G4double tolerance)
const G4GenericPolycone & operator=(const G4GenericPolycone &source)
static const G4double kInfinity
Definition: geomdefs.hh:42
G4PolyconeSideRZ * corners
CLHEP::Hep3Vector G4ThreeVector
G4double Amin() const
G4double z
Definition: TRTMaterials.hh:39
G4String name
Definition: TRTMaterials.hh:40
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double a
Definition: TRTMaterials.hh:39
G4bool MustBeOutside(const G4ThreeVector &p) const
int G4int
Definition: G4Types.hh:78
G4GenericPolycone(const G4String &name, G4double phiStart, G4double phiTotal, G4int numRZ, const G4double r[], const G4double z[])
G4GeometryType GetEntityType() const
G4Polyhedron * CreatePolyhedron() const
G4bool RemoveDuplicateVertices(G4double tolerance)
G4ThreeVector GetPointOnSurface() const
G4EnclosingCylinder * enclosingCylinder
void Create(G4double phiStart, G4double phiTotal, G4ReduciblePolygon *rz)
G4bool RemoveRedundantVertices(G4double tolerance)
G4ThreeVector GetPointOnSurfaceGeneric() const
G4int NumVertices() const
void CopyStuff(const G4GenericPolycone &source)
bool G4bool
Definition: G4Types.hh:79
std::ostream & StreamInfo(std::ostream &os) const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
const G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
static const G4double A[nN]
G4VCSGface ** faces
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4double ab
EInside Inside(const G4ThreeVector &p) const
EInside
Definition: geomdefs.hh:58
static const double degree
Definition: G4SIunits.hh:125
G4bool BisectedBy(G4double a1, G4double b1, G4double a2, G4double b2, G4double tolerance)
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:305
double G4double
Definition: G4Types.hh:76
G4VSolid * Clone() const
virtual EInside Inside(const G4ThreeVector &p) const
G4bool ShouldMiss(const G4ThreeVector &p, const G4ThreeVector &v) const