G4IntersectingCone.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: G4IntersectingCone.cc 67011 2013-01-29 16:17:41Z gcosmo $
00028 //
00029 // 
00030 // --------------------------------------------------------------------
00031 // GEANT 4 class source file
00032 //
00033 //
00034 // G4IntersectingCone.cc
00035 //
00036 // Implementation of a utility class which calculates the intersection
00037 // of an arbitrary line with a fixed cone
00038 // --------------------------------------------------------------------
00039 
00040 #include "G4IntersectingCone.hh"
00041 #include "G4GeometryTolerance.hh"
00042 
00043 //
00044 // Constructor
00045 //
00046 G4IntersectingCone::G4IntersectingCone( const G4double r[2],
00047                                         const G4double z[2] )
00048 { 
00049   static const G4double halfCarTolerance
00050     = 0.5 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00051 
00052   //
00053   // What type of cone are we?
00054   //
00055   type1 = (std::fabs(z[1]-z[0]) > std::fabs(r[1]-r[0]));
00056   
00057   if (type1)
00058   {
00059     B = (r[1]-r[0])/(z[1]-z[0]);      // tube like
00060     A = 0.5*( r[1]+r[0] - B*(z[1]+z[0]) );
00061   }
00062   else
00063   {
00064     B = (z[1]-z[0])/(r[1]-r[0]);      // disk like
00065     A = 0.5*( z[1]+z[0] - B*(r[1]+r[0]) );
00066   }
00067   //
00068   // Calculate extent
00069   //
00070   if (r[0] < r[1])
00071   {
00072     rLo = r[0]-halfCarTolerance; rHi = r[1]+halfCarTolerance;
00073   }
00074   else
00075   {
00076     rLo = r[1]-halfCarTolerance; rHi = r[0]+halfCarTolerance;
00077   }
00078   
00079   if (z[0] < z[1])
00080   {
00081     zLo = z[0]-halfCarTolerance; zHi = z[1]+halfCarTolerance;
00082   }
00083   else
00084   {
00085     zLo = z[1]-halfCarTolerance; zHi = z[0]+halfCarTolerance;
00086   }
00087 }
00088 
00089 
00090 //
00091 // Fake default constructor - sets only member data and allocates memory
00092 //                            for usage restricted to object persistency.
00093 //
00094 G4IntersectingCone::G4IntersectingCone( __void__& )
00095   : zLo(0.), zHi(0.), rLo(0.), rHi(0.), type1(false), A(0.), B(0.)
00096 {
00097 }
00098 
00099 
00100 //
00101 // Destructor
00102 //
00103 G4IntersectingCone::~G4IntersectingCone()
00104 {
00105 }
00106 
00107 
00108 //
00109 // HitOn
00110 //
00111 // Check r or z extent, as appropriate, to see if the point is possibly
00112 // on the cone.
00113 //
00114 G4bool G4IntersectingCone::HitOn( const G4double r,
00115                                   const G4double z )
00116 {
00117   //
00118   // Be careful! The inequalities cannot be "<=" and ">=" here without
00119   // punching a tiny hole in our shape!
00120   //
00121   if (type1)
00122   {
00123     if (z < zLo || z > zHi) return false;
00124   }
00125   else
00126   {
00127     if (r < rLo || r > rHi) return false;
00128   }
00129 
00130   return true;
00131 }
00132 
00133 
00134 //
00135 // LineHitsCone
00136 //
00137 // Calculate the intersection of a line with our conical surface, ignoring
00138 // any phi division
00139 //
00140 G4int G4IntersectingCone::LineHitsCone( const G4ThreeVector &p,
00141                                         const G4ThreeVector &v,
00142                                               G4double *s1, G4double *s2 )
00143 {
00144   if (type1)
00145   {
00146     return LineHitsCone1( p, v, s1, s2 );
00147   }
00148   else
00149   {
00150     return LineHitsCone2( p, v, s1, s2 );
00151   }
00152 }
00153 
00154 
00155 //
00156 // LineHitsCone1
00157 //
00158 // Calculate the intersections of a line with a conical surface. Only
00159 // suitable if zPlane[0] != zPlane[1].
00160 //
00161 // Equation of a line:
00162 //
00163 //       x = x0 + s*tx      y = y0 + s*ty      z = z0 + s*tz
00164 //
00165 // Equation of a conical surface:
00166 //
00167 //       x**2 + y**2 = (A + B*z)**2
00168 //
00169 // Solution is quadratic:
00170 //
00171 //  a*s**2 + b*s + c = 0
00172 //
00173 // where:
00174 //
00175 //  a = x0**2 + y0**2 - (A + B*z0)**2
00176 //
00177 //  b = 2*( x0*tx + y0*ty - (A*B - B*B*z0)*tz)
00178 //
00179 //  c = tx**2 + ty**2 - (B*tz)**2
00180 //
00181 // Notice, that if a < 0, this indicates that the two solutions (assuming
00182 // they exist) are in opposite cones (that is, given z0 = -A/B, one z < z0
00183 // and the other z > z0). For our shapes, the invalid solution is one
00184 // which produces A + Bz < 0, or the one where Bz is smallest (most negative).
00185 // Since Bz = B*s*tz, if B*tz > 0, we want the largest s, otherwise,
00186 // the smaller.
00187 //
00188 // If there are two solutions on one side of the cone, we want to make
00189 // sure that they are on the "correct" side, that is A + B*z0 + s*B*tz >= 0.
00190 //
00191 // If a = 0, we have a linear problem: s = c/b, which again gives one solution.
00192 // This should be rare.
00193 //
00194 // For b*b - 4*a*c = 0, we also have one solution, which is almost always
00195 // a line just grazing the surface of a the cone, which we want to ignore. 
00196 // However, there are two other, very rare, possibilities:
00197 // a line intersecting the z axis and either:
00198 //       1. At the same angle std::atan(B) to just miss one side of the cone, or
00199 //       2. Intersecting the cone apex (0,0,-A/B)
00200 // We *don't* want to miss these! How do we identify them? Well, since
00201 // this case is rare, we can at least swallow a little more CPU than we would
00202 // normally be comfortable with. Intersection with the z axis means
00203 // x0*ty - y0*tx = 0. Case (1) means a==0, and we've already dealt with that
00204 // above. Case (2) means a < 0.
00205 //
00206 // Now: x0*tx + y0*ty = 0 in terms of roundoff error. We can write:
00207 //             Delta = x0*tx + y0*ty
00208 //             b = 2*( Delta - (A*B + B*B*z0)*tz )
00209 // For:
00210 //             b*b - 4*a*c = epsilon
00211 // where epsilon is small, then:
00212 //             Delta = epsilon/2/B
00213 // 
00214 G4int G4IntersectingCone::LineHitsCone1( const G4ThreeVector &p,
00215                                          const G4ThreeVector &v,
00216                                                G4double *s1, G4double *s2 )
00217 {
00218   G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
00219   G4double tx = v.x(), ty = v.y(), tz = v.z();
00220 
00221   G4double a = tx*tx + ty*ty - sqr(B*tz);
00222   G4double b = 2*( x0*tx + y0*ty - (A*B + B*B*z0)*tz);
00223   G4double c = x0*x0 + y0*y0 - sqr(A + B*z0);
00224   
00225   G4double radical = b*b - 4*a*c;
00226  
00227   if (radical < -1E-6*std::fabs(b))  { return 0; }    // No solution
00228   
00229   if (radical < 1E-6*std::fabs(b))
00230   {
00231     //
00232     // The radical is roughly zero: check for special, very rare, cases
00233     //
00234     if (std::fabs(a) > 1/kInfinity)
00235       {
00236       if(B==0.) { return 0; }
00237       if ( std::fabs(x0*ty - y0*tx) < std::fabs(1E-6/B) )
00238       {
00239          *s1 = -0.5*b/a;
00240          return 1;
00241       }
00242       return 0;
00243     }
00244   }
00245   else
00246   {
00247     radical = std::sqrt(radical);
00248   }
00249   
00250   if (a > 1/kInfinity)
00251   {
00252     G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
00253     sa = q/a;
00254     sb = c/q;
00255     if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
00256     if (A + B*(z0+(*s1)*tz) < 0)  { return 0; }
00257     return 2;
00258   }
00259   else if (a < -1/kInfinity)
00260   {
00261     G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
00262     sa = q/a;
00263     sb = c/q;
00264     *s1 = (B*tz > 0)^(sa > sb) ? sb : sa;
00265     return 1;
00266   }
00267   else if (std::fabs(b) < 1/kInfinity)
00268   {
00269     return 0;
00270   }
00271   else
00272   {
00273     *s1 = -c/b;
00274     if (A + B*(z0+(*s1)*tz) < 0)  { return 0; }
00275     return 1;
00276   }
00277 }
00278 
00279   
00280 //
00281 // LineHitsCone2
00282 //
00283 // See comments under LineHitsCone1. In this routine, case2, we have:
00284 //
00285 //   Z = A + B*R
00286 //
00287 // The solution is still quadratic:
00288 //
00289 //  a = tz**2 - B*B*(tx**2 + ty**2)
00290 //
00291 //  b = 2*( (z0-A)*tz - B*B*(x0*tx+y0*ty) )
00292 //
00293 //  c = ( (z0-A)**2 - B*B*(x0**2 + y0**2) )
00294 //
00295 // The rest is much the same, except some details.
00296 //
00297 // a > 0 now means we intersect only once in the correct hemisphere.
00298 //
00299 // a > 0 ? We only want solution which produces R > 0. 
00300 // since R = (z0+s*tz-A)/B, for tz/B > 0, this is the largest s
00301 //          for tz/B < 0, this is the smallest s
00302 // thus, same as in case 1 ( since sign(tz/B) = sign(tz*B) )
00303 //
00304 G4int G4IntersectingCone::LineHitsCone2( const G4ThreeVector &p,
00305                                          const G4ThreeVector &v,
00306                                                G4double *s1, G4double *s2 )
00307 {
00308   G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
00309   G4double tx = v.x(), ty = v.y(), tz = v.z();
00310   
00311   
00312   // Special case which might not be so rare: B = 0 (precisely)
00313   //
00314   if (B==0)
00315   {
00316     if (std::fabs(tz) < 1/kInfinity)  { return 0; }
00317     
00318     *s1 = (A-z0)/tz;
00319     return 1;
00320   }
00321 
00322   G4double B2 = B*B;
00323 
00324   G4double a = tz*tz - B2*(tx*tx + ty*ty);
00325   G4double b = 2*( (z0-A)*tz - B2*(x0*tx + y0*ty) );
00326   G4double c = sqr(z0-A) - B2*( x0*x0 + y0*y0 );
00327   
00328   G4double radical = b*b - 4*a*c;
00329  
00330   if (radical < -1E-6*std::fabs(b)) { return 0; }   // No solution
00331   
00332   if (radical < 1E-6*std::fabs(b))
00333   {
00334     //
00335     // The radical is roughly zero: check for special, very rare, cases
00336     //
00337     if (std::fabs(a) > 1/kInfinity)
00338     {
00339       if ( std::fabs(x0*ty - y0*tx) < std::fabs(1E-6/B) )
00340       {
00341         *s1 = -0.5*b/a;
00342         return 1;
00343       }
00344       return 0;
00345     }
00346   }
00347   else
00348   {
00349     radical = std::sqrt(radical);
00350   }
00351   
00352   if (a < -1/kInfinity)
00353   {
00354     G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
00355     sa = q/a;
00356     sb = c/q;
00357     if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
00358     if ((z0 + (*s1)*tz  - A)/B < 0)  { return 0; }
00359     return 2;
00360   }
00361   else if (a > 1/kInfinity)
00362   {
00363     G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
00364     sa = q/a;
00365     sb = c/q;
00366     *s1 = (tz*B > 0)^(sa > sb) ? sb : sa;
00367     return 1;
00368   }
00369   else if (std::fabs(b) < 1/kInfinity)
00370   {
00371     return 0;
00372   }
00373   else
00374   {
00375     *s1 = -c/b;
00376     if ((z0 + (*s1)*tz  - A)/B < 0)  { return 0; }
00377     return 1;
00378   }
00379 }

Generated on Mon May 27 17:48:38 2013 for Geant4 by  doxygen 1.4.7