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