PhotonUtil.cxx

Go to the documentation of this file.

#include "PhotonUtil.h"
#include <TMath.h>

#include "MessageService/MsgService.h"

CVSID(" $Id: PhotonUtil.cxx,v 1.5 2005/02/02 16:44:33 tagg Exp $ " );

ClassImp(PhotonUtil)

Bool_t PhotonUtil::Refract(const TVector3& v1,
                           Double_t n1, 
                           Double_t n2, 
                           const TVector3& norm,
                           TVector3& v2)
{
  // 3-d version of snell's law. 
  // Refracts ray with direction v1 as it goes through surface with normal norm
  // if v1 is in medium with refractive index n1 and v2 is in refractive
  // index n2. Returns v2, the final direction.
  // v1 and norm must be unit vectors.
  // Return true for refraction, false for total internal reflection.

  // Decompose v1 into two components: along normal, and in the plane of the
  // refraction:
  double ns = 1;
  double b1 = norm.Dot(v1);
  if(b1<0) {
    b1 = -b1;
    ns = -1;
  };

  if(b1>1.0) {
    MSG("Photon",Msg::kWarning) << "Bad arguments to Refract()" << endl;
    return false;
  }
  double a1 = sqrt(1.0-b1*b1);

  // The vector tangent to the surface, and in the plane of the refraction:
  TVector3        a = ((1./a1)*v1) - (((b1/a1)*ns)*norm);
  //const TVector3& b = norm; // For completeness.

  // After refraction:
  double a2 = a1*n1/n2;
  if(a2>=1.0) {
    // Reflection.
    // For completeness, actually do the reflection:

    v2 = v1 - 2.0*b1*ns*norm;
    return false;
  }
  
  double b2 = sqrt(1.0-a2*a2);
  
  // Finally, compose the final direction:
  v2 = a2*a + (b2*ns)*norm;
  return true;
}


Bool_t PhotonUtil::IntersectCylinder(const TVector3& x,  // Position, with cylinder centered on x-axis
                                     const TVector3& v,  // Direction
                                     Double_t R, // Radius.
                                     Double_t &d1, // Distance to first intersection
                                     Double_t &d2  // Distance to second intesection
                                     )
{
  // Find the intersection of a line with a cylinder
  // given that the cylinder is radius R parallel to the x-axis
  // centered on the origin.
  // Sets d1 and d2, the distances to the first and second positions respectively.
  // d2-d1 is the chord.
  // Return false if there is no intersection.
  
  d1=d2=0;

  double v2  = (v.y()*v.y() + v.z()*v.z());
  if(v2==0) return false; // going exactly along the x-axis.
  
  double vdotdx = (v.y() * x.y()) + (v.z() * x.z());
  double insqrt =  (vdotdx*vdotdx) - 
    (v2)*((x.y()*x.y()) + (x.z()*x.z()) - (R*R));

  // The track does not intersect the fibre ever.
  if(insqrt<0) return false;

  double a = -vdotdx/v2;
  double b = TMath::Sqrt(insqrt)/v2;

  // the first solution:
  d1 = a-b;

  // the second solution:
  d2 = a+b;

  return true;
}


Int_t PhotonUtil::InteractWithCylinder(  TVector3& x, // pos (at cylinder!)
                                          TVector3& v, // direction
                                          Double_t R,  // cylinder radius
                                          Double_t n1,
                                          Double_t n2)
{
  // Move the photon along v to the first intesection point.
  // Reflect or refract off the surface.
  // Return 0 for refraction, 1 for reflection, -1 for no hit.
  // (i.e. return 0 if it goes inside)
  // Change x and v to match the new position, direction.

  double d1, d2;
  bool intersect = PhotonUtil::IntersectCylinder( x, 
                                                  v, 
                                                  R, 
                                                  d1,d2);

  if(!intersect) return -1;
  
  // Otherwise, move to the closest intercept point.
  // If we are on the cylinder or have not yet hit the cylinder, this 
  // should always be the d1.
  x += v*d1;

  // x is now at the surface of the cylinder. Find the surface 
  // normal at this point.

  TVector3 normal(0,x.y(),x.z());

  if(normal.Mag2()==0.0) {
    // Hmm... just skimmed it.
    // Assume it's a reflection and return nothing.
    return 1;
  }
  normal.SetMag(1.0);
  
  bool refract = PhotonUtil::Refract(v,n1,n2,normal,v);
  if(refract) return 0;

  return 1; // Reflection
}

---