G4TwistedTubs.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id: G4TwistedTubs.cc 67011 2013-01-29 16:17:41Z gcosmo $
00028 //
00029 // 
00030 // --------------------------------------------------------------------
00031 // GEANT 4 class source file
00032 //
00033 //
00034 // G4TwistTubsSide.cc
00035 //
00036 // Author: 
00037 //   01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
00038 //
00039 // History:
00040 //   13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
00041 //                 from original version in Jupiter-2.5.02 application.
00042 // --------------------------------------------------------------------
00043 
00044 #include "G4TwistedTubs.hh"
00045 
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4GeometryTolerance.hh"
00049 #include "G4VoxelLimits.hh"
00050 #include "G4AffineTransform.hh"
00051 #include "G4SolidExtentList.hh"
00052 #include "G4ClippablePolygon.hh"
00053 #include "G4VPVParameterisation.hh"
00054 #include "meshdefs.hh"
00055 
00056 #include "G4VGraphicsScene.hh"
00057 #include "G4Polyhedron.hh"
00058 #include "G4VisExtent.hh"
00059 #include "G4NURBS.hh"
00060 #include "G4NURBStube.hh"
00061 #include "G4NURBScylinder.hh"
00062 #include "G4NURBStubesector.hh"
00063 
00064 #include "Randomize.hh"
00065 
00066 //=====================================================================
00067 //* constructors ------------------------------------------------------
00068 
00069 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
00070                                    G4double  twistedangle,
00071                                    G4double  endinnerrad,
00072                                    G4double  endouterrad,
00073                                    G4double  halfzlen,
00074                                    G4double  dphi)
00075    : G4VSolid(pname), fDPhi(dphi), 
00076      fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
00077      fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
00078      fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00079 {
00080    if (endinnerrad < DBL_MIN)
00081    {
00082       G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00083                   FatalErrorInArgument, "Invalid end-inner-radius!");
00084    }
00085             
00086    G4double sinhalftwist = std::sin(0.5 * twistedangle);
00087 
00088    G4double endinnerradX = endinnerrad * sinhalftwist;
00089    G4double innerrad     = std::sqrt( endinnerrad * endinnerrad
00090                                  - endinnerradX * endinnerradX );
00091 
00092    G4double endouterradX = endouterrad * sinhalftwist;
00093    G4double outerrad     = std::sqrt( endouterrad * endouterrad
00094                                  - endouterradX * endouterradX );
00095    
00096    // temporary treatment!!
00097    SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
00098    CreateSurfaces();
00099 }
00100 
00101 G4TwistedTubs::G4TwistedTubs(const G4String &pname,     
00102                                    G4double  twistedangle,    
00103                                    G4double  endinnerrad,  
00104                                    G4double  endouterrad,  
00105                                    G4double  halfzlen,
00106                                    G4int     nseg,
00107                                    G4double  totphi)
00108    : G4VSolid(pname),
00109      fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
00110      fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
00111      fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00112 {
00113 
00114    if (!nseg)
00115    {
00116       std::ostringstream message;
00117       message << "Invalid number of segments." << G4endl
00118               << "        nseg = " << nseg;
00119       G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00120                   FatalErrorInArgument, message);
00121    }
00122    if (totphi == DBL_MIN || endinnerrad < DBL_MIN)
00123    {
00124       G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00125                 FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
00126    }
00127          
00128    G4double sinhalftwist = std::sin(0.5 * twistedangle);
00129 
00130    G4double endinnerradX = endinnerrad * sinhalftwist;
00131    G4double innerrad     = std::sqrt( endinnerrad * endinnerrad
00132                                  - endinnerradX * endinnerradX );
00133 
00134    G4double endouterradX = endouterrad * sinhalftwist;
00135    G4double outerrad     = std::sqrt( endouterrad * endouterrad
00136                                  - endouterradX * endouterradX );
00137    
00138    // temporary treatment!!
00139    fDPhi = totphi / nseg;
00140    SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
00141    CreateSurfaces();
00142 }
00143 
00144 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
00145                                    G4double  twistedangle,
00146                                    G4double  innerrad,
00147                                    G4double  outerrad,
00148                                    G4double  negativeEndz,
00149                                    G4double  positiveEndz,
00150                                    G4double  dphi)
00151    : G4VSolid(pname), fDPhi(dphi),
00152      fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
00153      fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
00154      fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00155 {
00156    if (innerrad < DBL_MIN)
00157    {
00158       G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00159                   FatalErrorInArgument, "Invalid end-inner-radius!");
00160    }
00161                  
00162    SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
00163    CreateSurfaces();
00164 }
00165 
00166 G4TwistedTubs::G4TwistedTubs(const G4String &pname,     
00167                                    G4double  twistedangle,    
00168                                    G4double  innerrad,  
00169                                    G4double  outerrad,  
00170                                    G4double  negativeEndz,
00171                                    G4double  positiveEndz,
00172                                    G4int     nseg,
00173                                    G4double  totphi)
00174    : G4VSolid(pname),
00175      fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
00176      fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
00177      fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00178 {
00179    if (!nseg)
00180    {
00181       std::ostringstream message;
00182       message << "Invalid number of segments." << G4endl
00183               << "        nseg = " << nseg;
00184       G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00185                   FatalErrorInArgument, message);
00186    }
00187    if (totphi == DBL_MIN || innerrad < DBL_MIN)
00188    {
00189       G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00190                 FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
00191    }
00192             
00193    fDPhi = totphi / nseg;
00194    SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
00195    CreateSurfaces();
00196 }
00197 
00198 //=====================================================================
00199 //* Fake default constructor ------------------------------------------
00200 
00201 G4TwistedTubs::G4TwistedTubs( __void__& a )
00202   : G4VSolid(a), fPhiTwist(0.), fInnerRadius(0.), fOuterRadius(0.), fDPhi(0.),
00203     fZHalfLength(0.), fInnerStereo(0.), fOuterStereo(0.), fTanInnerStereo(0.),
00204     fTanOuterStereo(0.), fKappa(0.), fInnerRadius2(0.), fOuterRadius2(0.),
00205     fTanInnerStereo2(0.), fTanOuterStereo2(0.), fLowerEndcap(0), fUpperEndcap(0),
00206     fLatterTwisted(0), fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
00207     fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00208 {
00209   fEndZ[0] = 0.; fEndZ[1] = 0.;
00210   fEndInnerRadius[0] = 0.; fEndInnerRadius[1] = 0.;
00211   fEndOuterRadius[0] = 0.; fEndOuterRadius[1] = 0.;
00212   fEndPhi[0] = 0.; fEndPhi[1] = 0.;
00213   fEndZ2[0] = 0.; fEndZ2[1] = 0.;
00214 }
00215 
00216 //=====================================================================
00217 //* destructor --------------------------------------------------------
00218 
00219 G4TwistedTubs::~G4TwistedTubs()
00220 {
00221    if (fLowerEndcap)   { delete fLowerEndcap;   }
00222    if (fUpperEndcap)   { delete fUpperEndcap;   }
00223    if (fLatterTwisted) { delete fLatterTwisted; }
00224    if (fFormerTwisted) { delete fFormerTwisted; }
00225    if (fInnerHype)     { delete fInnerHype;     }
00226    if (fOuterHype)     { delete fOuterHype;     }
00227    if (fpPolyhedron)   { delete fpPolyhedron;   }
00228 }
00229 
00230 //=====================================================================
00231 //* Copy constructor --------------------------------------------------
00232 
00233 G4TwistedTubs::G4TwistedTubs(const G4TwistedTubs& rhs)
00234   : G4VSolid(rhs), fPhiTwist(rhs.fPhiTwist),
00235     fInnerRadius(rhs.fInnerRadius), fOuterRadius(rhs.fOuterRadius),
00236     fDPhi(rhs.fDPhi), fZHalfLength(rhs.fZHalfLength),
00237     fInnerStereo(rhs.fInnerStereo), fOuterStereo(rhs.fOuterStereo),
00238     fTanInnerStereo(rhs.fTanInnerStereo), fTanOuterStereo(rhs.fTanOuterStereo),
00239     fKappa(rhs.fKappa), fInnerRadius2(rhs.fInnerRadius2), 
00240     fOuterRadius2(rhs.fOuterRadius2), fTanInnerStereo2(rhs.fTanInnerStereo2),
00241     fTanOuterStereo2(rhs.fTanOuterStereo2),
00242     fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0), fFormerTwisted(0),
00243     fInnerHype(0), fOuterHype(0),
00244     fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
00245     fpPolyhedron(0), fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
00246     fLastDistanceToIn(rhs.fLastDistanceToIn),
00247     fLastDistanceToOut(rhs.fLastDistanceToOut),
00248     fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
00249     fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
00250 {
00251   for (size_t i=0; i<2; ++i)
00252   {
00253     fEndZ[i] = rhs.fEndZ[i];
00254     fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
00255     fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
00256     fEndPhi[i] = rhs.fEndPhi[i];
00257     fEndZ2[i] = rhs.fEndZ2[i];
00258   }
00259   CreateSurfaces();
00260 }
00261 
00262 
00263 //=====================================================================
00264 //* Assignment operator -----------------------------------------------
00265 
00266 G4TwistedTubs& G4TwistedTubs::operator = (const G4TwistedTubs& rhs) 
00267 {
00268    // Check assignment to self
00269    //
00270    if (this == &rhs)  { return *this; }
00271 
00272    // Copy base class data
00273    //
00274    G4VSolid::operator=(rhs);
00275 
00276    // Copy data
00277    //
00278    fPhiTwist= rhs.fPhiTwist;
00279    fInnerRadius= rhs.fInnerRadius; fOuterRadius= rhs.fOuterRadius;
00280    fDPhi= rhs.fDPhi; fZHalfLength= rhs.fZHalfLength;
00281    fInnerStereo= rhs.fInnerStereo; fOuterStereo= rhs.fOuterStereo;
00282    fTanInnerStereo= rhs.fTanInnerStereo; fTanOuterStereo= rhs.fTanOuterStereo;
00283    fKappa= rhs.fKappa; fInnerRadius2= rhs.fInnerRadius2; 
00284    fOuterRadius2= rhs.fOuterRadius2; fTanInnerStereo2= rhs.fTanInnerStereo2;
00285    fTanOuterStereo2= rhs.fTanOuterStereo2;
00286    fLowerEndcap= fUpperEndcap= fLatterTwisted= fFormerTwisted= 0;
00287    fInnerHype= fOuterHype= 0;
00288    fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
00289    fpPolyhedron= 0;
00290    fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
00291    fLastDistanceToIn= rhs.fLastDistanceToIn;
00292    fLastDistanceToOut= rhs.fLastDistanceToOut;
00293    fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
00294    fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
00295  
00296    for (size_t i=0; i<2; ++i)
00297    {
00298      fEndZ[i] = rhs.fEndZ[i];
00299      fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
00300      fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
00301      fEndPhi[i] = rhs.fEndPhi[i];
00302      fEndZ2[i] = rhs.fEndZ2[i];
00303    }
00304  
00305    CreateSurfaces();
00306 
00307    return *this;
00308 }
00309 
00310 //=====================================================================
00311 //* ComputeDimensions -------------------------------------------------
00312 
00313 void G4TwistedTubs::ComputeDimensions(G4VPVParameterisation* /* p */ ,
00314                                       const G4int            /* n  */ ,
00315                                       const G4VPhysicalVolume* /* pRep */ )
00316 {
00317   G4Exception("G4TwistedTubs::ComputeDimensions()",
00318               "GeomSolids0001", FatalException,
00319               "G4TwistedTubs does not support Parameterisation.");
00320 }
00321 
00322 
00323 //=====================================================================
00324 //* CalculateExtent ---------------------------------------------------
00325 
00326 G4bool G4TwistedTubs::CalculateExtent( const EAxis              axis,
00327                                        const G4VoxelLimits     &voxelLimit,
00328                                        const G4AffineTransform &transform,
00329                                              G4double          &min,
00330                                              G4double          &max ) const
00331 {
00332 
00333   G4SolidExtentList  extentList( axis, voxelLimit );
00334   G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
00335                              fEndOuterRadius[0] : fEndOuterRadius[1]);
00336   G4double maxEndInnerRad = (fEndInnerRadius[0] > fEndInnerRadius[1] ?
00337                              fEndInnerRadius[0] : fEndInnerRadius[1]);
00338   G4double maxphi         = (std::fabs(fEndPhi[0]) > std::fabs(fEndPhi[1]) ?
00339                              std::fabs(fEndPhi[0]) : std::fabs(fEndPhi[1]));
00340    
00341   //
00342   // Choose phi size of our segment(s) based on constants as
00343   // defined in meshdefs.hh
00344   //
00345   // G4int numPhi = kMaxMeshSections;
00346   G4double sigPhi = 2*maxphi + fDPhi;
00347   G4double rFudge = 1.0/std::cos(0.5*sigPhi);
00348   G4double fudgeEndOuterRad = rFudge * maxEndOuterRad;
00349   
00350   //
00351   // We work around in phi building polygons along the way.
00352   // As a reasonable compromise between accuracy and
00353   // complexity (=cpu time), the following facets are chosen:
00354   //
00355   //   1. If fOuterRadius/maxEndOuterRad > 0.95, approximate
00356   //      the outer surface as a cylinder, and use one
00357   //      rectangular polygon (0-1) to build its mesh.
00358   //
00359   //      Otherwise, use two trapazoidal polygons that 
00360   //      meet at z = 0 (0-4-1)
00361   //
00362   //   2. If there is no inner surface, then use one
00363   //      polygon for each entire endcap.  (0) and (1)
00364   //
00365   //      Otherwise, use a trapazoidal polygon for each
00366   //      phi segment of each endcap.    (0-2) and (1-3)
00367   //
00368   //   3. For the inner surface, if fInnerRadius/maxEndInnerRad > 0.95,
00369   //      approximate the inner surface as a cylinder of
00370   //      radius fInnerRadius and use one rectangular polygon
00371   //      to build each phi segment of its mesh.   (2-3)
00372   //
00373   //      Otherwise, use one rectangular polygon centered
00374   //      at z = 0 (5-6) and two connecting trapazoidal polygons
00375   //      for each phi segment (2-5) and (3-6).
00376   //
00377 
00378   G4bool splitOuter = (fOuterRadius/maxEndOuterRad < 0.95);
00379   G4bool splitInner = (fInnerRadius/maxEndInnerRad < 0.95);
00380   
00381   //
00382   // Vertex assignments (v and w arrays)
00383   // [0] and [1] are mandatory
00384   // the rest are optional
00385   //
00386   //     +                     -
00387   //      [0]------[4]------[1]      <--- outer radius
00388   //       |                 |       
00389   //       |                 |       
00390   //      [2]---[5]---[6]---[3]      <--- inner radius
00391   //
00392 
00393   G4ClippablePolygon endPoly1, endPoly2;
00394   
00395   G4double phimax   = maxphi + 0.5*fDPhi;
00396   if ( phimax > pi/2)  { phimax = pi-phimax; }
00397   G4double phimin   = - phimax;
00398 
00399   G4ThreeVector v0, v1, v2, v3, v4, v5, v6;   // -ve phi verticies for polygon
00400   G4ThreeVector w0, w1, w2, w3, w4, w5, w6;   // +ve phi verticies for polygon
00401 
00402   //
00403   // decide verticies of -ve phi boundary
00404   //
00405   
00406   G4double cosPhi = std::cos(phimin);
00407   G4double sinPhi = std::sin(phimin);
00408  
00409   // Outer hyperbolic surface  
00410 
00411   v0 = transform.TransformPoint( 
00412        G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 
00413                      + fZHalfLength));
00414   v1 = transform.TransformPoint( 
00415        G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 
00416                      - fZHalfLength));
00417   if (splitOuter)
00418   {
00419      v4 = transform.TransformPoint(
00420           G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 0));
00421   }
00422   
00423   // Inner hyperbolic surface  
00424 
00425   G4double zInnerSplit = 0.;
00426   if (splitInner)
00427   {
00428      v2 = transform.TransformPoint( 
00429           G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, 
00430                         + fZHalfLength));
00431      v3 = transform.TransformPoint( 
00432           G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, 
00433                         - fZHalfLength));
00434       
00435      // Find intersection of tangential line of inner
00436      // surface at z = fZHalfLength and line r=fInnerRadius.
00437      G4double dr = fZHalfLength * fTanInnerStereo2;
00438      G4double dz = maxEndInnerRad;
00439      zInnerSplit = fZHalfLength + (fInnerRadius - maxEndInnerRad) * dz / dr;
00440 
00441      // Build associated vertices
00442      v5 = transform.TransformPoint( 
00443           G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, 
00444                         + zInnerSplit));
00445      v6 = transform.TransformPoint( 
00446           G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, 
00447                         - zInnerSplit));
00448   }
00449   else
00450   {
00451      v2 = transform.TransformPoint(
00452           G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, 
00453                         + fZHalfLength));
00454      v3 = transform.TransformPoint(
00455           G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, 
00456                         - fZHalfLength));
00457   }
00458 
00459   //
00460   // decide vertices of +ve phi boundary
00461   // 
00462 
00463   cosPhi = std::cos(phimax);
00464   sinPhi = std::sin(phimax);
00465    
00466   // Outer hyperbolic surface  
00467   
00468   w0 = transform.TransformPoint(
00469        G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
00470                      + fZHalfLength));
00471   w1 = transform.TransformPoint(
00472        G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
00473                      - fZHalfLength));
00474   if (splitOuter)
00475   {
00476      G4double r = rFudge*fOuterRadius;
00477      
00478      w4 = transform.TransformPoint(G4ThreeVector( r*cosPhi, r*sinPhi, 0 ));
00479       
00480      AddPolyToExtent( v0, v4, w4, w0, voxelLimit, axis, extentList );
00481      AddPolyToExtent( v4, v1, w1, w4, voxelLimit, axis, extentList );
00482   }
00483   else
00484   {
00485      AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList );
00486   }
00487   
00488   // Inner hyperbolic surface
00489   
00490   if (splitInner)
00491   {
00492      w2 = transform.TransformPoint(
00493           G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, 
00494                         + fZHalfLength));
00495      w3 = transform.TransformPoint(
00496           G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, 
00497                         - fZHalfLength));
00498           
00499      w5 = transform.TransformPoint(
00500           G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00501                         + zInnerSplit));
00502      w6 = transform.TransformPoint(
00503           G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00504                         - zInnerSplit));
00505                                    
00506      AddPolyToExtent( v3, v6, w6, w3, voxelLimit, axis, extentList );
00507      AddPolyToExtent( v6, v5, w5, w6, voxelLimit, axis, extentList );
00508      AddPolyToExtent( v5, v2, w2, w5, voxelLimit, axis, extentList );
00509      
00510   }
00511   else
00512   {
00513      w2 = transform.TransformPoint(
00514           G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00515                         + fZHalfLength));
00516      w3 = transform.TransformPoint(
00517           G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00518                         - fZHalfLength));
00519 
00520      AddPolyToExtent( v3, v2, w2, w3, voxelLimit, axis, extentList );
00521   }
00522 
00523   //
00524   // Endplate segments
00525   //
00526   AddPolyToExtent( v1, v3, w3, w1, voxelLimit, axis, extentList );
00527   AddPolyToExtent( v2, v0, w0, w2, voxelLimit, axis, extentList );
00528   
00529   //
00530   // Return min/max value
00531   //
00532   return extentList.GetExtent( min, max );
00533 }
00534 
00535 
00536 //=====================================================================
00537 //* AddPolyToExtent ---------------------------------------------------
00538 
00539 void G4TwistedTubs::AddPolyToExtent( const G4ThreeVector &v0,
00540                                      const G4ThreeVector &v1,
00541                                      const G4ThreeVector &w1,
00542                                      const G4ThreeVector &w0,
00543                                      const G4VoxelLimits &voxelLimit,
00544                                      const EAxis          axis,
00545                                      G4SolidExtentList   &extentList ) 
00546 {
00547     // Utility function for CalculateExtent
00548     //
00549     G4ClippablePolygon phiPoly;
00550 
00551     phiPoly.AddVertexInOrder( v0 );
00552     phiPoly.AddVertexInOrder( v1 );
00553     phiPoly.AddVertexInOrder( w1 );
00554     phiPoly.AddVertexInOrder( w0 );
00555 
00556     if (phiPoly.PartialClip( voxelLimit, axis ))
00557     {
00558         phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
00559         extentList.AddSurface( phiPoly );
00560     }
00561 }
00562 
00563 
00564 //=====================================================================
00565 //* Inside ------------------------------------------------------------
00566 
00567 EInside G4TwistedTubs::Inside(const G4ThreeVector& p) const
00568 {
00569 
00570    static const G4double halftol
00571      = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00572    // static G4int timerid = -1;
00573    // G4Timer timer(timerid, "G4TwistedTubs", "Inside");
00574    // timer.Start();
00575    
00576 
00577 
00578    G4ThreeVector *tmpp;
00579    EInside       *tmpinside;
00580    if (fLastInside.p == p)
00581    {
00582      return fLastInside.inside;
00583    }
00584    else
00585    {
00586       tmpp      = const_cast<G4ThreeVector*>(&(fLastInside.p));
00587       tmpinside = const_cast<EInside*>(&(fLastInside.inside));
00588       tmpp->set(p.x(), p.y(), p.z());
00589    }
00590    
00591    EInside  outerhypearea = ((G4TwistTubsHypeSide *)fOuterHype)->Inside(p);
00592    G4double innerhyperho  = ((G4TwistTubsHypeSide *)fInnerHype)->GetRhoAtPZ(p);
00593    G4double distanceToOut = p.getRho() - innerhyperho; // +ve: inside
00594 
00595    if ((outerhypearea == kOutside) || (distanceToOut < -halftol))
00596    {
00597       *tmpinside = kOutside;
00598    }
00599    else if (outerhypearea == kSurface)
00600    {
00601       *tmpinside = kSurface;
00602    }
00603    else
00604    {
00605       if (distanceToOut <= halftol)
00606       {
00607          *tmpinside = kSurface;
00608       }
00609       else
00610       {
00611          *tmpinside = kInside;
00612       }
00613    }
00614 
00615    return fLastInside.inside;
00616 }
00617 
00618 //=====================================================================
00619 //* SurfaceNormal -----------------------------------------------------
00620 
00621 G4ThreeVector G4TwistedTubs::SurfaceNormal(const G4ThreeVector& p) const
00622 {
00623    //
00624    // return the normal unit vector to the Hyperbolical Surface at a point 
00625    // p on (or nearly on) the surface
00626    //
00627    // Which of the three or four surfaces are we closest to?
00628    //
00629 
00630    if (fLastNormal.p == p)
00631    {
00632       return fLastNormal.vec;
00633    }    
00634    G4ThreeVector *tmpp          =
00635      const_cast<G4ThreeVector*>(&(fLastNormal.p));
00636    G4ThreeVector *tmpnormal     =
00637      const_cast<G4ThreeVector*>(&(fLastNormal.vec));
00638    G4VTwistSurface **tmpsurface =
00639      const_cast<G4VTwistSurface**>(fLastNormal.surface);
00640    tmpp->set(p.x(), p.y(), p.z());
00641 
00642    G4double      distance = kInfinity;
00643 
00644    G4VTwistSurface *surfaces[6];
00645    surfaces[0] = fLatterTwisted;
00646    surfaces[1] = fFormerTwisted;
00647    surfaces[2] = fInnerHype;
00648    surfaces[3] = fOuterHype;
00649    surfaces[4] = fLowerEndcap;
00650    surfaces[5] = fUpperEndcap;
00651 
00652    G4ThreeVector xx;
00653    G4ThreeVector bestxx;
00654    G4int i;
00655    G4int besti = -1;
00656    for (i=0; i< 6; i++)
00657    {
00658       G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
00659       if (tmpdistance < distance)
00660       {
00661          distance = tmpdistance;
00662          bestxx = xx;
00663          besti = i; 
00664       }
00665    }
00666 
00667    tmpsurface[0] = surfaces[besti];
00668    *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
00669    
00670    return fLastNormal.vec;
00671 }
00672 
00673 //=====================================================================
00674 //* DistanceToIn (p, v) -----------------------------------------------
00675 
00676 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p,
00677                                       const G4ThreeVector& v ) const
00678 {
00679 
00680    // DistanceToIn (p, v):
00681    // Calculate distance to surface of shape from `outside' 
00682    // along with the v, allowing for tolerance.
00683    // The function returns kInfinity if no intersection or
00684    // just grazing within tolerance.
00685 
00686    //
00687    // checking last value
00688    //
00689    
00690    G4ThreeVector *tmpp;
00691    G4ThreeVector *tmpv;
00692    G4double      *tmpdist;
00693    if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v))
00694    {
00695      return fLastDistanceToIn.value;
00696    }
00697    else
00698    {
00699       tmpp    = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
00700       tmpv    = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
00701       tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
00702       tmpp->set(p.x(), p.y(), p.z());
00703       tmpv->set(v.x(), v.y(), v.z());
00704    }
00705 
00706    //
00707    // Calculate DistanceToIn(p,v)
00708    //
00709    
00710    EInside currentside = Inside(p);
00711 
00712    if (currentside == kInside)
00713    {
00714    }
00715    else
00716    {
00717      if (currentside == kSurface)
00718      {
00719        // particle is just on a boundary.
00720        // If the particle is entering to the volume, return 0.
00721        //
00722        G4ThreeVector normal = SurfaceNormal(p);
00723        if (normal*v < 0)
00724        {
00725          *tmpdist = 0;
00726          return fLastDistanceToInWithV.value;
00727        } 
00728      }
00729    }
00730       
00731    // now, we can take smallest positive distance.
00732    
00733    // Initialize
00734    //
00735    G4double      distance = kInfinity;   
00736 
00737    // find intersections and choose nearest one.
00738    //
00739    G4VTwistSurface *surfaces[6];
00740    surfaces[0] = fLowerEndcap;
00741    surfaces[1] = fUpperEndcap;
00742    surfaces[2] = fLatterTwisted;
00743    surfaces[3] = fFormerTwisted;
00744    surfaces[4] = fInnerHype;
00745    surfaces[5] = fOuterHype;
00746    
00747    G4ThreeVector xx;
00748    G4ThreeVector bestxx;
00749    G4int i;
00750    for (i=0; i< 6; i++)
00751    {
00752       G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
00753       if (tmpdistance < distance)
00754       {
00755          distance = tmpdistance;
00756          bestxx = xx;
00757       }
00758    }
00759    *tmpdist = distance;
00760 
00761    return fLastDistanceToInWithV.value;
00762 }
00763  
00764 //=====================================================================
00765 //* DistanceToIn (p) --------------------------------------------------
00766 
00767 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p) const
00768 {
00769    // DistanceToIn(p):
00770    // Calculate distance to surface of shape from `outside',
00771    // allowing for tolerance
00772    
00773    //
00774    // checking last value
00775    //
00776    
00777    G4ThreeVector *tmpp;
00778    G4double      *tmpdist;
00779    if (fLastDistanceToIn.p == p)
00780    {
00781      return fLastDistanceToIn.value;
00782    }
00783    else
00784    {
00785       tmpp    = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
00786       tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
00787       tmpp->set(p.x(), p.y(), p.z());
00788    }
00789 
00790    //
00791    // Calculate DistanceToIn(p) 
00792    //
00793    
00794    EInside currentside = Inside(p);
00795 
00796    switch (currentside)
00797    {
00798       case (kInside) :
00799       {}
00800       case (kSurface) :
00801       {
00802          *tmpdist = 0.;
00803          return fLastDistanceToIn.value;
00804       }
00805       case (kOutside) :
00806       {
00807          // Initialize
00808          G4double      distance = kInfinity;   
00809 
00810          // find intersections and choose nearest one.
00811          G4VTwistSurface *surfaces[6];
00812          surfaces[0] = fLowerEndcap;
00813          surfaces[1] = fUpperEndcap;
00814          surfaces[2] = fLatterTwisted;
00815          surfaces[3] = fFormerTwisted;
00816          surfaces[4] = fInnerHype;
00817          surfaces[5] = fOuterHype;
00818 
00819          G4int i;
00820          G4ThreeVector xx;
00821          G4ThreeVector bestxx;
00822          for (i=0; i< 6; i++)
00823          {
00824             G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
00825             if (tmpdistance < distance)
00826             {
00827                distance = tmpdistance;
00828                bestxx = xx;
00829             }
00830          }
00831          *tmpdist = distance;
00832          return fLastDistanceToIn.value;
00833       }
00834       default :
00835       {
00836          G4Exception("G4TwistedTubs::DistanceToIn(p)", "GeomSolids0003",
00837                      FatalException, "Unknown point location!");
00838       }
00839    } // switch end
00840 
00841    return kInfinity;
00842 }
00843 
00844 //=====================================================================
00845 //* DistanceToOut (p, v) ----------------------------------------------
00846 
00847 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p,
00848                                        const G4ThreeVector& v,
00849                                        const G4bool calcNorm,
00850                                        G4bool *validNorm,
00851                                        G4ThreeVector *norm ) const
00852 {
00853    // DistanceToOut (p, v):
00854    // Calculate distance to surface of shape from `inside'
00855    // along with the v, allowing for tolerance.
00856    // The function returns kInfinity if no intersection or
00857    // just grazing within tolerance.
00858 
00859    //
00860    // checking last value
00861    //
00862    
00863    G4ThreeVector *tmpp;
00864    G4ThreeVector *tmpv;
00865    G4double      *tmpdist;
00866    if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) )
00867    {
00868      return fLastDistanceToOutWithV.value;
00869    }
00870    else
00871    {
00872       tmpp    = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
00873       tmpv    = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
00874       tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
00875       tmpp->set(p.x(), p.y(), p.z());
00876       tmpv->set(v.x(), v.y(), v.z());
00877    }
00878 
00879    //
00880    // Calculate DistanceToOut(p,v)
00881    //
00882    
00883    EInside currentside = Inside(p);
00884 
00885    if (currentside == kOutside)
00886    {
00887    }
00888    else
00889    {
00890      if (currentside == kSurface)
00891      {
00892        // particle is just on a boundary.
00893        // If the particle is exiting from the volume, return 0.
00894        //
00895        G4ThreeVector normal = SurfaceNormal(p);
00896        G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
00897        if (normal*v > 0)
00898        {
00899          if (calcNorm)
00900          {
00901            *norm = (blockedsurface->GetNormal(p, true));
00902            *validNorm = blockedsurface->IsValidNorm();
00903          }
00904          *tmpdist = 0.;
00905          return fLastDistanceToOutWithV.value;
00906        }
00907      }
00908    }
00909       
00910    // now, we can take smallest positive distance.
00911    
00912    // Initialize
00913    //
00914    G4double      distance = kInfinity;
00915        
00916    // find intersections and choose nearest one.
00917    //
00918    G4VTwistSurface *surfaces[6];
00919    surfaces[0] = fLatterTwisted;
00920    surfaces[1] = fFormerTwisted;
00921    surfaces[2] = fInnerHype;
00922    surfaces[3] = fOuterHype;
00923    surfaces[4] = fLowerEndcap;
00924    surfaces[5] = fUpperEndcap;
00925    
00926    G4int i;
00927    G4int besti = -1;
00928    G4ThreeVector xx;
00929    G4ThreeVector bestxx;
00930    for (i=0; i< 6; i++)
00931    {
00932       G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
00933       if (tmpdistance < distance)
00934       {
00935          distance = tmpdistance;
00936          bestxx = xx; 
00937          besti = i;
00938       }
00939    }
00940 
00941    if (calcNorm)
00942    {
00943       if (besti != -1)
00944       {
00945          *norm = (surfaces[besti]->GetNormal(p, true));
00946          *validNorm = surfaces[besti]->IsValidNorm();
00947       }
00948    }
00949 
00950    *tmpdist = distance;
00951 
00952    return fLastDistanceToOutWithV.value;
00953 }
00954 
00955 
00956 //=====================================================================
00957 //* DistanceToOut (p) ----------------------------------------------
00958 
00959 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p ) const
00960 {
00961    // DistanceToOut(p):
00962    // Calculate distance to surface of shape from `inside', 
00963    // allowing for tolerance
00964 
00965    //
00966    // checking last value
00967    //
00968    
00969    G4ThreeVector *tmpp;
00970    G4double      *tmpdist;
00971    if (fLastDistanceToOut.p == p)
00972    {
00973       return fLastDistanceToOut.value;
00974    }
00975    else
00976    {
00977       tmpp    = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
00978       tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
00979       tmpp->set(p.x(), p.y(), p.z());
00980    }
00981    
00982    //
00983    // Calculate DistanceToOut(p)
00984    //
00985    
00986    EInside currentside = Inside(p);
00987 
00988    switch (currentside)
00989    {
00990       case (kOutside) :
00991       {
00992       }
00993       case (kSurface) :
00994       {
00995         *tmpdist = 0.;
00996          return fLastDistanceToOut.value;
00997       }
00998       case (kInside) :
00999       {
01000          // Initialize
01001          G4double      distance = kInfinity;
01002    
01003          // find intersections and choose nearest one.
01004          G4VTwistSurface *surfaces[6];
01005          surfaces[0] = fLatterTwisted;
01006          surfaces[1] = fFormerTwisted;
01007          surfaces[2] = fInnerHype;
01008          surfaces[3] = fOuterHype;
01009          surfaces[4] = fLowerEndcap;
01010          surfaces[5] = fUpperEndcap;
01011 
01012          G4int i;
01013          G4ThreeVector xx;
01014          G4ThreeVector bestxx;
01015          for (i=0; i< 6; i++)
01016          {
01017             G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
01018             if (tmpdistance < distance)
01019             {
01020                distance = tmpdistance;
01021                bestxx = xx;
01022             }
01023          }
01024          *tmpdist = distance;
01025    
01026          return fLastDistanceToOut.value;
01027       }
01028       default :
01029       {
01030          G4Exception("G4TwistedTubs::DistanceToOut(p)", "GeomSolids0003",
01031                      FatalException, "Unknown point location!");
01032       }
01033    } // switch end
01034 
01035    return 0;
01036 }
01037 
01038 //=====================================================================
01039 //* StreamInfo --------------------------------------------------------
01040 
01041 std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
01042 {
01043   //
01044   // Stream object contents to an output stream
01045   //
01046   G4int oldprc = os.precision(16);
01047   os << "-----------------------------------------------------------\n"
01048      << "    *** Dump for solid - " << GetName() << " ***\n"
01049      << "    ===================================================\n"
01050      << " Solid type: G4TwistedTubs\n"
01051      << " Parameters: \n"
01052      << "    -ve end Z              : " << fEndZ[0]/mm << " mm \n"
01053      << "    +ve end Z              : " << fEndZ[1]/mm << " mm \n"
01054      << "    inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
01055      << "    inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
01056      << "    outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
01057      << "    outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
01058      << "    inner radius (z=0)     : " << fInnerRadius/mm << " mm \n"
01059      << "    outer radius (z=0)     : " << fOuterRadius/mm << " mm \n"
01060      << "    twisted angle          : " << fPhiTwist/degree << " degrees \n"
01061      << "    inner stereo angle     : " << fInnerStereo/degree << " degrees \n"
01062      << "    outer stereo angle     : " << fOuterStereo/degree << " degrees \n"
01063      << "    phi-width of a piece   : " << fDPhi/degree << " degrees \n"
01064      << "-----------------------------------------------------------\n";
01065   os.precision(oldprc);
01066 
01067   return os;
01068 }
01069 
01070 
01071 //=====================================================================
01072 //* DiscribeYourselfTo ------------------------------------------------
01073 
01074 void G4TwistedTubs::DescribeYourselfTo (G4VGraphicsScene& scene) const 
01075 {
01076   scene.AddSolid (*this);
01077 }
01078 
01079 //=====================================================================
01080 //* GetExtent ---------------------------------------------------------
01081 
01082 G4VisExtent G4TwistedTubs::GetExtent() const 
01083 {
01084   // Define the sides of the box into which the G4Tubs instance would fit.
01085 
01086   G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
01087   return G4VisExtent( -maxEndOuterRad, maxEndOuterRad, 
01088                       -maxEndOuterRad, maxEndOuterRad, 
01089                       -fZHalfLength, fZHalfLength );
01090 }
01091 
01092 //=====================================================================
01093 //* CreatePolyhedron --------------------------------------------------
01094 
01095 G4Polyhedron* G4TwistedTubs::CreatePolyhedron () const 
01096 {
01097   // number of meshes
01098   //
01099   G4double dA = std::max(fDPhi,fPhiTwist);
01100   const G4int k =
01101     G4int(G4Polyhedron::GetNumberOfRotationSteps() * dA / twopi) + 2;
01102   const G4int n =
01103     G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2;
01104 
01105   const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
01106   const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
01107 
01108   G4Polyhedron *ph=new G4Polyhedron;
01109   typedef G4double G4double3[3];
01110   typedef G4int G4int4[4];
01111   G4double3* xyz = new G4double3[nnodes];  // number of nodes 
01112   G4int4*  faces = new G4int4[nfaces] ;    // number of faces
01113   fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
01114   fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
01115   fInnerHype->GetFacets(k,n,xyz,faces,2) ;
01116   fFormerTwisted->GetFacets(k,n,xyz,faces,3) ;
01117   fOuterHype->GetFacets(k,n,xyz,faces,4) ;
01118   fLatterTwisted->GetFacets(k,n,xyz,faces,5) ;
01119 
01120   ph->createPolyhedron(nnodes,nfaces,xyz,faces);
01121 
01122   delete[] xyz;
01123   delete[] faces;
01124 
01125   return ph;
01126 }
01127 
01128 //=====================================================================
01129 //* CreateNUBS --------------------------------------------------------
01130 
01131 G4NURBS* G4TwistedTubs::CreateNURBS () const 
01132 {
01133    G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
01134    G4double maxEndInnerRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
01135    return new G4NURBStube(maxEndInnerRad, maxEndOuterRad, fZHalfLength); 
01136    // Tube for now!!!
01137 }
01138 
01139 //=====================================================================
01140 //* GetPolyhedron -----------------------------------------------------
01141 
01142 G4Polyhedron* G4TwistedTubs::GetPolyhedron () const
01143 {
01144   if ((!fpPolyhedron) ||
01145       (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
01146        fpPolyhedron->GetNumberOfRotationSteps()))
01147   {
01148     delete fpPolyhedron;
01149     fpPolyhedron = CreatePolyhedron();
01150   }
01151   return fpPolyhedron;
01152 }
01153 
01154 //=====================================================================
01155 //* CreateSurfaces ----------------------------------------------------
01156 
01157 void G4TwistedTubs::CreateSurfaces() 
01158 {
01159    // create 6 surfaces of TwistedTub
01160 
01161    G4ThreeVector x0(0, 0, fEndZ[0]);
01162    G4ThreeVector n (0, 0, -1);
01163 
01164    fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
01165                                     fEndInnerRadius, fEndOuterRadius,
01166                                     fDPhi, fEndPhi, fEndZ, -1) ;
01167 
01168    fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",  
01169                                     fEndInnerRadius, fEndOuterRadius,
01170                                     fDPhi, fEndPhi, fEndZ, 1) ;
01171 
01172    G4RotationMatrix    rotHalfDPhi;
01173    rotHalfDPhi.rotateZ(0.5*fDPhi);
01174 
01175    fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
01176                                          fEndInnerRadius, fEndOuterRadius,
01177                                          fDPhi, fEndPhi, fEndZ, 
01178                                          fInnerRadius, fOuterRadius, fKappa,
01179                                          1 ) ; 
01180    fFormerTwisted = new G4TwistTubsSide("FormerTwisted", 
01181                                          fEndInnerRadius, fEndOuterRadius,
01182                                          fDPhi, fEndPhi, fEndZ, 
01183                                          fInnerRadius, fOuterRadius, fKappa,
01184                                          -1 ) ; 
01185 
01186    fInnerHype = new G4TwistTubsHypeSide("InnerHype",
01187                                         fEndInnerRadius, fEndOuterRadius,
01188                                         fDPhi, fEndPhi, fEndZ, 
01189                                         fInnerRadius, fOuterRadius,fKappa,
01190                                         fTanInnerStereo, fTanOuterStereo, -1) ;
01191    fOuterHype = new G4TwistTubsHypeSide("OuterHype", 
01192                                         fEndInnerRadius, fEndOuterRadius,
01193                                         fDPhi, fEndPhi, fEndZ, 
01194                                         fInnerRadius, fOuterRadius,fKappa,
01195                                         fTanInnerStereo, fTanOuterStereo, 1) ;
01196 
01197 
01198    // set neighbour surfaces
01199    //
01200    fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
01201                                fOuterHype, fFormerTwisted);
01202    fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
01203                                fOuterHype, fFormerTwisted);
01204    fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
01205                                  fOuterHype, fUpperEndcap);
01206    fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
01207                                  fOuterHype, fUpperEndcap);
01208    fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
01209                              fFormerTwisted, fUpperEndcap);
01210    fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
01211                              fFormerTwisted, fUpperEndcap);
01212 }
01213 
01214 
01215 //=====================================================================
01216 //* GetEntityType -----------------------------------------------------
01217 
01218 G4GeometryType G4TwistedTubs::GetEntityType() const
01219 {
01220   return G4String("G4TwistedTubs");
01221 }
01222 
01223 //=====================================================================
01224 //* Clone -------------------------------------------------------------
01225 
01226 G4VSolid* G4TwistedTubs::Clone() const
01227 {
01228   return new G4TwistedTubs(*this);
01229 }
01230 
01231 //=====================================================================
01232 //* GetCubicVolume ----------------------------------------------------
01233 
01234 G4double G4TwistedTubs::GetCubicVolume()
01235 {
01236   if(fCubicVolume != 0.) {;}
01237   else   { fCubicVolume = fDPhi*fZHalfLength*(fOuterRadius*fOuterRadius
01238                                              -fInnerRadius*fInnerRadius); }
01239   return fCubicVolume;
01240 }
01241 
01242 //=====================================================================
01243 //* GetSurfaceArea ----------------------------------------------------
01244 
01245 G4double G4TwistedTubs::GetSurfaceArea()
01246 {
01247   if(fSurfaceArea != 0.) {;}
01248   else   { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
01249   return fSurfaceArea;
01250 }
01251 
01252 //=====================================================================
01253 //* GetPointOnSurface -------------------------------------------------
01254 
01255 G4ThreeVector G4TwistedTubs::GetPointOnSurface() const
01256 {
01257 
01258   G4double  z = G4RandFlat::shoot(fEndZ[0],fEndZ[1]);
01259   G4double phi , phimin, phimax ;
01260   G4double x   , xmin,   xmax ;
01261   G4double r   , rmin,   rmax ;
01262 
01263   G4double a1 = fOuterHype->GetSurfaceArea() ;
01264   G4double a2 = fInnerHype->GetSurfaceArea() ;
01265   G4double a3 = fLatterTwisted->GetSurfaceArea() ;
01266   G4double a4 = fFormerTwisted->GetSurfaceArea() ;
01267   G4double a5 = fLowerEndcap->GetSurfaceArea()  ;
01268   G4double a6 = fUpperEndcap->GetSurfaceArea() ;
01269 
01270   G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
01271 
01272   if(chose < a1)
01273   {
01274 
01275     phimin = fOuterHype->GetBoundaryMin(z) ;
01276     phimax = fOuterHype->GetBoundaryMax(z) ;
01277     phi = G4RandFlat::shoot(phimin,phimax) ;
01278 
01279     return fOuterHype->SurfacePoint(phi,z,true) ;
01280 
01281   }
01282   else if ( (chose >= a1) && (chose < a1 + a2 ) )
01283   {
01284 
01285     phimin = fInnerHype->GetBoundaryMin(z) ;
01286     phimax = fInnerHype->GetBoundaryMax(z) ;
01287     phi = G4RandFlat::shoot(phimin,phimax) ;
01288 
01289     return fInnerHype->SurfacePoint(phi,z,true) ;
01290 
01291   }
01292   else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) ) 
01293   {
01294 
01295     xmin = fLatterTwisted->GetBoundaryMin(z) ; 
01296     xmax = fLatterTwisted->GetBoundaryMax(z) ;
01297     x = G4RandFlat::shoot(xmin,xmax) ;
01298     
01299     return fLatterTwisted->SurfacePoint(x,z,true) ;
01300 
01301   }
01302   else if ( (chose >= a1 + a2 + a3  ) && (chose < a1 + a2 + a3 + a4  ) )
01303   {
01304 
01305     xmin = fFormerTwisted->GetBoundaryMin(z) ; 
01306     xmax = fFormerTwisted->GetBoundaryMax(z) ;
01307     x = G4RandFlat::shoot(xmin,xmax) ;
01308 
01309     return fFormerTwisted->SurfacePoint(x,z,true) ;
01310    }
01311   else if( (chose >= a1 + a2 + a3 + a4  )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
01312   {
01313     rmin = GetEndInnerRadius(0) ;
01314     rmax = GetEndOuterRadius(0) ;
01315     r = std::sqrt(G4RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
01316 
01317     phimin = fLowerEndcap->GetBoundaryMin(r) ; 
01318     phimax = fLowerEndcap->GetBoundaryMax(r) ;
01319     phi    = G4RandFlat::shoot(phimin,phimax) ;
01320 
01321     return fLowerEndcap->SurfacePoint(phi,r,true) ;
01322   }
01323   else
01324   {
01325     rmin = GetEndInnerRadius(1) ;
01326     rmax = GetEndOuterRadius(1) ;
01327     r = rmin + (rmax-rmin)*std::sqrt(G4RandFlat::shoot());
01328 
01329     phimin = fUpperEndcap->GetBoundaryMin(r) ; 
01330     phimax = fUpperEndcap->GetBoundaryMax(r) ;
01331     phi    = G4RandFlat::shoot(phimin,phimax) ;
01332 
01333     return fUpperEndcap->SurfacePoint(phi,r,true) ;
01334   }
01335 }

Generated on Mon May 27 17:50:04 2013 for Geant4 by  doxygen 1.4.7