include/Geant4/G4VTwistSurface.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4VTwistSurface class inline methods
0027 //
0028 // Author: Kotoyo Hoshina (Chiba University), 01.08.2002 - Created.
0029 //         Oliver Link (CERN), 13.11.2003 - Integration in Geant4
0030 //               from original version in Jupiter-2.5.02 application.
0031 // --------------------------------------------------------------------
0032
0033 struct Intersection
0034 {
0035   G4double phi ;  // parameter phi
0036   G4double u ;    // parameter u
0037   G4ThreeVector xx ;   // intersection point in cartesian
0038   G4double distance ;  // distance to intersection
0039   G4int areacode ;     // the areacode of the intersection
0040   G4bool isvalid ;     // valid intersection ??
0041
0042 };
0043
0044 inline
0045 G4bool DistanceSort( const Intersection& a, const Intersection& b)
0046 {
0047   return a.distance < b.distance ;
0048 }
0049
0050 inline
0051 G4bool EqualIntersection( const Intersection& a, const Intersection& b)
0052 {
0053   return ( ( a.xx - b.xx ).mag() < 1E-9*CLHEP::mm ) ;
0054 }
0055
0056 //=====================================================================
0057 //* DistanceToPlaneWithV ----------------------------------------------
0058
0059 inline
0060 G4double G4VTwistSurface::DistanceToPlaneWithV(const G4ThreeVector& p,
0061                                                const G4ThreeVector& v,
0062                                                const G4ThreeVector& x0,
0063                                                const G4ThreeVector& n0,
0064                                                      G4ThreeVector& xx)
0065 {
0066    G4double q = n0 * v;
0067    G4double t = kInfinity;
0068    if (q != 0.0) { t = (n0 * (x0 - p)) / q; }
0069    xx = p + t * v;
0070    return t;
0071 }
0072
0073 //=====================================================================
0074 //* DistanceToPlane ---------------------------------------------------
0075
0076 inline
0077 G4double G4VTwistSurface::DistanceToPlane(const G4ThreeVector& p,
0078                                      const G4ThreeVector& x0,
0079                                      const G4ThreeVector& n0,
0080                                            G4ThreeVector& xx)
0081 {
0082    // DistanceToPlane :
0083    // Calculate distance to plane in local coordinate,
0084    // then return distance and global intersection points.
0085    //
0086    // p          - location of flying particle
0087    // x0         - reference point of surface
0088    // xx         - a foot of perpendicular line from p to the plane
0089    // t          - distance from xx to p
0090    // n          - a unit normal of this plane from plane to p.
0091    //
0092    // equation of plane:
0093    //      n*(x - x0) = 0;
0094    //
0095    // vector to xx:
0096    //      xx = p - t*n
0097    //
0098    //         where
0099    //         t = n * (p - x0) / std::fabs(n)
0100    //
0101    G4double t;
0102    G4ThreeVector n = n0.unit();
0103    t = n * (p - x0);
0104    xx = p - t * n;
0105    return t;
0106 }
0107
0108 //=====================================================================
0109 //* DistanceToPlane ---------------------------------------------------
0110
0111 inline
0112 G4double G4VTwistSurface::DistanceToPlane(const G4ThreeVector& p,
0113                                           const G4ThreeVector& x0,
0114                                           const G4ThreeVector& t1,
0115                                           const G4ThreeVector& t2,
0116                                                 G4ThreeVector& xx,
0117                                                 G4ThreeVector& n)
0118 {
0119    // DistanceToPlane :
0120    // Calculate distance to plane in local coordinate,
0121    // then return distance and global intersection points.
0122    // t1         - 1st. vector lying on the plane
0123    // t2         - 2nd. vector lying on the plane
0124
0125    n = (t1.cross(t2)).unit();
0126    return DistanceToPlane(p, x0, n, xx);
0127 }
0128
0129 //=====================================================================
0130 //* DistanceToLine ----------------------------------------------------
0131
0132 inline
0133 G4double G4VTwistSurface::DistanceToLine(const G4ThreeVector& p,
0134                                          const G4ThreeVector& x0,
0135                                          const G4ThreeVector& d,
0136                                                G4ThreeVector& xx)
0137 {
0138    // DistanceToLine :
0139    // Calculate distance to line,
0140    // then return distance and global intersection points.
0141    //
0142    // p          - location of flying particle
0143    // x0         - reference point of line
0144    // d          - direction vector of line
0145    // xx         - a foot of perpendicular line from p to the plane
0146    // t          - distance from xx to p
0147    //
0148    // Equation
0149    //
0150    //    distance^2 = |(xx - p)|^2
0151    //    with
0152    //       xx = x0 + t*d
0153    //
0154    //   (d/dt)distance^2 = (d/dt)|((x0 + t*d) - p)|^2
0155    //                    = 2*t*|d|^2 + 2*d*(x0 - p)
0156    //                    = 0  // smallest distance
0157    //   then
0158    //      t = - d*(x0 - p) / |d|^2
0159    //
0160
0161    G4double t;
0162    G4ThreeVector dir = d.unit();
0163    t  = - dir * (x0 - p);      // |dir|^2 = 1.
0164    xx = x0 + t * dir;
0165
0166    G4ThreeVector dist = xx - p;
0167    return dist.mag();
0168 }
0169
0170 //=====================================================================
0171 //* IsAxis0 -----------------------------------------------------------
0172
0173 inline
0174 G4bool G4VTwistSurface::IsAxis0(G4int areacode) const
0175 {
0176    return (areacode & sAxis0) != 0;
0177 }
0178
0179 //=====================================================================
0180 //* IsAxis1 -----------------------------------------------------------
0181
0182 inline
0183 G4bool G4VTwistSurface::IsAxis1(G4int areacode) const
0184 {
0185    return (areacode & sAxis1) != 0;
0186 }
0187
0188 //=====================================================================
0189 //* IsOutside ---------------------------------------------------------
0190
0191 inline
0192 G4bool G4VTwistSurface::IsOutside(G4int areacode) const
0193 {
0194    return (areacode & sInside) == 0;
0195 }
0196
0197 //=====================================================================
0198 //* IsInside ----------------------------------------------------------
0199
0200 inline
0201 G4bool G4VTwistSurface::IsInside(G4int areacode, G4bool testbitmode) const
0202 {
0203   if ((areacode & sInside) != 0)
0204   {
0205     if (testbitmode) { return true; }
0206     if (((areacode & sBoundary) == 0)
0207      && ((areacode & sCorner) == 0)) { return true; }
0208   }
0209   return false;
0210 }
0211
0212 //=====================================================================
0213 //* IsBoundary --------------------------------------------------------
0214
0215 inline
0216 G4bool G4VTwistSurface::IsBoundary(G4int areacode, G4bool testbitmode) const
0217 {
0218   if ((areacode & sBoundary) == sBoundary)
0219   {
0220     if (testbitmode) { return true; }
0221     if ((areacode & sInside) == sInside) { return true; }
0222   }
0223   return false;
0224 }
0225
0226 //=====================================================================
0227 //* IsCorner ----------------------------------------------------------
0228
0229 inline
0230 G4bool G4VTwistSurface::IsCorner(G4int areacode, G4bool testbitmode) const
0231 {
0232   if ((areacode & sCorner) == sCorner)
0233   {
0234     if (testbitmode) { return true; }
0235     if ((areacode & sInside) == sInside) { return true; }
0236   }
0237   return false;
0238 }
0239
0240 //=====================================================================
0241 //* GetAxisType -------------------------------------------------------
0242
0243 inline
0244 G4int G4VTwistSurface::GetAxisType(G4int areacode, G4int whichaxis) const
0245 {
0246   G4int axiscode = areacode & sAxisMask & whichaxis;
0247
0248   if (axiscode == (sAxisX & sAxis0) ||
0249       axiscode == (sAxisX & sAxis1))
0250   {
0251     return sAxisX;
0252   }
0253   if (axiscode == (sAxisY & sAxis0) ||
0254       axiscode == (sAxisY & sAxis1))
0255   {
0256     return sAxisY;
0257   }
0258   if (axiscode == (sAxisZ & sAxis0) ||
0259       axiscode == (sAxisZ & sAxis1))
0260   {
0261     return sAxisZ;
0262   }
0263   if (axiscode == (sAxisRho & sAxis0) ||
0264       axiscode == (sAxisRho & sAxis1))
0265   {
0266     return sAxisRho;
0267   }
0268   if (axiscode == (sAxisPhi & sAxis0) ||
0269       axiscode == (sAxisPhi & sAxis1))
0270   {
0271     return sAxisPhi;
0272   }
0273
0274   std::ostringstream message;
0275   message << "Configuration not supported." << G4endl
0276           << "        areacode = " << areacode;
0277   G4Exception("G4VTwistSurface::GetAxisType()","GeomSolids0001",
0278               FatalException, message);
0279   return 1;
0280 }
0281
0282 //=====================================================================
0283 //* ComputeGlobalPoint ------------------------------------------------
0284
0285 inline
0286 G4ThreeVector G4VTwistSurface::ComputeGlobalPoint(const G4ThreeVector& lp) const
0287 {
0288    return fRot * G4ThreeVector(lp) + fTrans;
0289 }
0290
0291 //=====================================================================
0292 //* ComputeGlobalPoint ------------------------------------------------
0293
0294 inline
0295 G4ThreeVector G4VTwistSurface::ComputeLocalPoint(const G4ThreeVector& gp) const
0296 {
0297    return fRot.inverse() * ( G4ThreeVector(gp) - fTrans ) ;
0298 }
0299
0300 //=====================================================================
0301 //* ComputeGlobalDirection --------------------------------------------
0302
0303 inline G4ThreeVector
0304 G4VTwistSurface::ComputeGlobalDirection(const G4ThreeVector& lp) const
0305 {
0306    return fRot * G4ThreeVector(lp);
0307 }
0308
0309 //=====================================================================
0310 //* ComputeLocalDirection ---------------------------------------------
0311
0312 inline G4ThreeVector
0313 G4VTwistSurface::ComputeLocalDirection(const G4ThreeVector& gp) const
0314 {
0315    return fRot.inverse() * G4ThreeVector(gp);
0316 }
0317
0318 //=====================================================================
0319 //* SetNeighbours -----------------------------------------------------
0320
0321 inline void
0322 G4VTwistSurface::SetNeighbours(G4VTwistSurface* ax0min, G4VTwistSurface* ax1min,
0323                                G4VTwistSurface* ax0max, G4VTwistSurface* ax1max)
0324 {
0325    fNeighbours[0] = ax0min;
0326    fNeighbours[1] = ax1min;
0327    fNeighbours[2] = ax0max;
0328    fNeighbours[3] = ax1max;
0329 }
0330
0331 //=====================================================================
0332 //* GetNeighbours -----------------------------------------------------
0333
0334 inline G4int
0335 G4VTwistSurface::GetNeighbours(G4int areacode, G4VTwistSurface** surfaces)
0336 {
0337
0338   G4int sAxis0Min = sAxis0 & sAxisMin ;
0339   G4int sAxis1Min = sAxis1 & sAxisMin ;
0340   G4int sAxis0Max = sAxis0 & sAxisMax ;
0341   G4int sAxis1Max = sAxis1 & sAxisMax ;
0342
0343   G4int i = 0;
0344
0345   if ( (areacode & sAxis0Min ) == sAxis0Min )
0346   {
0347     surfaces[i] = fNeighbours[0] ;
0348     ++i ;
0349   }
0350
0351   if ( ( areacode & sAxis1Min ) == sAxis1Min )
0352   {
0353      surfaces[i] = fNeighbours[1] ;
0354     ++i ;
0355     if ( i == 2 ) { return i ; }
0356   }
0357
0358   if ( ( areacode & sAxis0Max ) == sAxis0Max )
0359   {
0360     surfaces[i] = fNeighbours[2] ;
0361     ++i ;
0362     if ( i == 2 ) { return i ; }
0363   }
0364
0365   if ( ( areacode & sAxis1Max ) == sAxis1Max )
0366   {
0367     surfaces[i] = fNeighbours[3] ;
0368     ++i ;
0369     if ( i == 2 ) { return i ; }
0370   }
0371
0372   return i ;
0373 }
0374
0375 //=====================================================================
0376 //* GetCorner ---------------------------------------------------------
0377
0378 inline
0379 G4ThreeVector G4VTwistSurface::GetCorner(G4int areacode) const
0380 {
0381   if ((areacode & sCorner) == 0)
0382   {
0383     std::ostringstream message;
0384     message << "Area code must represent corner." << G4endl
0385             << "        areacode = " << areacode;
0386     G4Exception("G4VTwistSurface::GetCorner()","GeomSolids0002",
0387                 FatalException, message);
0388   }
0389
0390   if ((areacode & sC0Min1Min) == sC0Min1Min)
0391   {
0392     return fCorners[0];
0393   }
0394   if ((areacode & sC0Max1Min) == sC0Max1Min)
0395   {
0396     return fCorners[1];
0397   }
0398   if ((areacode & sC0Max1Max) == sC0Max1Max)
0399   {
0400     return fCorners[2];
0401   }
0402   if ((areacode & sC0Min1Max) == sC0Min1Max)
0403   {
0404     return fCorners[3];
0405   }
0406
0407   std::ostringstream message;
0408   message << "Configuration not supported." << G4endl
0409           << "        areacode = " << areacode;
0410   G4Exception("G4VTwistSurface::GetCorner()", "GeomSolids0001",
0411               FatalException, message);
0412
0413   return fCorners[0];
0414 }