Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TwistedTubs.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: G4TwistedTubs.cc 67011 2013-01-29 16:17:41Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4TwistTubsSide.cc
35 //
36 // Author:
37 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
38 //
39 // History:
40 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
41 // from original version in Jupiter-2.5.02 application.
42 // --------------------------------------------------------------------
43 
44 #include "G4TwistedTubs.hh"
45 
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4GeometryTolerance.hh"
49 #include "G4VoxelLimits.hh"
50 #include "G4AffineTransform.hh"
51 #include "G4SolidExtentList.hh"
52 #include "G4ClippablePolygon.hh"
53 #include "G4VPVParameterisation.hh"
54 #include "meshdefs.hh"
55 
56 #include "G4VGraphicsScene.hh"
57 #include "G4Polyhedron.hh"
58 #include "G4VisExtent.hh"
59 #include "G4NURBS.hh"
60 #include "G4NURBStube.hh"
61 #include "G4NURBScylinder.hh"
62 #include "G4NURBStubesector.hh"
63 
64 #include "Randomize.hh"
65 
66 //=====================================================================
67 //* constructors ------------------------------------------------------
68 
69 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
70  G4double twistedangle,
71  G4double endinnerrad,
72  G4double endouterrad,
73  G4double halfzlen,
74  G4double dphi)
75  : G4VSolid(pname), fDPhi(dphi),
76  fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
77  fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
78  fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
79 {
80  if (endinnerrad < DBL_MIN)
81  {
82  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
83  FatalErrorInArgument, "Invalid end-inner-radius!");
84  }
85 
86  G4double sinhalftwist = std::sin(0.5 * twistedangle);
87 
88  G4double endinnerradX = endinnerrad * sinhalftwist;
89  G4double innerrad = std::sqrt( endinnerrad * endinnerrad
90  - endinnerradX * endinnerradX );
91 
92  G4double endouterradX = endouterrad * sinhalftwist;
93  G4double outerrad = std::sqrt( endouterrad * endouterrad
94  - endouterradX * endouterradX );
95 
96  // temporary treatment!!
97  SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
98  CreateSurfaces();
99 }
100 
101 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
102  G4double twistedangle,
103  G4double endinnerrad,
104  G4double endouterrad,
105  G4double halfzlen,
106  G4int nseg,
107  G4double totphi)
108  : G4VSolid(pname),
109  fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
110  fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
111  fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
112 {
113 
114  if (!nseg)
115  {
116  std::ostringstream message;
117  message << "Invalid number of segments." << G4endl
118  << " nseg = " << nseg;
119  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
120  FatalErrorInArgument, message);
121  }
122  if (totphi == DBL_MIN || endinnerrad < DBL_MIN)
123  {
124  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
125  FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
126  }
127 
128  G4double sinhalftwist = std::sin(0.5 * twistedangle);
129 
130  G4double endinnerradX = endinnerrad * sinhalftwist;
131  G4double innerrad = std::sqrt( endinnerrad * endinnerrad
132  - endinnerradX * endinnerradX );
133 
134  G4double endouterradX = endouterrad * sinhalftwist;
135  G4double outerrad = std::sqrt( endouterrad * endouterrad
136  - endouterradX * endouterradX );
137 
138  // temporary treatment!!
139  fDPhi = totphi / nseg;
140  SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
141  CreateSurfaces();
142 }
143 
144 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
145  G4double twistedangle,
146  G4double innerrad,
147  G4double outerrad,
148  G4double negativeEndz,
149  G4double positiveEndz,
150  G4double dphi)
151  : G4VSolid(pname), fDPhi(dphi),
152  fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
153  fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
154  fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
155 {
156  if (innerrad < DBL_MIN)
157  {
158  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
159  FatalErrorInArgument, "Invalid end-inner-radius!");
160  }
161 
162  SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
163  CreateSurfaces();
164 }
165 
166 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
167  G4double twistedangle,
168  G4double innerrad,
169  G4double outerrad,
170  G4double negativeEndz,
171  G4double positiveEndz,
172  G4int nseg,
173  G4double totphi)
174  : G4VSolid(pname),
175  fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
176  fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
177  fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
178 {
179  if (!nseg)
180  {
181  std::ostringstream message;
182  message << "Invalid number of segments." << G4endl
183  << " nseg = " << nseg;
184  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
185  FatalErrorInArgument, message);
186  }
187  if (totphi == DBL_MIN || innerrad < DBL_MIN)
188  {
189  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
190  FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
191  }
192 
193  fDPhi = totphi / nseg;
194  SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
195  CreateSurfaces();
196 }
197 
198 //=====================================================================
199 //* Fake default constructor ------------------------------------------
200 
201 G4TwistedTubs::G4TwistedTubs( __void__& a )
202  : G4VSolid(a), fPhiTwist(0.), fInnerRadius(0.), fOuterRadius(0.), fDPhi(0.),
203  fZHalfLength(0.), fInnerStereo(0.), fOuterStereo(0.), fTanInnerStereo(0.),
204  fTanOuterStereo(0.), fKappa(0.), fInnerRadius2(0.), fOuterRadius2(0.),
205  fTanInnerStereo2(0.), fTanOuterStereo2(0.), fLowerEndcap(0), fUpperEndcap(0),
206  fLatterTwisted(0), fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
207  fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
208 {
209  fEndZ[0] = 0.; fEndZ[1] = 0.;
210  fEndInnerRadius[0] = 0.; fEndInnerRadius[1] = 0.;
211  fEndOuterRadius[0] = 0.; fEndOuterRadius[1] = 0.;
212  fEndPhi[0] = 0.; fEndPhi[1] = 0.;
213  fEndZ2[0] = 0.; fEndZ2[1] = 0.;
214 }
215 
216 //=====================================================================
217 //* destructor --------------------------------------------------------
218 
219 G4TwistedTubs::~G4TwistedTubs()
220 {
221  if (fLowerEndcap) { delete fLowerEndcap; }
222  if (fUpperEndcap) { delete fUpperEndcap; }
223  if (fLatterTwisted) { delete fLatterTwisted; }
224  if (fFormerTwisted) { delete fFormerTwisted; }
225  if (fInnerHype) { delete fInnerHype; }
226  if (fOuterHype) { delete fOuterHype; }
227  if (fpPolyhedron) { delete fpPolyhedron; }
228 }
229 
230 //=====================================================================
231 //* Copy constructor --------------------------------------------------
232 
233 G4TwistedTubs::G4TwistedTubs(const G4TwistedTubs& rhs)
234  : G4VSolid(rhs), fPhiTwist(rhs.fPhiTwist),
235  fInnerRadius(rhs.fInnerRadius), fOuterRadius(rhs.fOuterRadius),
236  fDPhi(rhs.fDPhi), fZHalfLength(rhs.fZHalfLength),
237  fInnerStereo(rhs.fInnerStereo), fOuterStereo(rhs.fOuterStereo),
238  fTanInnerStereo(rhs.fTanInnerStereo), fTanOuterStereo(rhs.fTanOuterStereo),
239  fKappa(rhs.fKappa), fInnerRadius2(rhs.fInnerRadius2),
240  fOuterRadius2(rhs.fOuterRadius2), fTanInnerStereo2(rhs.fTanInnerStereo2),
241  fTanOuterStereo2(rhs.fTanOuterStereo2),
242  fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0), fFormerTwisted(0),
243  fInnerHype(0), fOuterHype(0),
244  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
245  fpPolyhedron(0), fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
246  fLastDistanceToIn(rhs.fLastDistanceToIn),
247  fLastDistanceToOut(rhs.fLastDistanceToOut),
248  fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
249  fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
250 {
251  for (size_t i=0; i<2; ++i)
252  {
253  fEndZ[i] = rhs.fEndZ[i];
254  fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
255  fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
256  fEndPhi[i] = rhs.fEndPhi[i];
257  fEndZ2[i] = rhs.fEndZ2[i];
258  }
259  CreateSurfaces();
260 }
261 
262 
263 //=====================================================================
264 //* Assignment operator -----------------------------------------------
265 
266 G4TwistedTubs& G4TwistedTubs::operator = (const G4TwistedTubs& rhs)
267 {
268  // Check assignment to self
269  //
270  if (this == &rhs) { return *this; }
271 
272  // Copy base class data
273  //
274  G4VSolid::operator=(rhs);
275 
276  // Copy data
277  //
278  fPhiTwist= rhs.fPhiTwist;
279  fInnerRadius= rhs.fInnerRadius; fOuterRadius= rhs.fOuterRadius;
280  fDPhi= rhs.fDPhi; fZHalfLength= rhs.fZHalfLength;
281  fInnerStereo= rhs.fInnerStereo; fOuterStereo= rhs.fOuterStereo;
282  fTanInnerStereo= rhs.fTanInnerStereo; fTanOuterStereo= rhs.fTanOuterStereo;
283  fKappa= rhs.fKappa; fInnerRadius2= rhs.fInnerRadius2;
284  fOuterRadius2= rhs.fOuterRadius2; fTanInnerStereo2= rhs.fTanInnerStereo2;
285  fTanOuterStereo2= rhs.fTanOuterStereo2;
286  fLowerEndcap= fUpperEndcap= fLatterTwisted= fFormerTwisted= 0;
287  fInnerHype= fOuterHype= 0;
288  fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
289  fpPolyhedron= 0;
290  fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
291  fLastDistanceToIn= rhs.fLastDistanceToIn;
292  fLastDistanceToOut= rhs.fLastDistanceToOut;
293  fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
294  fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
295 
296  for (size_t i=0; i<2; ++i)
297  {
298  fEndZ[i] = rhs.fEndZ[i];
299  fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
300  fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
301  fEndPhi[i] = rhs.fEndPhi[i];
302  fEndZ2[i] = rhs.fEndZ2[i];
303  }
304 
305  CreateSurfaces();
306 
307  return *this;
308 }
309 
310 //=====================================================================
311 //* ComputeDimensions -------------------------------------------------
312 
313 void G4TwistedTubs::ComputeDimensions(G4VPVParameterisation* /* p */ ,
314  const G4int /* n */ ,
315  const G4VPhysicalVolume* /* pRep */ )
316 {
317  G4Exception("G4TwistedTubs::ComputeDimensions()",
318  "GeomSolids0001", FatalException,
319  "G4TwistedTubs does not support Parameterisation.");
320 }
321 
322 
323 //=====================================================================
324 //* CalculateExtent ---------------------------------------------------
325 
326 G4bool G4TwistedTubs::CalculateExtent( const EAxis axis,
327  const G4VoxelLimits &voxelLimit,
328  const G4AffineTransform &transform,
329  G4double &min,
330  G4double &max ) const
331 {
332 
333  G4SolidExtentList extentList( axis, voxelLimit );
334  G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
335  fEndOuterRadius[0] : fEndOuterRadius[1]);
336  G4double maxEndInnerRad = (fEndInnerRadius[0] > fEndInnerRadius[1] ?
337  fEndInnerRadius[0] : fEndInnerRadius[1]);
338  G4double maxphi = (std::fabs(fEndPhi[0]) > std::fabs(fEndPhi[1]) ?
339  std::fabs(fEndPhi[0]) : std::fabs(fEndPhi[1]));
340 
341  //
342  // Choose phi size of our segment(s) based on constants as
343  // defined in meshdefs.hh
344  //
345  // G4int numPhi = kMaxMeshSections;
346  G4double sigPhi = 2*maxphi + fDPhi;
347  G4double rFudge = 1.0/std::cos(0.5*sigPhi);
348  G4double fudgeEndOuterRad = rFudge * maxEndOuterRad;
349 
350  //
351  // We work around in phi building polygons along the way.
352  // As a reasonable compromise between accuracy and
353  // complexity (=cpu time), the following facets are chosen:
354  //
355  // 1. If fOuterRadius/maxEndOuterRad > 0.95, approximate
356  // the outer surface as a cylinder, and use one
357  // rectangular polygon (0-1) to build its mesh.
358  //
359  // Otherwise, use two trapazoidal polygons that
360  // meet at z = 0 (0-4-1)
361  //
362  // 2. If there is no inner surface, then use one
363  // polygon for each entire endcap. (0) and (1)
364  //
365  // Otherwise, use a trapazoidal polygon for each
366  // phi segment of each endcap. (0-2) and (1-3)
367  //
368  // 3. For the inner surface, if fInnerRadius/maxEndInnerRad > 0.95,
369  // approximate the inner surface as a cylinder of
370  // radius fInnerRadius and use one rectangular polygon
371  // to build each phi segment of its mesh. (2-3)
372  //
373  // Otherwise, use one rectangular polygon centered
374  // at z = 0 (5-6) and two connecting trapazoidal polygons
375  // for each phi segment (2-5) and (3-6).
376  //
377 
378  G4bool splitOuter = (fOuterRadius/maxEndOuterRad < 0.95);
379  G4bool splitInner = (fInnerRadius/maxEndInnerRad < 0.95);
380 
381  //
382  // Vertex assignments (v and w arrays)
383  // [0] and [1] are mandatory
384  // the rest are optional
385  //
386  // + -
387  // [0]------[4]------[1] <--- outer radius
388  // | |
389  // | |
390  // [2]---[5]---[6]---[3] <--- inner radius
391  //
392 
393  G4ClippablePolygon endPoly1, endPoly2;
394 
395  G4double phimax = maxphi + 0.5*fDPhi;
396  if ( phimax > pi/2) { phimax = pi-phimax; }
397  G4double phimin = - phimax;
398 
399  G4ThreeVector v0, v1, v2, v3, v4, v5, v6; // -ve phi verticies for polygon
400  G4ThreeVector w0, w1, w2, w3, w4, w5, w6; // +ve phi verticies for polygon
401 
402  //
403  // decide verticies of -ve phi boundary
404  //
405 
406  G4double cosPhi = std::cos(phimin);
407  G4double sinPhi = std::sin(phimin);
408 
409  // Outer hyperbolic surface
410 
411  v0 = transform.TransformPoint(
412  G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
413  + fZHalfLength));
414  v1 = transform.TransformPoint(
415  G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
416  - fZHalfLength));
417  if (splitOuter)
418  {
419  v4 = transform.TransformPoint(
420  G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 0));
421  }
422 
423  // Inner hyperbolic surface
424 
425  G4double zInnerSplit = 0.;
426  if (splitInner)
427  {
428  v2 = transform.TransformPoint(
429  G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
430  + fZHalfLength));
431  v3 = transform.TransformPoint(
432  G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
433  - fZHalfLength));
434 
435  // Find intersection of tangential line of inner
436  // surface at z = fZHalfLength and line r=fInnerRadius.
437  G4double dr = fZHalfLength * fTanInnerStereo2;
438  G4double dz = maxEndInnerRad;
439  zInnerSplit = fZHalfLength + (fInnerRadius - maxEndInnerRad) * dz / dr;
440 
441  // Build associated vertices
442  v5 = transform.TransformPoint(
443  G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
444  + zInnerSplit));
445  v6 = transform.TransformPoint(
446  G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
447  - zInnerSplit));
448  }
449  else
450  {
451  v2 = transform.TransformPoint(
452  G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
453  + fZHalfLength));
454  v3 = transform.TransformPoint(
455  G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
456  - fZHalfLength));
457  }
458 
459  //
460  // decide vertices of +ve phi boundary
461  //
462 
463  cosPhi = std::cos(phimax);
464  sinPhi = std::sin(phimax);
465 
466  // Outer hyperbolic surface
467 
468  w0 = transform.TransformPoint(
469  G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
470  + fZHalfLength));
471  w1 = transform.TransformPoint(
472  G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
473  - fZHalfLength));
474  if (splitOuter)
475  {
476  G4double r = rFudge*fOuterRadius;
477 
478  w4 = transform.TransformPoint(G4ThreeVector( r*cosPhi, r*sinPhi, 0 ));
479 
480  AddPolyToExtent( v0, v4, w4, w0, voxelLimit, axis, extentList );
481  AddPolyToExtent( v4, v1, w1, w4, voxelLimit, axis, extentList );
482  }
483  else
484  {
485  AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList );
486  }
487 
488  // Inner hyperbolic surface
489 
490  if (splitInner)
491  {
492  w2 = transform.TransformPoint(
493  G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
494  + fZHalfLength));
495  w3 = transform.TransformPoint(
496  G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
497  - fZHalfLength));
498 
499  w5 = transform.TransformPoint(
500  G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
501  + zInnerSplit));
502  w6 = transform.TransformPoint(
503  G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
504  - zInnerSplit));
505 
506  AddPolyToExtent( v3, v6, w6, w3, voxelLimit, axis, extentList );
507  AddPolyToExtent( v6, v5, w5, w6, voxelLimit, axis, extentList );
508  AddPolyToExtent( v5, v2, w2, w5, voxelLimit, axis, extentList );
509 
510  }
511  else
512  {
513  w2 = transform.TransformPoint(
514  G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
515  + fZHalfLength));
516  w3 = transform.TransformPoint(
517  G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
518  - fZHalfLength));
519 
520  AddPolyToExtent( v3, v2, w2, w3, voxelLimit, axis, extentList );
521  }
522 
523  //
524  // Endplate segments
525  //
526  AddPolyToExtent( v1, v3, w3, w1, voxelLimit, axis, extentList );
527  AddPolyToExtent( v2, v0, w0, w2, voxelLimit, axis, extentList );
528 
529  //
530  // Return min/max value
531  //
532  return extentList.GetExtent( min, max );
533 }
534 
535 
536 //=====================================================================
537 //* AddPolyToExtent ---------------------------------------------------
538 
539 void G4TwistedTubs::AddPolyToExtent( const G4ThreeVector &v0,
540  const G4ThreeVector &v1,
541  const G4ThreeVector &w1,
542  const G4ThreeVector &w0,
543  const G4VoxelLimits &voxelLimit,
544  const EAxis axis,
545  G4SolidExtentList &extentList )
546 {
547  // Utility function for CalculateExtent
548  //
549  G4ClippablePolygon phiPoly;
550 
551  phiPoly.AddVertexInOrder( v0 );
552  phiPoly.AddVertexInOrder( v1 );
553  phiPoly.AddVertexInOrder( w1 );
554  phiPoly.AddVertexInOrder( w0 );
555 
556  if (phiPoly.PartialClip( voxelLimit, axis ))
557  {
558  phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
559  extentList.AddSurface( phiPoly );
560  }
561 }
562 
563 
564 //=====================================================================
565 //* Inside ------------------------------------------------------------
566 
567 EInside G4TwistedTubs::Inside(const G4ThreeVector& p) const
568 {
569 
570  static const G4double halftol
571  = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
572  // static G4int timerid = -1;
573  // G4Timer timer(timerid, "G4TwistedTubs", "Inside");
574  // timer.Start();
575 
576 
577 
578  G4ThreeVector *tmpp;
579  EInside *tmpinside;
580  if (fLastInside.p == p)
581  {
582  return fLastInside.inside;
583  }
584  else
585  {
586  tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
587  tmpinside = const_cast<EInside*>(&(fLastInside.inside));
588  tmpp->set(p.x(), p.y(), p.z());
589  }
590 
591  EInside outerhypearea = ((G4TwistTubsHypeSide *)fOuterHype)->Inside(p);
592  G4double innerhyperho = ((G4TwistTubsHypeSide *)fInnerHype)->GetRhoAtPZ(p);
593  G4double distanceToOut = p.getRho() - innerhyperho; // +ve: inside
594 
595  if ((outerhypearea == kOutside) || (distanceToOut < -halftol))
596  {
597  *tmpinside = kOutside;
598  }
599  else if (outerhypearea == kSurface)
600  {
601  *tmpinside = kSurface;
602  }
603  else
604  {
605  if (distanceToOut <= halftol)
606  {
607  *tmpinside = kSurface;
608  }
609  else
610  {
611  *tmpinside = kInside;
612  }
613  }
614 
615  return fLastInside.inside;
616 }
617 
618 //=====================================================================
619 //* SurfaceNormal -----------------------------------------------------
620 
621 G4ThreeVector G4TwistedTubs::SurfaceNormal(const G4ThreeVector& p) const
622 {
623  //
624  // return the normal unit vector to the Hyperbolical Surface at a point
625  // p on (or nearly on) the surface
626  //
627  // Which of the three or four surfaces are we closest to?
628  //
629 
630  if (fLastNormal.p == p)
631  {
632  return fLastNormal.vec;
633  }
634  G4ThreeVector *tmpp =
635  const_cast<G4ThreeVector*>(&(fLastNormal.p));
636  G4ThreeVector *tmpnormal =
637  const_cast<G4ThreeVector*>(&(fLastNormal.vec));
638  G4VTwistSurface **tmpsurface =
639  const_cast<G4VTwistSurface**>(fLastNormal.surface);
640  tmpp->set(p.x(), p.y(), p.z());
641 
642  G4double distance = kInfinity;
643 
644  G4VTwistSurface *surfaces[6];
645  surfaces[0] = fLatterTwisted;
646  surfaces[1] = fFormerTwisted;
647  surfaces[2] = fInnerHype;
648  surfaces[3] = fOuterHype;
649  surfaces[4] = fLowerEndcap;
650  surfaces[5] = fUpperEndcap;
651 
652  G4ThreeVector xx;
653  G4ThreeVector bestxx;
654  G4int i;
655  G4int besti = -1;
656  for (i=0; i< 6; i++)
657  {
658  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
659  if (tmpdistance < distance)
660  {
661  distance = tmpdistance;
662  bestxx = xx;
663  besti = i;
664  }
665  }
666 
667  tmpsurface[0] = surfaces[besti];
668  *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
669 
670  return fLastNormal.vec;
671 }
672 
673 //=====================================================================
674 //* DistanceToIn (p, v) -----------------------------------------------
675 
676 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p,
677  const G4ThreeVector& v ) const
678 {
679 
680  // DistanceToIn (p, v):
681  // Calculate distance to surface of shape from `outside'
682  // along with the v, allowing for tolerance.
683  // The function returns kInfinity if no intersection or
684  // just grazing within tolerance.
685 
686  //
687  // checking last value
688  //
689 
690  G4ThreeVector *tmpp;
691  G4ThreeVector *tmpv;
692  G4double *tmpdist;
693  if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v))
694  {
695  return fLastDistanceToIn.value;
696  }
697  else
698  {
699  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
700  tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
701  tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
702  tmpp->set(p.x(), p.y(), p.z());
703  tmpv->set(v.x(), v.y(), v.z());
704  }
705 
706  //
707  // Calculate DistanceToIn(p,v)
708  //
709 
710  EInside currentside = Inside(p);
711 
712  if (currentside == kInside)
713  {
714  }
715  else
716  {
717  if (currentside == kSurface)
718  {
719  // particle is just on a boundary.
720  // If the particle is entering to the volume, return 0.
721  //
722  G4ThreeVector normal = SurfaceNormal(p);
723  if (normal*v < 0)
724  {
725  *tmpdist = 0;
726  return fLastDistanceToInWithV.value;
727  }
728  }
729  }
730 
731  // now, we can take smallest positive distance.
732 
733  // Initialize
734  //
735  G4double distance = kInfinity;
736 
737  // find intersections and choose nearest one.
738  //
739  G4VTwistSurface *surfaces[6];
740  surfaces[0] = fLowerEndcap;
741  surfaces[1] = fUpperEndcap;
742  surfaces[2] = fLatterTwisted;
743  surfaces[3] = fFormerTwisted;
744  surfaces[4] = fInnerHype;
745  surfaces[5] = fOuterHype;
746 
747  G4ThreeVector xx;
748  G4ThreeVector bestxx;
749  G4int i;
750  for (i=0; i< 6; i++)
751  {
752  G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
753  if (tmpdistance < distance)
754  {
755  distance = tmpdistance;
756  bestxx = xx;
757  }
758  }
759  *tmpdist = distance;
760 
761  return fLastDistanceToInWithV.value;
762 }
763 
764 //=====================================================================
765 //* DistanceToIn (p) --------------------------------------------------
766 
767 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p) const
768 {
769  // DistanceToIn(p):
770  // Calculate distance to surface of shape from `outside',
771  // allowing for tolerance
772 
773  //
774  // checking last value
775  //
776 
777  G4ThreeVector *tmpp;
778  G4double *tmpdist;
779  if (fLastDistanceToIn.p == p)
780  {
781  return fLastDistanceToIn.value;
782  }
783  else
784  {
785  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
786  tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
787  tmpp->set(p.x(), p.y(), p.z());
788  }
789 
790  //
791  // Calculate DistanceToIn(p)
792  //
793 
794  EInside currentside = Inside(p);
795 
796  switch (currentside)
797  {
798  case (kInside) :
799  {}
800  case (kSurface) :
801  {
802  *tmpdist = 0.;
803  return fLastDistanceToIn.value;
804  }
805  case (kOutside) :
806  {
807  // Initialize
808  G4double distance = kInfinity;
809 
810  // find intersections and choose nearest one.
811  G4VTwistSurface *surfaces[6];
812  surfaces[0] = fLowerEndcap;
813  surfaces[1] = fUpperEndcap;
814  surfaces[2] = fLatterTwisted;
815  surfaces[3] = fFormerTwisted;
816  surfaces[4] = fInnerHype;
817  surfaces[5] = fOuterHype;
818 
819  G4int i;
820  G4ThreeVector xx;
821  G4ThreeVector bestxx;
822  for (i=0; i< 6; i++)
823  {
824  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
825  if (tmpdistance < distance)
826  {
827  distance = tmpdistance;
828  bestxx = xx;
829  }
830  }
831  *tmpdist = distance;
832  return fLastDistanceToIn.value;
833  }
834  default :
835  {
836  G4Exception("G4TwistedTubs::DistanceToIn(p)", "GeomSolids0003",
837  FatalException, "Unknown point location!");
838  }
839  } // switch end
840 
841  return kInfinity;
842 }
843 
844 //=====================================================================
845 //* DistanceToOut (p, v) ----------------------------------------------
846 
847 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p,
848  const G4ThreeVector& v,
849  const G4bool calcNorm,
850  G4bool *validNorm,
851  G4ThreeVector *norm ) const
852 {
853  // DistanceToOut (p, v):
854  // Calculate distance to surface of shape from `inside'
855  // along with the v, allowing for tolerance.
856  // The function returns kInfinity if no intersection or
857  // just grazing within tolerance.
858 
859  //
860  // checking last value
861  //
862 
863  G4ThreeVector *tmpp;
864  G4ThreeVector *tmpv;
865  G4double *tmpdist;
866  if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) )
867  {
868  return fLastDistanceToOutWithV.value;
869  }
870  else
871  {
872  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
873  tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
874  tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
875  tmpp->set(p.x(), p.y(), p.z());
876  tmpv->set(v.x(), v.y(), v.z());
877  }
878 
879  //
880  // Calculate DistanceToOut(p,v)
881  //
882 
883  EInside currentside = Inside(p);
884 
885  if (currentside == kOutside)
886  {
887  }
888  else
889  {
890  if (currentside == kSurface)
891  {
892  // particle is just on a boundary.
893  // If the particle is exiting from the volume, return 0.
894  //
895  G4ThreeVector normal = SurfaceNormal(p);
896  G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
897  if (normal*v > 0)
898  {
899  if (calcNorm)
900  {
901  *norm = (blockedsurface->GetNormal(p, true));
902  *validNorm = blockedsurface->IsValidNorm();
903  }
904  *tmpdist = 0.;
905  return fLastDistanceToOutWithV.value;
906  }
907  }
908  }
909 
910  // now, we can take smallest positive distance.
911 
912  // Initialize
913  //
914  G4double distance = kInfinity;
915 
916  // find intersections and choose nearest one.
917  //
918  G4VTwistSurface *surfaces[6];
919  surfaces[0] = fLatterTwisted;
920  surfaces[1] = fFormerTwisted;
921  surfaces[2] = fInnerHype;
922  surfaces[3] = fOuterHype;
923  surfaces[4] = fLowerEndcap;
924  surfaces[5] = fUpperEndcap;
925 
926  G4int i;
927  G4int besti = -1;
928  G4ThreeVector xx;
929  G4ThreeVector bestxx;
930  for (i=0; i< 6; i++)
931  {
932  G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
933  if (tmpdistance < distance)
934  {
935  distance = tmpdistance;
936  bestxx = xx;
937  besti = i;
938  }
939  }
940 
941  if (calcNorm)
942  {
943  if (besti != -1)
944  {
945  *norm = (surfaces[besti]->GetNormal(p, true));
946  *validNorm = surfaces[besti]->IsValidNorm();
947  }
948  }
949 
950  *tmpdist = distance;
951 
952  return fLastDistanceToOutWithV.value;
953 }
954 
955 
956 //=====================================================================
957 //* DistanceToOut (p) ----------------------------------------------
958 
959 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p ) const
960 {
961  // DistanceToOut(p):
962  // Calculate distance to surface of shape from `inside',
963  // allowing for tolerance
964 
965  //
966  // checking last value
967  //
968 
969  G4ThreeVector *tmpp;
970  G4double *tmpdist;
971  if (fLastDistanceToOut.p == p)
972  {
973  return fLastDistanceToOut.value;
974  }
975  else
976  {
977  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
978  tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
979  tmpp->set(p.x(), p.y(), p.z());
980  }
981 
982  //
983  // Calculate DistanceToOut(p)
984  //
985 
986  EInside currentside = Inside(p);
987 
988  switch (currentside)
989  {
990  case (kOutside) :
991  {
992  }
993  case (kSurface) :
994  {
995  *tmpdist = 0.;
996  return fLastDistanceToOut.value;
997  }
998  case (kInside) :
999  {
1000  // Initialize
1001  G4double distance = kInfinity;
1002 
1003  // find intersections and choose nearest one.
1004  G4VTwistSurface *surfaces[6];
1005  surfaces[0] = fLatterTwisted;
1006  surfaces[1] = fFormerTwisted;
1007  surfaces[2] = fInnerHype;
1008  surfaces[3] = fOuterHype;
1009  surfaces[4] = fLowerEndcap;
1010  surfaces[5] = fUpperEndcap;
1011 
1012  G4int i;
1013  G4ThreeVector xx;
1014  G4ThreeVector bestxx;
1015  for (i=0; i< 6; i++)
1016  {
1017  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
1018  if (tmpdistance < distance)
1019  {
1020  distance = tmpdistance;
1021  bestxx = xx;
1022  }
1023  }
1024  *tmpdist = distance;
1025 
1026  return fLastDistanceToOut.value;
1027  }
1028  default :
1029  {
1030  G4Exception("G4TwistedTubs::DistanceToOut(p)", "GeomSolids0003",
1031  FatalException, "Unknown point location!");
1032  }
1033  } // switch end
1034 
1035  return 0;
1036 }
1037 
1038 //=====================================================================
1039 //* StreamInfo --------------------------------------------------------
1040 
1041 std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
1042 {
1043  //
1044  // Stream object contents to an output stream
1045  //
1046  G4int oldprc = os.precision(16);
1047  os << "-----------------------------------------------------------\n"
1048  << " *** Dump for solid - " << GetName() << " ***\n"
1049  << " ===================================================\n"
1050  << " Solid type: G4TwistedTubs\n"
1051  << " Parameters: \n"
1052  << " -ve end Z : " << fEndZ[0]/mm << " mm \n"
1053  << " +ve end Z : " << fEndZ[1]/mm << " mm \n"
1054  << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
1055  << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
1056  << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
1057  << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
1058  << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n"
1059  << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n"
1060  << " twisted angle : " << fPhiTwist/degree << " degrees \n"
1061  << " inner stereo angle : " << fInnerStereo/degree << " degrees \n"
1062  << " outer stereo angle : " << fOuterStereo/degree << " degrees \n"
1063  << " phi-width of a piece : " << fDPhi/degree << " degrees \n"
1064  << "-----------------------------------------------------------\n";
1065  os.precision(oldprc);
1066 
1067  return os;
1068 }
1069 
1070 
1071 //=====================================================================
1072 //* DiscribeYourselfTo ------------------------------------------------
1073 
1074 void G4TwistedTubs::DescribeYourselfTo (G4VGraphicsScene& scene) const
1075 {
1076  scene.AddSolid (*this);
1077 }
1078 
1079 //=====================================================================
1080 //* GetExtent ---------------------------------------------------------
1081 
1082 G4VisExtent G4TwistedTubs::GetExtent() const
1083 {
1084  // Define the sides of the box into which the G4Tubs instance would fit.
1085 
1086  G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
1087  return G4VisExtent( -maxEndOuterRad, maxEndOuterRad,
1088  -maxEndOuterRad, maxEndOuterRad,
1089  -fZHalfLength, fZHalfLength );
1090 }
1091 
1092 //=====================================================================
1093 //* CreatePolyhedron --------------------------------------------------
1094 
1095 G4Polyhedron* G4TwistedTubs::CreatePolyhedron () const
1096 {
1097  // number of meshes
1098  //
1099  G4double dA = std::max(fDPhi,fPhiTwist);
1100  const G4int k =
1101  G4int(G4Polyhedron::GetNumberOfRotationSteps() * dA / twopi) + 2;
1102  const G4int n =
1103  G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2;
1104 
1105  const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
1106  const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
1107 
1108  G4Polyhedron *ph=new G4Polyhedron;
1109  typedef G4double G4double3[3];
1110  typedef G4int G4int4[4];
1111  G4double3* xyz = new G4double3[nnodes]; // number of nodes
1112  G4int4* faces = new G4int4[nfaces] ; // number of faces
1113  fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
1114  fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
1115  fInnerHype->GetFacets(k,n,xyz,faces,2) ;
1116  fFormerTwisted->GetFacets(k,n,xyz,faces,3) ;
1117  fOuterHype->GetFacets(k,n,xyz,faces,4) ;
1118  fLatterTwisted->GetFacets(k,n,xyz,faces,5) ;
1119 
1120  ph->createPolyhedron(nnodes,nfaces,xyz,faces);
1121 
1122  delete[] xyz;
1123  delete[] faces;
1124 
1125  return ph;
1126 }
1127 
1128 //=====================================================================
1129 //* CreateNUBS --------------------------------------------------------
1130 
1131 G4NURBS* G4TwistedTubs::CreateNURBS () const
1132 {
1133  G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
1134  G4double maxEndInnerRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
1135  return new G4NURBStube(maxEndInnerRad, maxEndOuterRad, fZHalfLength);
1136  // Tube for now!!!
1137 }
1138 
1139 //=====================================================================
1140 //* GetPolyhedron -----------------------------------------------------
1141 
1142 G4Polyhedron* G4TwistedTubs::GetPolyhedron () const
1143 {
1144  if ((!fpPolyhedron) ||
1145  (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1146  fpPolyhedron->GetNumberOfRotationSteps()))
1147  {
1148  delete fpPolyhedron;
1149  fpPolyhedron = CreatePolyhedron();
1150  }
1151  return fpPolyhedron;
1152 }
1153 
1154 //=====================================================================
1155 //* CreateSurfaces ----------------------------------------------------
1156 
1157 void G4TwistedTubs::CreateSurfaces()
1158 {
1159  // create 6 surfaces of TwistedTub
1160 
1161  G4ThreeVector x0(0, 0, fEndZ[0]);
1162  G4ThreeVector n (0, 0, -1);
1163 
1164  fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
1165  fEndInnerRadius, fEndOuterRadius,
1166  fDPhi, fEndPhi, fEndZ, -1) ;
1167 
1168  fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",
1169  fEndInnerRadius, fEndOuterRadius,
1170  fDPhi, fEndPhi, fEndZ, 1) ;
1171 
1172  G4RotationMatrix rotHalfDPhi;
1173  rotHalfDPhi.rotateZ(0.5*fDPhi);
1174 
1175  fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
1176  fEndInnerRadius, fEndOuterRadius,
1177  fDPhi, fEndPhi, fEndZ,
1178  fInnerRadius, fOuterRadius, fKappa,
1179  1 ) ;
1180  fFormerTwisted = new G4TwistTubsSide("FormerTwisted",
1181  fEndInnerRadius, fEndOuterRadius,
1182  fDPhi, fEndPhi, fEndZ,
1183  fInnerRadius, fOuterRadius, fKappa,
1184  -1 ) ;
1185 
1186  fInnerHype = new G4TwistTubsHypeSide("InnerHype",
1187  fEndInnerRadius, fEndOuterRadius,
1188  fDPhi, fEndPhi, fEndZ,
1189  fInnerRadius, fOuterRadius,fKappa,
1190  fTanInnerStereo, fTanOuterStereo, -1) ;
1191  fOuterHype = new G4TwistTubsHypeSide("OuterHype",
1192  fEndInnerRadius, fEndOuterRadius,
1193  fDPhi, fEndPhi, fEndZ,
1194  fInnerRadius, fOuterRadius,fKappa,
1195  fTanInnerStereo, fTanOuterStereo, 1) ;
1196 
1197 
1198  // set neighbour surfaces
1199  //
1200  fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
1201  fOuterHype, fFormerTwisted);
1202  fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
1203  fOuterHype, fFormerTwisted);
1204  fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
1205  fOuterHype, fUpperEndcap);
1206  fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
1207  fOuterHype, fUpperEndcap);
1208  fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
1209  fFormerTwisted, fUpperEndcap);
1210  fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
1211  fFormerTwisted, fUpperEndcap);
1212 }
1213 
1214 
1215 //=====================================================================
1216 //* GetEntityType -----------------------------------------------------
1217 
1218 G4GeometryType G4TwistedTubs::GetEntityType() const
1219 {
1220  return G4String("G4TwistedTubs");
1221 }
1222 
1223 //=====================================================================
1224 //* Clone -------------------------------------------------------------
1225 
1226 G4VSolid* G4TwistedTubs::Clone() const
1227 {
1228  return new G4TwistedTubs(*this);
1229 }
1230 
1231 //=====================================================================
1232 //* GetCubicVolume ----------------------------------------------------
1233 
1234 G4double G4TwistedTubs::GetCubicVolume()
1235 {
1236  if(fCubicVolume != 0.) {;}
1237  else { fCubicVolume = fDPhi*fZHalfLength*(fOuterRadius*fOuterRadius
1238  -fInnerRadius*fInnerRadius); }
1239  return fCubicVolume;
1240 }
1241 
1242 //=====================================================================
1243 //* GetSurfaceArea ----------------------------------------------------
1244 
1245 G4double G4TwistedTubs::GetSurfaceArea()
1246 {
1247  if(fSurfaceArea != 0.) {;}
1248  else { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
1249  return fSurfaceArea;
1250 }
1251 
1252 //=====================================================================
1253 //* GetPointOnSurface -------------------------------------------------
1254 
1255 G4ThreeVector G4TwistedTubs::GetPointOnSurface() const
1256 {
1257 
1258  G4double z = G4RandFlat::shoot(fEndZ[0],fEndZ[1]);
1259  G4double phi , phimin, phimax ;
1260  G4double x , xmin, xmax ;
1261  G4double r , rmin, rmax ;
1262 
1263  G4double a1 = fOuterHype->GetSurfaceArea() ;
1264  G4double a2 = fInnerHype->GetSurfaceArea() ;
1265  G4double a3 = fLatterTwisted->GetSurfaceArea() ;
1266  G4double a4 = fFormerTwisted->GetSurfaceArea() ;
1267  G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1268  G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1269 
1270  G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1271 
1272  if(chose < a1)
1273  {
1274 
1275  phimin = fOuterHype->GetBoundaryMin(z) ;
1276  phimax = fOuterHype->GetBoundaryMax(z) ;
1277  phi = G4RandFlat::shoot(phimin,phimax) ;
1278 
1279  return fOuterHype->SurfacePoint(phi,z,true) ;
1280 
1281  }
1282  else if ( (chose >= a1) && (chose < a1 + a2 ) )
1283  {
1284 
1285  phimin = fInnerHype->GetBoundaryMin(z) ;
1286  phimax = fInnerHype->GetBoundaryMax(z) ;
1287  phi = G4RandFlat::shoot(phimin,phimax) ;
1288 
1289  return fInnerHype->SurfacePoint(phi,z,true) ;
1290 
1291  }
1292  else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1293  {
1294 
1295  xmin = fLatterTwisted->GetBoundaryMin(z) ;
1296  xmax = fLatterTwisted->GetBoundaryMax(z) ;
1297  x = G4RandFlat::shoot(xmin,xmax) ;
1298 
1299  return fLatterTwisted->SurfacePoint(x,z,true) ;
1300 
1301  }
1302  else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1303  {
1304 
1305  xmin = fFormerTwisted->GetBoundaryMin(z) ;
1306  xmax = fFormerTwisted->GetBoundaryMax(z) ;
1307  x = G4RandFlat::shoot(xmin,xmax) ;
1308 
1309  return fFormerTwisted->SurfacePoint(x,z,true) ;
1310  }
1311  else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1312  {
1313  rmin = GetEndInnerRadius(0) ;
1314  rmax = GetEndOuterRadius(0) ;
1315  r = std::sqrt(G4RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
1316 
1317  phimin = fLowerEndcap->GetBoundaryMin(r) ;
1318  phimax = fLowerEndcap->GetBoundaryMax(r) ;
1319  phi = G4RandFlat::shoot(phimin,phimax) ;
1320 
1321  return fLowerEndcap->SurfacePoint(phi,r,true) ;
1322  }
1323  else
1324  {
1325  rmin = GetEndInnerRadius(1) ;
1326  rmax = GetEndOuterRadius(1) ;
1327  r = rmin + (rmax-rmin)*std::sqrt(G4RandFlat::shoot());
1328 
1329  phimin = fUpperEndcap->GetBoundaryMin(r) ;
1330  phimax = fUpperEndcap->GetBoundaryMax(r) ;
1331  phi = G4RandFlat::shoot(phimin,phimax) ;
1332 
1333  return fUpperEndcap->SurfacePoint(phi,r,true) ;
1334  }
1335 }
Generated on Sat May 25 2013 14:33:16 for Geant4 by   doxygen 1.8.4