Geant4  10.00.p03
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  fRebuildPolyhedron = false;
330  fpPolyhedron = 0;
331 }
332 
333 
334 //
335 // Reset
336 //
338 {
339 
340  std::ostringstream message;
341  message << "Solid " << GetName() << " built using generic construct."
342  << G4endl << "Not applicable to the generic construct !";
343  G4Exception("G4GenericPolycone::Reset()", "GeomSolids1001",
344  JustWarning, message, "Parameters NOT resetted.");
345  return 1;
346 
347 }
348 
349 
350 //
351 // Inside
352 //
353 // This is an override of G4VCSGfaceted::Inside, created in order
354 // to speed things up by first checking with G4EnclosingCylinder.
355 //
357 {
358  //
359  // Quick test
360  //
361  if (enclosingCylinder->MustBeOutside(p)) return kOutside;
362 
363  //
364  // Long answer
365  //
366  return G4VCSGfaceted::Inside(p);
367 }
368 
369 
370 //
371 // DistanceToIn
372 //
373 // This is an override of G4VCSGfaceted::Inside, created in order
374 // to speed things up by first checking with G4EnclosingCylinder.
375 //
377  const G4ThreeVector &v ) const
378 {
379  //
380  // Quick test
381  //
382  if (enclosingCylinder->ShouldMiss(p,v))
383  return kInfinity;
384 
385  //
386  // Long answer
387  //
388  return G4VCSGfaceted::DistanceToIn( p, v );
389 }
390 
391 
392 //
393 // DistanceToIn
394 //
396 {
397  return G4VCSGfaceted::DistanceToIn(p);
398 }
399 
400 
401 //
402 // ComputeDimensions
403 //
404 /*void G4GenericPolycone::ComputeDimensions( G4VPVParameterisation* p,
405  const G4int n,
406  const G4VPhysicalVolume* pRep )
407 {
408  p->ComputeDimensions(*this,n,pRep);
409 }
410 */
411 //
412 // GetEntityType
413 //
415 {
416  return G4String("G4GenericPolycone");
417 }
418 
419 
420 //
421 // Make a clone of the object
422 //
424 {
425  return new G4GenericPolycone(*this);
426 }
427 
428 //
429 // Stream object contents to an output stream
430 //
431 std::ostream& G4GenericPolycone::StreamInfo( std::ostream& os ) const
432 {
433  G4int oldprc = os.precision(16);
434  os << "-----------------------------------------------------------\n"
435  << " *** Dump for solid - " << GetName() << " ***\n"
436  << " ===================================================\n"
437  << " Solid type: G4GenericPolycone\n"
438  << " Parameters: \n"
439  << " starting phi angle : " << startPhi/degree << " degrees \n"
440  << " ending phi angle : " << endPhi/degree << " degrees \n";
441  G4int i=0;
442 
443  os << " number of RZ points: " << numCorner << "\n"
444  << " RZ values (corners): \n";
445  for (i=0; i<numCorner; i++)
446  {
447  os << " "
448  << corners[i].r << ", " << corners[i].z << "\n";
449  }
450  os << "-----------------------------------------------------------\n";
451  os.precision(oldprc);
452 
453  return os;
454 }
455 
456 
457 
458 //
459 // GetPointOnSurface
460 //
462 {
463  return GetPointOnSurfaceGeneric();
464 
465 }
466 
467 //
468 // CreatePolyhedron
469 //
471 {
472  // The following code prepares for:
473  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
474  // const double xyz[][3],
475  // const int faces_vec[][4])
476  // Here is an extract from the header file HepPolyhedron.h:
494  const G4int numSide =
495  G4int(G4Polyhedron::GetNumberOfRotationSteps()
496  * (endPhi - startPhi) / twopi) + 1;
497  G4int nNodes;
498  G4int nFaces;
499  typedef G4double double3[3];
500  double3* xyz;
501  typedef G4int int4[4];
502  int4* faces_vec;
503  if (phiIsOpen)
504  {
505  // Triangulate open ends. Simple ear-chopping algorithm...
506  // I'm not sure how robust this algorithm is (J.Allison).
507  //
508  std::vector<G4bool> chopped(numCorner, false);
509  std::vector<G4int*> triQuads;
510  G4int remaining = numCorner;
511  G4int iStarter = 0;
512  while (remaining >= 3)
513  {
514  // Find unchopped corners...
515  //
516  G4int A = -1, B = -1, C = -1;
517  G4int iStepper = iStarter;
518  do
519  {
520  if (A < 0) { A = iStepper; }
521  else if (B < 0) { B = iStepper; }
522  else if (C < 0) { C = iStepper; }
523  do
524  {
525  if (++iStepper >= numCorner) { iStepper = 0; }
526  }
527  while (chopped[iStepper]);
528  }
529  while (C < 0 && iStepper != iStarter);
530 
531  // Check triangle at B is pointing outward (an "ear").
532  // Sign of z cross product determines...
533  //
534  G4double BAr = corners[A].r - corners[B].r;
535  G4double BAz = corners[A].z - corners[B].z;
536  G4double BCr = corners[C].r - corners[B].r;
537  G4double BCz = corners[C].z - corners[B].z;
538  if (BAr * BCz - BAz * BCr < kCarTolerance)
539  {
540  G4int* tq = new G4int[3];
541  tq[0] = A + 1;
542  tq[1] = B + 1;
543  tq[2] = C + 1;
544  triQuads.push_back(tq);
545  chopped[B] = true;
546  --remaining;
547  }
548  else
549  {
550  do
551  {
552  if (++iStarter >= numCorner) { iStarter = 0; }
553  }
554  while (chopped[iStarter]);
555  }
556  }
557  // Transfer to faces...
558  //
559  nNodes = (numSide + 1) * numCorner;
560  nFaces = numSide * numCorner + 2 * triQuads.size();
561  faces_vec = new int4[nFaces];
562  G4int iface = 0;
563  G4int addition = numCorner * numSide;
564  G4int d = numCorner - 1;
565  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
566  {
567  for (size_t i = 0; i < triQuads.size(); ++i)
568  {
569  // Negative for soft/auxiliary/normally invisible edges...
570  //
571  G4int a, b, c;
572  if (iEnd == 0)
573  {
574  a = triQuads[i][0];
575  b = triQuads[i][1];
576  c = triQuads[i][2];
577  }
578  else
579  {
580  a = triQuads[i][0] + addition;
581  b = triQuads[i][2] + addition;
582  c = triQuads[i][1] + addition;
583  }
584  G4int ab = std::abs(b - a);
585  G4int bc = std::abs(c - b);
586  G4int ca = std::abs(a - c);
587  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
588  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
589  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
590  faces_vec[iface][3] = 0;
591  ++iface;
592  }
593  }
594 
595  // Continue with sides...
596 
597  xyz = new double3[nNodes];
598  const G4double dPhi = (endPhi - startPhi) / numSide;
599  G4double phi = startPhi;
600  G4int ixyz = 0;
601  for (G4int iSide = 0; iSide < numSide; ++iSide)
602  {
603  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
604  {
605  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
606  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
607  xyz[ixyz][2] = corners[iCorner].z;
608  if (iSide == 0) // startPhi
609  {
610  if (iCorner < numCorner - 1)
611  {
612  faces_vec[iface][0] = ixyz + 1;
613  faces_vec[iface][1] = -(ixyz + numCorner + 1);
614  faces_vec[iface][2] = ixyz + numCorner + 2;
615  faces_vec[iface][3] = ixyz + 2;
616  }
617  else
618  {
619  faces_vec[iface][0] = ixyz + 1;
620  faces_vec[iface][1] = -(ixyz + numCorner + 1);
621  faces_vec[iface][2] = ixyz + 2;
622  faces_vec[iface][3] = ixyz - numCorner + 2;
623  }
624  }
625  else if (iSide == numSide - 1) // endPhi
626  {
627  if (iCorner < numCorner - 1)
628  {
629  faces_vec[iface][0] = ixyz + 1;
630  faces_vec[iface][1] = ixyz + numCorner + 1;
631  faces_vec[iface][2] = ixyz + numCorner + 2;
632  faces_vec[iface][3] = -(ixyz + 2);
633  }
634  else
635  {
636  faces_vec[iface][0] = ixyz + 1;
637  faces_vec[iface][1] = ixyz + numCorner + 1;
638  faces_vec[iface][2] = ixyz + 2;
639  faces_vec[iface][3] = -(ixyz - numCorner + 2);
640  }
641  }
642  else
643  {
644  if (iCorner < numCorner - 1)
645  {
646  faces_vec[iface][0] = ixyz + 1;
647  faces_vec[iface][1] = -(ixyz + numCorner + 1);
648  faces_vec[iface][2] = ixyz + numCorner + 2;
649  faces_vec[iface][3] = -(ixyz + 2);
650  }
651  else
652  {
653  faces_vec[iface][0] = ixyz + 1;
654  faces_vec[iface][1] = -(ixyz + numCorner + 1);
655  faces_vec[iface][2] = ixyz + 2;
656  faces_vec[iface][3] = -(ixyz - numCorner + 2);
657  }
658  }
659  ++iface;
660  ++ixyz;
661  }
662  phi += dPhi;
663  }
664 
665  // Last corners...
666 
667  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
668  {
669  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
670  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
671  xyz[ixyz][2] = corners[iCorner].z;
672  ++ixyz;
673  }
674  }
675  else // !phiIsOpen - i.e., a complete 360 degrees.
676  {
677  nNodes = numSide * numCorner;
678  nFaces = numSide * numCorner;;
679  xyz = new double3[nNodes];
680  faces_vec = new int4[nFaces];
681  const G4double dPhi = (endPhi - startPhi) / numSide;
682  G4double phi = startPhi;
683  G4int ixyz = 0, iface = 0;
684  for (G4int iSide = 0; iSide < numSide; ++iSide)
685  {
686  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
687  {
688  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
689  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
690  xyz[ixyz][2] = corners[iCorner].z;
691 
692  if (iSide < numSide - 1)
693  {
694  if (iCorner < numCorner - 1)
695  {
696  faces_vec[iface][0] = ixyz + 1;
697  faces_vec[iface][1] = -(ixyz + numCorner + 1);
698  faces_vec[iface][2] = ixyz + numCorner + 2;
699  faces_vec[iface][3] = -(ixyz + 2);
700  }
701  else
702  {
703  faces_vec[iface][0] = ixyz + 1;
704  faces_vec[iface][1] = -(ixyz + numCorner + 1);
705  faces_vec[iface][2] = ixyz + 2;
706  faces_vec[iface][3] = -(ixyz - numCorner + 2);
707  }
708  }
709  else // Last side joins ends...
710  {
711  if (iCorner < numCorner - 1)
712  {
713  faces_vec[iface][0] = ixyz + 1;
714  faces_vec[iface][1] = -(ixyz + numCorner - nFaces + 1);
715  faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
716  faces_vec[iface][3] = -(ixyz + 2);
717  }
718  else
719  {
720  faces_vec[iface][0] = ixyz + 1;
721  faces_vec[iface][1] = -(ixyz - nFaces + numCorner + 1);
722  faces_vec[iface][2] = ixyz - nFaces + 2;
723  faces_vec[iface][3] = -(ixyz - numCorner + 2);
724  }
725  }
726  ++ixyz;
727  ++iface;
728  }
729  phi += dPhi;
730  }
731  }
732  G4Polyhedron* polyhedron = new G4Polyhedron;
733  G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
734  delete [] faces_vec;
735  delete [] xyz;
736  if (problem)
737  {
738  std::ostringstream message;
739  message << "Problem creating G4Polyhedron for: " << GetName();
740  G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002",
741  JustWarning, message);
742  delete polyhedron;
743  return 0;
744  }
745  else
746  {
747  return polyhedron;
748  }
749 }
750 
751 #endif
G4String GetName() const
G4bool CrossesItself(G4double tolerance)
G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
G4bool fRebuildPolyhedron
G4Polyhedron * fpPolyhedron
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
static const G4double A[nN]
G4VCSGface ** faces
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4GenericPolycone & operator=(const G4GenericPolycone &source)
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