Geant4_10
G4VTwistedFaceted.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 // $Id: G4VTwistedFaceted.cc 66356 2012-12-18 09:02:32Z gcosmo $
27 //
28 //
29 // --------------------------------------------------------------------
30 // GEANT 4 class source file
31 //
32 //
33 // G4VTwistedFaceted.cc
34 //
35 // Author:
36 //
37 // 04-Nov-2004 - O.Link (Oliver.Link@cern.ch)
38 //
39 // --------------------------------------------------------------------
40 
41 #include "G4VTwistedFaceted.hh"
42 
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4VoxelLimits.hh"
46 #include "G4AffineTransform.hh"
47 #include "G4SolidExtentList.hh"
48 #include "G4ClippablePolygon.hh"
49 #include "G4VPVParameterisation.hh"
50 #include "G4GeometryTolerance.hh"
51 #include "meshdefs.hh"
52 
53 #include "G4VGraphicsScene.hh"
54 #include "G4Polyhedron.hh"
55 #include "G4VisExtent.hh"
56 
57 #include "Randomize.hh"
58 
59 //=====================================================================
60 //* constructors ------------------------------------------------------
61 
63 G4VTwistedFaceted( const G4String &pname, // Name of instance
64  G4double PhiTwist, // twist angle
65  G4double pDz, // half z length
66  G4double pTheta, // direction between end planes
67  G4double pPhi, // defined by polar and azim. angles
68  G4double pDy1, // half y length at -pDz
69  G4double pDx1, // half x length at -pDz,-pDy
70  G4double pDx2, // half x length at -pDz,+pDy
71  G4double pDy2, // half y length at +pDz
72  G4double pDx3, // half x length at +pDz,-pDy
73  G4double pDx4, // half x length at +pDz,+pDy
74  G4double pAlph // tilt angle
75  )
76  : G4VSolid(pname),
77  fLowerEndcap(0), fUpperEndcap(0), fSide0(0),
78  fSide90(0), fSide180(0), fSide270(0),
79  fSurfaceArea(0.), fpPolyhedron(0)
80 {
81 
82  G4double pDytmp ;
83  G4double fDxUp ;
84  G4double fDxDown ;
85 
86  fDx1 = pDx1 ;
87  fDx2 = pDx2 ;
88  fDx3 = pDx3 ;
89  fDx4 = pDx4 ;
90  fDy1 = pDy1 ;
91  fDy2 = pDy2 ;
92  fDz = pDz ;
93 
94  G4double kAngTolerance
96 
97  // maximum values
98  //
99  fDxDown = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ;
100  fDxUp = ( fDx3 > fDx4 ? fDx3 : fDx4 ) ;
101  fDx = ( fDxUp > fDxDown ? fDxUp : fDxDown ) ;
102  fDy = ( fDy1 > fDy2 ? fDy1 : fDy2 ) ;
103 
104  // planarity check
105  //
106  if ( fDx1 != fDx2 && fDx3 != fDx4 )
107  {
108  pDytmp = fDy1 * ( fDx3 - fDx4 ) / ( fDx1 - fDx2 ) ;
109  if ( std::fabs(pDytmp - fDy2) > kCarTolerance )
110  {
111  std::ostringstream message;
112  message << "Not planar surface in untwisted Trapezoid: "
113  << GetName() << G4endl
114  << "fDy2 is " << fDy2 << " but should be "
115  << pDytmp << ".";
116  G4Exception("G4VTwistedFaceted::G4VTwistedFaceted()", "GeomSolids0002",
117  FatalErrorInArgument, message);
118  }
119  }
120 
121 #ifdef G4TWISTDEBUG
122  if ( fDx1 == fDx2 && fDx3 == fDx4 )
123  {
124  G4cout << "Trapezoid is a box" << G4endl ;
125  }
126 
127 #endif
128 
129  if ( ( fDx1 == fDx2 && fDx3 != fDx4 ) || ( fDx1 != fDx2 && fDx3 == fDx4 ) )
130  {
131  std::ostringstream message;
132  message << "Not planar surface in untwisted Trapezoid: "
133  << GetName() << G4endl
134  << "One endcap is rectangular, the other is a trapezoid." << G4endl
135  << "For planarity reasons they have to be rectangles or trapezoids "
136  << "on both sides.";
137  G4Exception("G4VTwistedFaceted::G4VTwistedFaceted()", "GeomSolids0002",
138  FatalErrorInArgument, message);
139  }
140 
141  // twist angle
142  //
143  fPhiTwist = PhiTwist ;
144 
145  // tilt angle
146  //
147  fAlph = pAlph ;
148  fTAlph = std::tan(fAlph) ;
149 
150  fTheta = pTheta ;
151  fPhi = pPhi ;
152 
153  // dx in surface equation
154  //
155  fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
156 
157  // dy in surface equation
158  //
159  fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
160 
161  if ( ! ( ( fDx1 > 2*kCarTolerance)
162  && ( fDx2 > 2*kCarTolerance)
163  && ( fDx3 > 2*kCarTolerance)
164  && ( fDx4 > 2*kCarTolerance)
165  && ( fDy1 > 2*kCarTolerance)
166  && ( fDy2 > 2*kCarTolerance)
167  && ( fDz > 2*kCarTolerance)
168  && ( std::fabs(fPhiTwist) > 2*kAngTolerance )
169  && ( std::fabs(fPhiTwist) < pi/2 )
170  && ( std::fabs(fAlph) < pi/2 )
171  && ( fTheta < pi/2 && fTheta >= 0 ) )
172  )
173  {
174  std::ostringstream message;
175  message << "Invalid dimensions. Too small, or twist angle too big: "
176  << GetName() << G4endl
177  << "fDx 1-4 = " << fDx1/cm << ", " << fDx2/cm << ", "
178  << fDx3/cm << ", " << fDx4/cm << " cm" << G4endl
179  << "fDy 1-2 = " << fDy1/cm << ", " << fDy2/cm << ", "
180  << " cm" << G4endl
181  << "fDz = " << fDz/cm << " cm" << G4endl
182  << " twistangle " << fPhiTwist/deg << " deg" << G4endl
183  << " phi,theta = " << fPhi/deg << ", " << fTheta/deg << " deg";
184  G4Exception("G4TwistedTrap::G4VTwistedFaceted()",
185  "GeomSolids0002", FatalErrorInArgument, message);
186  }
187  CreateSurfaces();
188  fCubicVolume = 2 * fDz * ( ( fDx1 + fDx2 ) * fDy1 + ( fDx3 + fDx4 ) * fDy2 );
189 }
190 
191 
192 //=====================================================================
193 //* Fake default constructor ------------------------------------------
194 
196  : G4VSolid(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
197  fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fDx(0.), fDy(0.), fAlph(0.),
198  fTAlph(0.), fdeltaX(0.), fdeltaY(0.), fPhiTwist(0.),
199  fLowerEndcap(0), fUpperEndcap(0), fSide0(0), fSide90(0), fSide180(0),
200  fSide270(0), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
201 {
202 }
203 
204 //=====================================================================
205 //* destructor --------------------------------------------------------
206 
208 {
209  if (fLowerEndcap) delete fLowerEndcap ;
210  if (fUpperEndcap) delete fUpperEndcap ;
211 
212  if (fSide0) delete fSide0 ;
213  if (fSide90) delete fSide90 ;
214  if (fSide180) delete fSide180 ;
215  if (fSide270) delete fSide270 ;
216  if (fpPolyhedron) delete fpPolyhedron;
217 }
218 
219 //=====================================================================
220 //* Copy constructor --------------------------------------------------
221 
223  : G4VSolid(rhs), fTheta(rhs.fTheta), fPhi(rhs.fPhi),
224  fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fDy2(rhs.fDy2),
225  fDx3(rhs.fDx3), fDx4(rhs.fDx4), fDz(rhs.fDz), fDx(rhs.fDx), fDy(rhs.fDy),
226  fAlph(rhs.fAlph), fTAlph(rhs.fTAlph), fdeltaX(rhs.fdeltaX),
227  fdeltaY(rhs.fdeltaY), fPhiTwist(rhs.fPhiTwist), fLowerEndcap(0),
228  fUpperEndcap(0), fSide0(0), fSide90(0), fSide180(0), fSide270(0),
229  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
230  fpPolyhedron(0),
231  fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
232  fLastDistanceToIn(rhs.fLastDistanceToIn),
233  fLastDistanceToOut(rhs.fLastDistanceToOut),
234  fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
235  fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
236 {
237  CreateSurfaces();
238 }
239 
240 
241 //=====================================================================
242 //* Assignment operator -----------------------------------------------
243 
245 {
246  // Check assignment to self
247  //
248  if (this == &rhs) { return *this; }
249 
250  // Copy base class data
251  //
252  G4VSolid::operator=(rhs);
253 
254  // Copy data
255  //
256  fTheta = rhs.fTheta; fPhi = rhs.fPhi;
257  fDy1= rhs.fDy1; fDx1= rhs.fDx1; fDx2= rhs.fDx2; fDy2= rhs.fDy2;
258  fDx3= rhs.fDx3; fDx4= rhs.fDx4; fDz= rhs.fDz; fDx= rhs.fDx; fDy= rhs.fDy;
259  fAlph= rhs.fAlph; fTAlph= rhs.fTAlph; fdeltaX= rhs.fdeltaX;
260  fdeltaY= rhs.fdeltaY; fPhiTwist= rhs.fPhiTwist; fLowerEndcap= 0;
261  fUpperEndcap= 0; fSide0= 0; fSide90= 0; fSide180= 0; fSide270= 0;
262  fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
263  fpPolyhedron= 0;
264  fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
265  fLastDistanceToIn= rhs.fLastDistanceToIn;
266  fLastDistanceToOut= rhs.fLastDistanceToOut;
267  fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
268  fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
269 
270  CreateSurfaces();
271 
272  return *this;
273 }
274 
275 //=====================================================================
276 //* ComputeDimensions -------------------------------------------------
277 
279  const G4int ,
280  const G4VPhysicalVolume* )
281 {
282  G4Exception("G4VTwistedFaceted::ComputeDimensions()",
283  "GeomSolids0001", FatalException,
284  "G4VTwistedFaceted does not support Parameterisation.");
285 }
286 
287 //=====================================================================
288 //* CalculateExtent ---------------------------------------------------
289 
290 G4bool
292  const G4VoxelLimits &pVoxelLimit,
293  const G4AffineTransform &pTransform,
294  G4double &pMin,
295  G4double &pMax ) const
296 {
297  G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
298 
299  if (!pTransform.IsRotated())
300  {
301  // Special case handling for unrotated boxes
302  // Compute x/y/z mins and maxs respecting limits, with early returns
303  // if outside limits. Then switch() on pAxis
304 
305  G4double xoffset,xMin,xMax;
306  G4double yoffset,yMin,yMax;
307  G4double zoffset,zMin,zMax;
308 
309  xoffset = pTransform.NetTranslation().x() ;
310  xMin = xoffset - maxRad ;
311  xMax = xoffset + maxRad ;
312 
313  if (pVoxelLimit.IsXLimited())
314  {
315  if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance ||
316  xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance ) return false;
317  else
318  {
319  if (xMin < pVoxelLimit.GetMinXExtent())
320  {
321  xMin = pVoxelLimit.GetMinXExtent() ;
322  }
323  if (xMax > pVoxelLimit.GetMaxXExtent())
324  {
325  xMax = pVoxelLimit.GetMaxXExtent() ;
326  }
327  }
328  }
329  yoffset = pTransform.NetTranslation().y() ;
330  yMin = yoffset - maxRad ;
331  yMax = yoffset + maxRad ;
332 
333  if (pVoxelLimit.IsYLimited())
334  {
335  if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance ||
336  yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance ) return false;
337  else
338  {
339  if (yMin < pVoxelLimit.GetMinYExtent())
340  {
341  yMin = pVoxelLimit.GetMinYExtent() ;
342  }
343  if (yMax > pVoxelLimit.GetMaxYExtent())
344  {
345  yMax = pVoxelLimit.GetMaxYExtent() ;
346  }
347  }
348  }
349  zoffset = pTransform.NetTranslation().z() ;
350  zMin = zoffset - fDz ;
351  zMax = zoffset + fDz ;
352 
353  if (pVoxelLimit.IsZLimited())
354  {
355  if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance ||
356  zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance ) return false;
357  else
358  {
359  if (zMin < pVoxelLimit.GetMinZExtent())
360  {
361  zMin = pVoxelLimit.GetMinZExtent() ;
362  }
363  if (zMax > pVoxelLimit.GetMaxZExtent())
364  {
365  zMax = pVoxelLimit.GetMaxZExtent() ;
366  }
367  }
368  }
369  switch (pAxis)
370  {
371  case kXAxis:
372  pMin = xMin ;
373  pMax = xMax ;
374  break ;
375  case kYAxis:
376  pMin=yMin;
377  pMax=yMax;
378  break;
379  case kZAxis:
380  pMin=zMin;
381  pMax=zMax;
382  break;
383  default:
384  break;
385  }
386  pMin -= kCarTolerance ;
387  pMax += kCarTolerance ;
388 
389  return true;
390  }
391  else // General rotated case - create and clip mesh to boundaries
392  {
393  G4bool existsAfterClip = false ;
394  G4ThreeVectorList* vertices ;
395 
396  pMin = +kInfinity ;
397  pMax = -kInfinity ;
398 
399  // Calculate rotated vertex coordinates
400 
401  vertices = CreateRotatedVertices(pTransform) ;
402  ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
403  ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax) ;
404  ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
405 
406  if (pVoxelLimit.IsLimited(pAxis) == false)
407  {
408  if ( pMin != kInfinity || pMax != -kInfinity )
409  {
410  existsAfterClip = true ;
411 
412  // Add 2*tolerance to avoid precision troubles
413 
414  pMin -= kCarTolerance;
415  pMax += kCarTolerance;
416  }
417  }
418  else
419  {
420  G4ThreeVector clipCentre(
421  ( pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
422  ( pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
423  ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
424 
425  if ( pMin != kInfinity || pMax != -kInfinity )
426  {
427  existsAfterClip = true ;
428 
429  // Check to see if endpoints are in the solid
430 
431  clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis);
432 
433  if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
434  != kOutside)
435  {
436  pMin = pVoxelLimit.GetMinExtent(pAxis);
437  }
438  else
439  {
440  pMin -= kCarTolerance;
441  }
442  clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis);
443 
444  if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
445  != kOutside)
446  {
447  pMax = pVoxelLimit.GetMaxExtent(pAxis);
448  }
449  else
450  {
451  pMax += kCarTolerance;
452  }
453  }
454 
455  // Check for case where completely enveloping clipping volume
456  // If point inside then we are confident that the solid completely
457  // envelopes the clipping volume. Hence set min/max extents according
458  // to clipping volume extents along the specified axis.
459 
460  else if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
461  != kOutside)
462  {
463  existsAfterClip = true ;
464  pMin = pVoxelLimit.GetMinExtent(pAxis) ;
465  pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
466  }
467  }
468  delete vertices;
469  return existsAfterClip;
470  }
471 
472 
473 }
474 
476 CreateRotatedVertices(const G4AffineTransform& pTransform) const
477 {
478 
479  G4ThreeVectorList* vertices = new G4ThreeVectorList();
480 
481  if (vertices)
482  {
483  vertices->reserve(8);
484 
485  G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
486 
487  G4ThreeVector vertex0(-maxRad,-maxRad,-fDz) ;
488  G4ThreeVector vertex1(maxRad,-maxRad,-fDz) ;
489  G4ThreeVector vertex2(maxRad,maxRad,-fDz) ;
490  G4ThreeVector vertex3(-maxRad,maxRad,-fDz) ;
491  G4ThreeVector vertex4(-maxRad,-maxRad,fDz) ;
492  G4ThreeVector vertex5(maxRad,-maxRad,fDz) ;
493  G4ThreeVector vertex6(maxRad,maxRad,fDz) ;
494  G4ThreeVector vertex7(-maxRad,maxRad,fDz) ;
495 
496  vertices->push_back(pTransform.TransformPoint(vertex0));
497  vertices->push_back(pTransform.TransformPoint(vertex1));
498  vertices->push_back(pTransform.TransformPoint(vertex2));
499  vertices->push_back(pTransform.TransformPoint(vertex3));
500  vertices->push_back(pTransform.TransformPoint(vertex4));
501  vertices->push_back(pTransform.TransformPoint(vertex5));
502  vertices->push_back(pTransform.TransformPoint(vertex6));
503  vertices->push_back(pTransform.TransformPoint(vertex7));
504  }
505  else
506  {
507  DumpInfo();
508  G4Exception("G4VTwistedFaceted::CreateRotatedVertices()",
509  "GeomSolids0003", FatalException,
510  "Error in allocation of vertices. Out of memory !");
511  }
512  return vertices;
513 }
514 
515 //=====================================================================
516 //* Inside ------------------------------------------------------------
517 
519 {
520 
521  G4ThreeVector *tmpp;
522  EInside *tmpin;
523  if (fLastInside.p == p) {
524  return fLastInside.inside;
525  } else {
526  tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
527  tmpin = const_cast<EInside*>(&(fLastInside.inside));
528  tmpp->set(p.x(), p.y(), p.z());
529  }
530 
531  *tmpin = kOutside ;
532 
533  G4double phi = p.z()/(2*fDz) * fPhiTwist ; // rotate the point to z=0
534  G4double cphi = std::cos(-phi) ;
535  G4double sphi = std::sin(-phi) ;
536 
537  G4double px = p.x() + fdeltaX * ( -phi/fPhiTwist) ; // shift
538  G4double py = p.y() + fdeltaY * ( -phi/fPhiTwist) ;
539  G4double pz = p.z() ;
540 
541  G4double posx = px * cphi - py * sphi ; // rotation
542  G4double posy = px * sphi + py * cphi ;
543  G4double posz = pz ;
544 
545  G4double xMin = Xcoef(posy,phi,fTAlph) - 2*Xcoef(posy,phi,0.) ;
546  G4double xMax = Xcoef(posy,phi,fTAlph) ;
547 
548  G4double yMax = GetValueB(phi)/2. ; // b(phi)/2 is limit
549  G4double yMin = -yMax ;
550 
551 #ifdef G4TWISTDEBUG
552 
553  G4cout << "inside called: p = " << p << G4endl ;
554  G4cout << "fDx1 = " << fDx1 << G4endl ;
555  G4cout << "fDx2 = " << fDx2 << G4endl ;
556  G4cout << "fDx3 = " << fDx3 << G4endl ;
557  G4cout << "fDx4 = " << fDx4 << G4endl ;
558 
559  G4cout << "fDy1 = " << fDy1 << G4endl ;
560  G4cout << "fDy2 = " << fDy2 << G4endl ;
561 
562  G4cout << "fDz = " << fDz << G4endl ;
563 
564  G4cout << "Tilt angle alpha = " << fAlph << G4endl ;
565  G4cout << "phi,theta = " << fPhi << " , " << fTheta << G4endl ;
566 
567  G4cout << "Twist angle = " << fPhiTwist << G4endl ;
568 
569  G4cout << "posx = " << posx << G4endl ;
570  G4cout << "posy = " << posy << G4endl ;
571  G4cout << "xMin = " << xMin << G4endl ;
572  G4cout << "xMax = " << xMax << G4endl ;
573  G4cout << "yMin = " << yMin << G4endl ;
574  G4cout << "yMax = " << yMax << G4endl ;
575 
576 #endif
577 
578 
579  if ( posx <= xMax - kCarTolerance*0.5
580  && posx >= xMin + kCarTolerance*0.5 )
581  {
582  if ( posy <= yMax - kCarTolerance*0.5
583  && posy >= yMin + kCarTolerance*0.5 )
584  {
585  if (std::fabs(posz) <= fDz - kCarTolerance*0.5 ) *tmpin = kInside ;
586  else if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) *tmpin = kSurface ;
587  }
588  else if ( posy <= yMax + kCarTolerance*0.5
589  && posy >= yMin - kCarTolerance*0.5 )
590  {
591  if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) *tmpin = kSurface ;
592  }
593  }
594  else if ( posx <= xMax + kCarTolerance*0.5
595  && posx >= xMin - kCarTolerance*0.5 )
596  {
597  if ( posy <= yMax + kCarTolerance*0.5
598  && posy >= yMin - kCarTolerance*0.5 )
599  {
600  if (std::fabs(posz) <= fDz + kCarTolerance*0.5) *tmpin = kSurface ;
601  }
602  }
603 
604 #ifdef G4TWISTDEBUG
605  G4cout << "inside = " << fLastInside.inside << G4endl ;
606 #endif
607 
608  return fLastInside.inside;
609 
610 }
611 
612 //=====================================================================
613 //* SurfaceNormal -----------------------------------------------------
614 
616 {
617  //
618  // return the normal unit vector to the Hyperbolical Surface at a point
619  // p on (or nearly on) the surface
620  //
621  // Which of the three or four surfaces are we closest to?
622  //
623 
624  if (fLastNormal.p == p)
625  {
626  return fLastNormal.vec;
627  }
628 
629  G4ThreeVector *tmpp = const_cast<G4ThreeVector*>(&(fLastNormal.p));
630  G4ThreeVector *tmpnormal = const_cast<G4ThreeVector*>(&(fLastNormal.vec));
631  G4VTwistSurface **tmpsurface = const_cast<G4VTwistSurface**>(fLastNormal.surface);
632  tmpp->set(p.x(), p.y(), p.z());
633 
634  G4double distance = kInfinity;
635 
636  G4VTwistSurface *surfaces[6];
637 
638  surfaces[0] = fSide0 ;
639  surfaces[1] = fSide90 ;
640  surfaces[2] = fSide180 ;
641  surfaces[3] = fSide270 ;
642  surfaces[4] = fLowerEndcap;
643  surfaces[5] = fUpperEndcap;
644 
646  G4ThreeVector bestxx;
647  G4int i;
648  G4int besti = -1;
649  for (i=0; i< 6; i++)
650  {
651  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
652  if (tmpdistance < distance)
653  {
654  distance = tmpdistance;
655  bestxx = xx;
656  besti = i;
657  }
658  }
659 
660  tmpsurface[0] = surfaces[besti];
661  *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
662 
663  return fLastNormal.vec;
664 }
665 
666 //=====================================================================
667 //* DistanceToIn (p, v) -----------------------------------------------
668 
670  const G4ThreeVector& v ) const
671 {
672 
673  // DistanceToIn (p, v):
674  // Calculate distance to surface of shape from `outside'
675  // along with the v, allowing for tolerance.
676  // The function returns kInfinity if no intersection or
677  // just grazing within tolerance.
678 
679  //
680  // checking last value
681  //
682 
683  G4ThreeVector *tmpp;
684  G4ThreeVector *tmpv;
685  G4double *tmpdist;
686  if (fLastDistanceToInWithV.p == p && fLastDistanceToInWithV.vec == v)
687  {
688  return fLastDistanceToIn.value;
689  }
690  else
691  {
692  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
693  tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
694  tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
695  tmpp->set(p.x(), p.y(), p.z());
696  tmpv->set(v.x(), v.y(), v.z());
697  }
698 
699  //
700  // Calculate DistanceToIn(p,v)
701  //
702 
703  EInside currentside = Inside(p);
704 
705  if (currentside == kInside)
706  {
707  }
708  else if (currentside == kSurface)
709  {
710  // particle is just on a boundary.
711  // if the particle is entering to the volume, return 0
712  //
713  G4ThreeVector normal = SurfaceNormal(p);
714  if (normal*v < 0)
715  {
716  *tmpdist = 0;
717  return fLastDistanceToInWithV.value;
718  }
719  }
720 
721  // now, we can take smallest positive distance.
722 
723  // Initialize
724  //
725  G4double distance = kInfinity;
726 
727  // Find intersections and choose nearest one
728  //
729  G4VTwistSurface *surfaces[6];
730 
731  surfaces[0] = fSide0;
732  surfaces[1] = fSide90 ;
733  surfaces[2] = fSide180 ;
734  surfaces[3] = fSide270 ;
735  surfaces[4] = fLowerEndcap;
736  surfaces[5] = fUpperEndcap;
737 
739  G4ThreeVector bestxx;
740  G4int i;
741  for (i=0; i < 6 ; i++)
742  {
743 #ifdef G4TWISTDEBUG
744  G4cout << G4endl << "surface " << i << ": " << G4endl << G4endl ;
745 #endif
746  G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
747 #ifdef G4TWISTDEBUG
748  G4cout << "Solid DistanceToIn : distance = " << tmpdistance << G4endl ;
749  G4cout << "intersection point = " << xx << G4endl ;
750 #endif
751  if (tmpdistance < distance)
752  {
753  distance = tmpdistance;
754  bestxx = xx;
755  }
756  }
757 
758 #ifdef G4TWISTDEBUG
759  G4cout << "best distance = " << distance << G4endl ;
760 #endif
761 
762  *tmpdist = distance;
763  // timer.Stop();
764  return fLastDistanceToInWithV.value;
765 }
766 
767 
769 {
770  // DistanceToIn(p):
771  // Calculate distance to surface of shape from `outside',
772  // allowing for tolerance
773  //
774 
775  //
776  // checking last value
777  //
778 
779  G4ThreeVector *tmpp;
780  G4double *tmpdist;
781  if (fLastDistanceToIn.p == p)
782  {
783  return fLastDistanceToIn.value;
784  }
785  else
786  {
787  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
788  tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
789  tmpp->set(p.x(), p.y(), p.z());
790  }
791 
792  //
793  // Calculate DistanceToIn(p)
794  //
795 
796  EInside currentside = Inside(p);
797 
798  switch (currentside)
799  {
800  case (kInside) :
801  {
802  }
803 
804  case (kSurface) :
805  {
806  *tmpdist = 0.;
807  return fLastDistanceToIn.value;
808  }
809 
810  case (kOutside) :
811  {
812  // Initialize
813  //
814  G4double distance = kInfinity;
815 
816  // Find intersections and choose nearest one
817  //
818  G4VTwistSurface *surfaces[6];
819 
820  surfaces[0] = fSide0;
821  surfaces[1] = fSide90 ;
822  surfaces[2] = fSide180 ;
823  surfaces[3] = fSide270 ;
824  surfaces[4] = fLowerEndcap;
825  surfaces[5] = fUpperEndcap;
826 
827  G4int i;
829  G4ThreeVector bestxx;
830  for (i=0; i< 6; i++)
831  {
832  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
833  if (tmpdistance < distance)
834  {
835  distance = tmpdistance;
836  bestxx = xx;
837  }
838  }
839  *tmpdist = distance;
840  return fLastDistanceToIn.value;
841  }
842 
843  default :
844  {
845  G4Exception("G4VTwistedFaceted::DistanceToIn(p)", "GeomSolids0003",
846  FatalException, "Unknown point location!");
847  }
848  } // switch end
849 
850  return 0;
851 }
852 
853 
854 //=====================================================================
855 //* DistanceToOut (p, v) ----------------------------------------------
856 
857 G4double
859  const G4ThreeVector& v,
860  const G4bool calcNorm,
861  G4bool *validNorm,
862  G4ThreeVector *norm ) const
863 {
864  // DistanceToOut (p, v):
865  // Calculate distance to surface of shape from `inside'
866  // along with the v, allowing for tolerance.
867  // The function returns kInfinity if no intersection or
868  // just grazing within tolerance.
869 
870  //
871  // checking last value
872  //
873 
874  G4ThreeVector *tmpp;
875  G4ThreeVector *tmpv;
876  G4double *tmpdist;
877  if (fLastDistanceToOutWithV.p == p && fLastDistanceToOutWithV.vec == v )
878  {
879  return fLastDistanceToOutWithV.value;
880  }
881  else
882  {
883  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
884  tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
885  tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
886  tmpp->set(p.x(), p.y(), p.z());
887  tmpv->set(v.x(), v.y(), v.z());
888  }
889 
890  //
891  // Calculate DistanceToOut(p,v)
892  //
893 
894  EInside currentside = Inside(p);
895 
896  if (currentside == kOutside)
897  {
898  }
899  else if (currentside == kSurface)
900  {
901  // particle is just on a boundary.
902  // if the particle is exiting from the volume, return 0
903  //
904  G4ThreeVector normal = SurfaceNormal(p);
905  G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
906  if (normal*v > 0)
907  {
908  if (calcNorm)
909  {
910  *norm = (blockedsurface->GetNormal(p, true));
911  *validNorm = blockedsurface->IsValidNorm();
912  }
913  *tmpdist = 0.;
914  // timer.Stop();
915  return fLastDistanceToOutWithV.value;
916  }
917  }
918 
919  // now, we can take smallest positive distance.
920 
921  // Initialize
922  G4double distance = kInfinity;
923 
924  // find intersections and choose nearest one.
925  G4VTwistSurface *surfaces[6];
926 
927  surfaces[0] = fSide0;
928  surfaces[1] = fSide90 ;
929  surfaces[2] = fSide180 ;
930  surfaces[3] = fSide270 ;
931  surfaces[4] = fLowerEndcap;
932  surfaces[5] = fUpperEndcap;
933 
934  G4int i;
935  G4int besti = -1;
937  G4ThreeVector bestxx;
938  for (i=0; i< 6 ; i++) {
939  G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
940  if (tmpdistance < distance)
941  {
942  distance = tmpdistance;
943  bestxx = xx;
944  besti = i;
945  }
946  }
947 
948  if (calcNorm)
949  {
950  if (besti != -1)
951  {
952  *norm = (surfaces[besti]->GetNormal(p, true));
953  *validNorm = surfaces[besti]->IsValidNorm();
954  }
955  }
956 
957  *tmpdist = distance;
958  // timer.Stop();
959  return fLastDistanceToOutWithV.value;
960 }
961 
962 
964 {
965  // DistanceToOut(p):
966  // Calculate distance to surface of shape from `inside',
967  // allowing for tolerance
968 
969  //
970  // checking last value
971  //
972 
973  G4ThreeVector *tmpp;
974  G4double *tmpdist;
975 
976  if (fLastDistanceToOut.p == p)
977  {
978  return fLastDistanceToOut.value;
979  }
980  else
981  {
982  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
983  tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
984  tmpp->set(p.x(), p.y(), p.z());
985  }
986 
987  //
988  // Calculate DistanceToOut(p)
989  //
990 
991  EInside currentside = Inside(p);
992  G4double retval = kInfinity;
993 
994  switch (currentside)
995  {
996  case (kOutside) :
997  {
998 #ifdef G4SPECSDEBUG
999  G4int oldprc = G4cout.precision(16) ;
1000  G4cout << G4endl ;
1001  DumpInfo();
1002  G4cout << "Position:" << G4endl << G4endl ;
1003  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1004  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1005  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1006  G4cout.precision(oldprc) ;
1007  G4Exception("G4VTwistedFaceted::DistanceToOut(p)", "GeomSolids1002",
1008  JustWarning, "Point p is outside !?" );
1009 #endif
1010  break;
1011  }
1012  case (kSurface) :
1013  {
1014  *tmpdist = 0.;
1015  retval = fLastDistanceToOut.value;
1016  break;
1017  }
1018 
1019  case (kInside) :
1020  {
1021  // Initialize
1022  //
1023  G4double distance = kInfinity;
1024 
1025  // find intersections and choose nearest one
1026  //
1027  G4VTwistSurface *surfaces[6];
1028 
1029  surfaces[0] = fSide0;
1030  surfaces[1] = fSide90 ;
1031  surfaces[2] = fSide180 ;
1032  surfaces[3] = fSide270 ;
1033  surfaces[4] = fLowerEndcap;
1034  surfaces[5] = fUpperEndcap;
1035 
1036  G4int i;
1037  G4ThreeVector xx;
1038  G4ThreeVector bestxx;
1039  for (i=0; i< 6; i++)
1040  {
1041  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
1042  if (tmpdistance < distance)
1043  {
1044  distance = tmpdistance;
1045  bestxx = xx;
1046  }
1047  }
1048  *tmpdist = distance;
1049 
1050  retval = fLastDistanceToOut.value;
1051  break;
1052  }
1053 
1054  default :
1055  {
1056  G4Exception("G4VTwistedFaceted::DistanceToOut(p)", "GeomSolids0003",
1057  FatalException, "Unknown point location!");
1058  break;
1059  }
1060  } // switch end
1061 
1062  return retval;
1063 }
1064 
1065 
1066 //=====================================================================
1067 //* StreamInfo --------------------------------------------------------
1068 
1069 std::ostream& G4VTwistedFaceted::StreamInfo(std::ostream& os) const
1070 {
1071  //
1072  // Stream object contents to an output stream
1073  //
1074  G4int oldprc = os.precision(16);
1075  os << "-----------------------------------------------------------\n"
1076  << " *** Dump for solid - " << GetName() << " ***\n"
1077  << " ===================================================\n"
1078  << " Solid type: G4VTwistedFaceted\n"
1079  << " Parameters: \n"
1080  << " polar angle theta = " << fTheta/degree << " deg" << G4endl
1081  << " azimuthal angle phi = " << fPhi/degree << " deg" << G4endl
1082  << " tilt angle alpha = " << fAlph/degree << " deg" << G4endl
1083  << " TWIST angle = " << fPhiTwist/degree << " deg" << G4endl
1084  << " Half length along y (lower endcap) = " << fDy1/cm << " cm"
1085  << G4endl
1086  << " Half length along x (lower endcap, bottom) = " << fDx1/cm << " cm"
1087  << G4endl
1088  << " Half length along x (lower endcap, top) = " << fDx2/cm << " cm"
1089  << G4endl
1090  << " Half length along y (upper endcap) = " << fDy2/cm << " cm"
1091  << G4endl
1092  << " Half length along x (upper endcap, bottom) = " << fDx3/cm << " cm"
1093  << G4endl
1094  << " Half length along x (upper endcap, top) = " << fDx4/cm << " cm"
1095  << G4endl
1096  << "-----------------------------------------------------------\n";
1097  os.precision(oldprc);
1098 
1099  return os;
1100 }
1101 
1102 
1103 //=====================================================================
1104 //* DiscribeYourselfTo ------------------------------------------------
1105 
1107 {
1108  scene.AddSolid (*this);
1109 }
1110 
1111 //=====================================================================
1112 //* GetExtent ---------------------------------------------------------
1113 
1115 {
1116  G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
1117 
1118  return G4VisExtent(-maxRad, maxRad ,
1119  -maxRad, maxRad ,
1120  -fDz, fDz );
1121 }
1122 
1123 
1124 //=====================================================================
1125 //* CreateSurfaces ----------------------------------------------------
1126 
1127 void G4VTwistedFaceted::CreateSurfaces()
1128 {
1129 
1130  // create 6 surfaces of TwistedTub.
1131 
1132  if ( fDx1 == fDx2 && fDx3 == fDx4 ) // special case : Box
1133  {
1134  fSide0 = new G4TwistBoxSide("0deg", fPhiTwist, fDz, fTheta, fPhi,
1135  fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 0.*deg);
1136  fSide180 = new G4TwistBoxSide("180deg", fPhiTwist, fDz, fTheta, fPhi+pi,
1137  fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 180.*deg);
1138  }
1139  else // default general case
1140  {
1141  fSide0 = new G4TwistTrapAlphaSide("0deg" ,fPhiTwist, fDz, fTheta,
1142  fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
1143  fSide180 = new G4TwistTrapAlphaSide("180deg", fPhiTwist, fDz, fTheta,
1144  fPhi+pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
1145  }
1146 
1147  // create parallel sides
1148  //
1149  fSide90 = new G4TwistTrapParallelSide("90deg", fPhiTwist, fDz, fTheta,
1150  fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
1151  fSide270 = new G4TwistTrapParallelSide("270deg", fPhiTwist, fDz, fTheta,
1152  fPhi+pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
1153 
1154  // create endcaps
1155  //
1156  fUpperEndcap = new G4TwistTrapFlatSide("UpperCap",fPhiTwist, fDx3, fDx4, fDy2,
1157  fDz, fAlph, fPhi, fTheta, 1 );
1158  fLowerEndcap = new G4TwistTrapFlatSide("LowerCap",fPhiTwist, fDx1, fDx2, fDy1,
1159  fDz, fAlph, fPhi, fTheta, -1 );
1160 
1161  // Set neighbour surfaces
1162 
1163  fSide0->SetNeighbours( fSide270 , fLowerEndcap , fSide90 , fUpperEndcap );
1164  fSide90->SetNeighbours( fSide0 , fLowerEndcap , fSide180 , fUpperEndcap );
1165  fSide180->SetNeighbours(fSide90 , fLowerEndcap , fSide270 , fUpperEndcap );
1166  fSide270->SetNeighbours(fSide180 , fLowerEndcap , fSide0 , fUpperEndcap );
1167  fUpperEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
1168  fLowerEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
1169 
1170 }
1171 
1172 
1173 //=====================================================================
1174 //* GetEntityType -----------------------------------------------------
1175 
1177 {
1178  return G4String("G4VTwistedFaceted");
1179 }
1180 
1181 
1182 //=====================================================================
1183 //* GetPolyhedron -----------------------------------------------------
1184 
1186 {
1187  if (!fpPolyhedron ||
1188  fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1189  fpPolyhedron->GetNumberOfRotationSteps())
1190  {
1191  delete fpPolyhedron;
1192  fpPolyhedron = CreatePolyhedron();
1193  }
1194 
1195  return fpPolyhedron;
1196 }
1197 
1198 
1199 //=====================================================================
1200 //* GetPointInSolid ---------------------------------------------------
1201 
1203 {
1204 
1205 
1206  // this routine is only used for a test
1207  // can be deleted ...
1208 
1209  if ( z == fDz ) z -= 0.1*fDz ;
1210  if ( z == -fDz ) z += 0.1*fDz ;
1211 
1212  G4double phi = z/(2*fDz)*fPhiTwist ;
1213 
1214  return G4ThreeVector(fdeltaX * phi/fPhiTwist, fdeltaY * phi/fPhiTwist, z ) ;
1215 }
1216 
1217 
1218 //=====================================================================
1219 //* GetPointOnSurface -------------------------------------------------
1220 
1222 {
1223 
1224  G4double phi = G4RandFlat::shoot(-fPhiTwist/2.,fPhiTwist/2.);
1225  G4double u , umin, umax ; // variable for twisted surfaces
1226  G4double y ; // variable for flat surface (top and bottom)
1227 
1228  // Compute the areas. Attention: Only correct for trapezoids
1229  // where the twisting is done along the z-axis. In the general case
1230  // the computed surface area is more difficult. However this simplification
1231  // does not affect the tracking through the solid.
1232 
1233  G4double a1 = fSide0->GetSurfaceArea();
1234  G4double a2 = fSide90->GetSurfaceArea();
1235  G4double a3 = fSide180->GetSurfaceArea() ;
1236  G4double a4 = fSide270->GetSurfaceArea() ;
1237  G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1238  G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1239 
1240 #ifdef G4TWISTDEBUG
1241  G4cout << "Surface 0 deg = " << a1 << G4endl ;
1242  G4cout << "Surface 90 deg = " << a2 << G4endl ;
1243  G4cout << "Surface 180 deg = " << a3 << G4endl ;
1244  G4cout << "Surface 270 deg = " << a4 << G4endl ;
1245  G4cout << "Surface Lower = " << a5 << G4endl ;
1246  G4cout << "Surface Upper = " << a6 << G4endl ;
1247 #endif
1248 
1249  G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1250 
1251  if(chose < a1)
1252  {
1253 
1254  umin = fSide0->GetBoundaryMin(phi) ;
1255  umax = fSide0->GetBoundaryMax(phi) ;
1256  u = G4RandFlat::shoot(umin,umax) ;
1257 
1258  return fSide0->SurfacePoint(phi, u, true) ; // point on 0deg surface
1259  }
1260 
1261  else if( (chose >= a1) && (chose < a1 + a2 ) )
1262  {
1263 
1264  umin = fSide90->GetBoundaryMin(phi) ;
1265  umax = fSide90->GetBoundaryMax(phi) ;
1266 
1267  u = G4RandFlat::shoot(umin,umax) ;
1268 
1269  return fSide90->SurfacePoint(phi, u, true); // point on 90deg surface
1270  }
1271 
1272  else if( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1273  {
1274 
1275  umin = fSide180->GetBoundaryMin(phi) ;
1276  umax = fSide180->GetBoundaryMax(phi) ;
1277  u = G4RandFlat::shoot(umin,umax) ;
1278 
1279  return fSide180->SurfacePoint(phi, u, true); // point on 180 deg surface
1280  }
1281 
1282  else if( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1283  {
1284 
1285  umin = fSide270->GetBoundaryMin(phi) ;
1286  umax = fSide270->GetBoundaryMax(phi) ;
1287  u = G4RandFlat::shoot(umin,umax) ;
1288 
1289  return fSide270->SurfacePoint(phi, u, true); // point on 270 deg surface
1290  }
1291 
1292  else if( (chose >= a1 + a2 + a3 + a4 ) && (chose < a1 + a2 + a3 + a4 + a5 ) )
1293  {
1294 
1295  y = G4RandFlat::shoot(-fDy1,fDy1) ;
1296  umin = fLowerEndcap->GetBoundaryMin(y) ;
1297  umax = fLowerEndcap->GetBoundaryMax(y) ;
1298  u = G4RandFlat::shoot(umin,umax) ;
1299 
1300  return fLowerEndcap->SurfacePoint(u,y,true); // point on lower endcap
1301  }
1302  else {
1303 
1304  y = G4RandFlat::shoot(-fDy2,fDy2) ;
1305  umin = fUpperEndcap->GetBoundaryMin(y) ;
1306  umax = fUpperEndcap->GetBoundaryMax(y) ;
1307  u = G4RandFlat::shoot(umin,umax) ;
1308 
1309  return fUpperEndcap->SurfacePoint(u,y,true) ; // point on upper endcap
1310 
1311  }
1312 }
1313 
1314 
1315 //=====================================================================
1316 //* CreatePolyhedron --------------------------------------------------
1317 
1319 {
1320  // number of meshes
1321  const G4int k =
1322  G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2;
1323  const G4int n = k;
1324 
1325  const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
1326  const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
1327 
1328  G4Polyhedron *ph=new G4Polyhedron;
1329  typedef G4double G4double3[3];
1330  typedef G4int G4int4[4];
1331  G4double3* xyz = new G4double3[nnodes]; // number of nodes
1332  G4int4* faces = new G4int4[nfaces] ; // number of faces
1333 
1334  fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
1335  fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
1336  fSide270->GetFacets(k,n,xyz,faces,2) ;
1337  fSide0->GetFacets(k,n,xyz,faces,3) ;
1338  fSide90->GetFacets(k,n,xyz,faces,4) ;
1339  fSide180->GetFacets(k,n,xyz,faces,5) ;
1340 
1341  ph->createPolyhedron(nnodes,nfaces,xyz,faces);
1342 
1343  return ph;
1344 }
void set(double x, double y, double z)
G4String GetName() const
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:347
ThreeVector shoot(const G4int Ap, const G4int Af)
virtual void ComputeDimensions(G4VPVParameterisation *, const G4int, const G4VPhysicalVolume *)
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
tuple a
Definition: test.py:11
G4double GetMinYExtent() const
G4VTwistedFaceted(const G4String &pname, G4double PhiTwist, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlph)
CLHEP::Hep3Vector G4ThreeVector
double x() const
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal)=0
G4AffineTransform Inverse() const
G4bool IsYLimited() const
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
virtual G4GeometryType GetEntityType() const
const char * p
Definition: xmltok.h:285
G4ThreeVector GetPointOnSurface() const
G4ThreeVector GetPointInSolid(G4double z) const
Float_t norm
G4bool IsRotated() const
virtual G4Polyhedron * CreatePolyhedron() const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Double_t xx
Definition: macro.C:10
G4ThreeVector NetTranslation() const
G4bool IsXLimited() const
G4VTwistedFaceted & operator=(const G4VTwistedFaceted &rhs)
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4bool IsValidNorm() const
double z() const
void DumpInfo() const
virtual G4double DistanceToIn(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
G4double Xcoef(G4double u, G4double phi, G4double ftg) const
Double_t y
Definition: plot.C:279
G4double GetMaxXExtent() const
virtual G4double GetSurfaceArea()=0
Char_t n[5]
G4double GetMinZExtent() const
G4GLOB_DLL std::ostream G4cout
G4bool IsLimited() const
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
void SetNeighbours(G4VTwistSurface *axis0min, G4VTwistSurface *axis1min, G4VTwistSurface *axis0max, G4VTwistSurface *axis1max)
virtual G4double DistanceTo(const G4ThreeVector &gp, G4ThreeVector &gxx)
tuple degree
Definition: hepunit.py:69
virtual G4double DistanceToOut(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)=0
bool G4bool
Definition: G4Types.hh:79
G4double GetValueB(G4double phi) const
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
virtual G4double GetBoundaryMax(G4double)=0
tuple v
Definition: test.py:18
string pname
Definition: eplot.py:33
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)=0
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcnorm=false, G4bool *validnorm=0, G4ThreeVector *n=0) const
virtual G4double GetBoundaryMin(G4double)=0
G4double GetMinXExtent() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double GetMaxZExtent() const
virtual G4Polyhedron * GetPolyhedron() const
static G4int GetNumberOfRotationSteps()
EInside
Definition: geomdefs.hh:58
virtual G4VisExtent GetExtent() const
virtual EInside Inside(const G4ThreeVector &p) const
EAxis
Definition: geomdefs.hh:54
double y() const
virtual void DescribeYourselfTo(G4VGraphicsScene &scene) const
tuple z
Definition: test.py:28
#define G4endl
Definition: G4ios.hh:61
G4double GetMaxYExtent() const
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:110
G4double kCarTolerance
Definition: G4VSolid.hh:305
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
double G4double
Definition: G4Types.hh:76
virtual std::ostream & StreamInfo(std::ostream &os) const
G4double GetMaxExtent(const EAxis pAxis) const
G4bool IsZLimited() const
G4double GetAngularTolerance() const
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:378
G4double GetMinExtent(const EAxis pAxis) const
static G4GeometryTolerance * GetInstance()