BRAT 2.4.5

//____________________________________________________________________
// 
//  BrLine3D
//
//  Line class for Geometry. A line is defined by an origin and a unit
//  vector defining its direction.
//

//____________________________________________________________________
//
// $Id: BrLine3D.cxx,v 1.3 2001/11/16 20:28:53 hagel Exp $
// $Author: hagel $
// $Date: 2001/11/16 20:28:53 $
// $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov>
//
#ifndef BRAT_BrLine3D
#include "BrLine3D.h"
#endif

#ifndef ROOT_TMath
#include "TMath.h"
#endif
#ifndef BRAT_BrModule
#include "BrModule.h"
#endif

//____________________________________________________________________
ClassImp(BrLine3D);

//____________________________________________________________________
 BrLine3D::BrLine3D()
  : TObject()
{
  // Default constructor - empty
}

//____________________________________________________________________
 BrLine3D::BrLine3D(const BrVector3D& origin, const BrVector3D& direction)
  : TObject()
{
  //
  //  Constructor
  //
  fOrigin = origin;
  fDirection = direction.Unit();
}

//  //____________________________________________________________________
//  BrLine3D::BrLine3D(const BrLine3D& line) 
//    : fOrigin(line.GetOrigin()), fDirection(line.GetDirection()), TObject()
//  {
//    //
//    //  Copy Constructor
//    //
//  }

//  //____________________________________________________________________
//  BrLine3D::~BrLine3D()
//  {
//    // Dtor - empty
//  }

//____________________________________________________________________
BrLine3D& BrLine3D::operator = (const BrLine3D& line){
  //
  // assignement operator
  //
  fOrigin = line.fOrigin;
  fDirection = line.fDirection;
  return *this;
}

//____________________________________________________________________
 BrVector3D BrLine3D::ClosestPoint(const BrVector3D& vector) const{
  //
  // Return the closest point on the line to the vector
  //
  Double_t value;
  BrVector3D vtmp;
  vtmp = vector - fOrigin;
  value = vtmp.Dot(fDirection);
  vtmp = fOrigin+value*fDirection;
  return vtmp;
}

//____________________________________________________________________
 Double_t BrLine3D::Distance(const BrVector3D& vector) const{
  //
  // distance from a vector (i.e. point) to this line.
  //
  Double_t value, value1;
  BrVector3D vtmp;
  vtmp = vector - fOrigin;
  value1 = vtmp.Dot(fDirection);
  value = vtmp.NormSq() - value1*value1; 
  return TMath::Sqrt(value);
}

//____________________________________________________________________
 BrLine3D BrLine3D::GetShortestLineBetween(const BrLine3D &line) const {
  //
  // Return the shortest line between the two lines.  The length of
  // the vector part is the shortest distance between the two lines.
  // put in according to 
  // http://www.mhri.edu.au/~pdb/geometry/lineline3d/ 1/31/99
  // but not tested yet?  Be careful before using in production mode.
  // dmnop = (xm-xn)(xo-xp) + (ym-yn)(yo-yp) + (zm-zn)(zo-zp)
  Double_t x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4;
  Double_t d1343,d4321,d1321,d4343,d2121;
  Double_t ua,ub;
  x1 = fOrigin[0];
  x2 = x1 + fDirection[0];
  y1 = fOrigin[1];
  y2 = y1 + fDirection[1];
  z1 = fOrigin[2];
  z2 = z1 + fDirection[2];
  x3 = line.GetOrigin()[0];
  x4 = x3 + line.GetDirection()[0];
  y3 = line.GetOrigin()[1];
  y4 = y3 + line.GetDirection()[1];
  z3 = line.GetOrigin()[2];
  z4 = z3 + line.GetDirection()[2];
  d1343 = (x1-x3)*(x4-x3) + (y1-y3)*(y4-y3) + (z1-z3)*(z4-z3);
  d4321 = (x4-x3)*(x2-x1) + (y4-y3)*(y2-y1) + (z4-z3)*(z2-z1);
  d1321 = (x1-x3)*(x2-x1) + (y1-y3)*(y2-y1) + (z1-z3)*(z2-z1);
  d4343 = (x4-x3)*(x4-x3) + (y4-y3)*(y4-y3) + (z4-z3)*(z4-z3);
  d2121 = (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1);
  ua = (d1343*d4321 - d1321*d4343) / (d2121*d4343 - d4321*d4321);
  ub = (d1343 + ua * d4321) / d4343;
  BrVector3D Pa = fOrigin + ua*fDirection.Unit();
  BrVector3D Pb = line.GetOrigin() + ub * line.GetDirection().Unit();
  BrLine3D shortline(Pa,Pa-Pb);
  return shortline;
}

//____________________________________________________________________
 Double_t BrLine3D::GetShortestDistanceBetween(const BrLine3D &line) const {
  //Get shortest distance between the two lines.  This should give the
  //same answer as the modulus of GetShortestLineBetween().  Coded
  //1-31-99 but not tested.
  BrVector3D N1 = fDirection;
  BrVector3D N2 = line.GetDirection();
  BrVector3D normal = N1.Cross(N2);

  const Double_t small = 1.e-10;
  //Check whether lines are parallel to each other by looking at
  //magnitude of cross product.
  if(TMath::Abs(normal.Norm()) < small) {
     //This means lines are parallel to each other
     //So get perpendicular distance between origin of line 1 to the other line
     return Distance(line.GetOrigin());
     }

  //BrVector3D unit = N1.Cross(N2).Unit();
  BrVector3D unit = normal.Unit();

  //should give vector between origins with proper length
  BrVector3D temp = fOrigin-line.GetOrigin();  
  return temp.Dot(unit);
}


//____________________________________________________________________
 void BrLine3D::Print(Option_t *option = "") const
{
  //
  // Simple output of this line
  //
  cout << "Origin   : ";
  fOrigin.Print();
  cout << "Diection : ";
  fDirection.Print();
}

//____________________________________________________________________
 Float_t BrLine3D::RelativeOverlap(const BrLine3D& line, Float_t zl, Float_t r)
{
  //
  // Calculate the relative overlap between "this" track and the given track
  //
  // Parameters:
  //   In:  track   pointer to the other track
  //        zl      length of detector (in local z direction)
  //        r       radius of cylinder around tracks
  //
  // Algorithm:
  //   A cylinder of radius r is put around each track.
  //   The relative overlap is defined as the ratio between
  //     the volume of the overlap region of two cylinders
  //   and
  //     the volume of the cylinder around "this" track is calculated.
  //   The cylinders are distorted in order to have a circular intersection
  //   between the cylinders and planes at arbitrary z.
  //
  // Requires:
  //   The tracks are parametrized as:
  //     x = ax * z + bx
  //     y = ay * z + by
  //   
  
  const Double_t slopeZthis = fDirection[2];
  const Double_t slopeZline = line.GetDirection()[2];
  if(slopeZthis == 0 || slopeZline == 0) {
    Warning("RelativeOverlap", "z slope == 0");
    return 0;
  }

  const Float_t twor  = 2*r;
  const Float_t halfl = zl/2;
  
  // Difference between the two tracks
  // Note that we normalise to slope z == 1
  const Float_t dax = fDirection[0]/slopeZthis - 
    line.GetDirection()[0]/slopeZline;
  const Float_t day = fDirection[1]/slopeZthis -
    line.GetDirection()[1]/slopeZline;
  const Float_t dbx = fOrigin[0] - line.GetOrigin()[0];
  const Float_t dby = fOrigin[1] - line.GetOrigin()[1];
  
  Float_t overlap = 0;
     
  const Float_t a = dax*dax + day*day;

  if (a == 0) {
    
    // Tracks are parallel
    const Float_t d = TMath::Sqrt(dbx*dbx+dby*dby)/twor;
    if (d < 1)
      overlap = TMath::ACos(d)-d*TMath::Sqrt(1-d*d);
  } else {

    const Float_t b = (dax*dbx+day*dby)/a;
    const Float_t c = (dbx*dbx+dby*dby)/a;
    Float_t d = b*b-c+twor*twor/a;
    if (d > 0) {
      // Cylinders intersect
      d = TMath::Sqrt(d);
      Float_t zn = -b - d;
      Float_t zx = -b + d;
      if (zn < -halfl) zn = -halfl;
      if (zx >  halfl) zx =  halfl;
      if (zn < zx) {

	// Cylinders intersect inside the detector
	//
	// Integrate
	//   acos(d(z)/(2r))-d(z)/(2r)*sqrt(1-(d(z)/(2r))**2)
	// where
	//   d(z)/(2r)=sqrt(a*((z+2*b)*z+c))/(2*r)
	// from zn to zx
	// 
	// Quick and dirty implementation: Simple N-point rule
	const Int_t   n  = 5;
	const Float_t dz = (zx-zn)/n;
	for (Int_t i = 0; i < n; i++) {

	  Float_t z = zn + dz*(i+0.5);
	  d = TMath::Sqrt(a*((z+2*b)*z+c))/twor;
	  overlap += TMath::ACos(d)-d*TMath::Sqrt(1-d*d);
	}

	overlap *= dz/zl;
      }
    }
  }
  
  return overlap/(TMath::Pi()/2);
}

//____________________________________________________________________
 void BrLine3D::SetDirection(const BrVector3D& vector)
{
  // set x, y, z
  const Double_t norm = vector.Norm();
  fDirection.SetX(vector.GetX() / norm);
  fDirection.SetY(vector.GetY() / norm);
  fDirection.SetZ(vector.GetZ() / norm);
}

//____________________________________________________________________
 void BrLine3D::SetOrigin(const BrVector3D& vector)
{
  // set x, y, z
  fOrigin.SetX(vector.GetX());
  fOrigin.SetY(vector.GetY());
  fOrigin.SetZ(vector.GetZ());
}

//____________________________________________________________________
ostream& operator<<(ostream& os, const BrLine3D& line) {
  //
  // Related Function (global non-member functions)
  // Output line parameters in default layout
  //
  os << line.GetOrigin() << "+ t*"<< line.GetDirection(); 
  return os; 
}

//____________________________________________________________________
ostream& operator<<(ostream& os, const BrLine3D* line) {
  os<<*line;
}

//____________________________________________________________________
//
// $Log: BrLine3D.cxx,v $
// Revision 1.3  2001/11/16 20:28:53  hagel
// Implemented iostream for BrLine3D
//
// Revision 1.2  2001/08/14 19:27:52  pchristi
// Removed copy constructors and destructors where the default was good enough.
//
// Revision 1.1.1.1  2001/06/21 14:55:20  hagel
// Initial revision of brat2
//
// Revision 1.4  2001/06/17 17:10:08  pchristi
// Cleaned up and modified the line 3D and vector3D classes.
// No major new stuff.
//
//
//  $Id: BrLine3D.cxx,v 1.3 2001/11/16 20:28:53 hagel Exp $
//
//  $Log: BrLine3D.cxx,v $
//  Revision 1.3  2001/11/16 20:28:53  hagel
//  Implemented iostream for BrLine3D
//
//  Revision 1.2  2001/08/14 19:27:52  pchristi
//  Removed copy constructors and destructors where the default was good enough.
//
//  Revision 1.1.1.1  2001/06/21 14:55:20  hagel
//  Initial revision of brat2
//
//  Revision 1.4  2001/06/17 17:10:08  pchristi
//  Cleaned up and modified the line 3D and vector3D classes.
//  No major new stuff.
//
//  Revision 1.3  2000/12/20 01:56:41  hagel
//  Added missing end-of-line per gcc-2.96 complaint
//
//  Revision 1.2  1999/01/31 21:41:22  hagel
//  Added Methods to get shortest distance between two lines
//
//  Revision 1.1  1998/12/01 20:50:11  videbaek
//  Added several classes to Geometry. Updated source and makefile(s)
//  - BrLine3D
//  - BrCoordinateSystem
//  - BrRotMatrix
//
This page automatically generated by script docBrat by Christian Holm
Validate HTML
Validate CSS