MucGeometron.cxx

Go to the documentation of this file.

//$id$
//
//$log$
 
/*
 *    2003/08/30   Zhengyun You     Peking University
 *
 *    2004/09/09   Zhengyun You     Peking University
 *                 transplanted to Gaudi framework
 */
 
using namespace std;
 
#include "MucGeomSvc/MucGeometron.h"
 
MucGeometron::MucGeometron() {
  // Constructor.
}
 
MucGeometron::~MucGeometron() {
  // Destructor.
}
 
bool MucGeometron::GetIntersectionLinePlane( const HepPoint3D pLine, const Hep3Vector vectLine,
                                             const HepPlane3D plane, HepPoint3D& cross ) {
  // Given a line by basepoint linePoint and direction lineDir,
  // find the intersection with the plane.
  HepPoint3D linePoint = pLine;
  Hep3Vector lineDir   = vectLine;
 
  HepPoint3D planePoint = plane.point();
  Hep3Vector planeDir( plane.a(), plane.b(), plane.c() );
 
  lineDir.setMag( 1.0 );
  planeDir.setMag( 1.0 );
 
  // Vector connecting the basepoint.
  HepPoint3D basePoint = planePoint - linePoint;
 
  double distance, denominator; // d , cos(theta).
 
  denominator = planeDir.dot( lineDir );
 
  if ( denominator != 0 )
  {
    // Line and gap are not parallel.
    distance = ( planeDir.dot( basePoint ) ) / denominator;
    cross    = linePoint + lineDir * distance;
    return 1;
  }
  else
  {
    // Line and gap are parrellel.
    // cout << " Line is parrellel to plane! No intersect point found!" << endl;
    return 0;
  }
}
 
bool MucGeometron::GetIntersectionLinePlaneWithSigma(
    const HepPoint3D pLine, const Hep3Vector vectLine, const HepPoint3D pLineSigma,
    const Hep3Vector vectLineSigma, const HepPlane3D plane, HepPoint3D& cross,
    HepPoint3D& crossSigma ) {
  HepPoint3D linePoint = pLine;
  Hep3Vector lineDir   = vectLine;
 
  HepPoint3D planePoint = plane.point();
  Hep3Vector planeDir( plane.a(), plane.b(), plane.c() );
 
  lineDir.setMag( 1.0 );
  planeDir.setMag( 1.0 );
  HepPoint3D basePoint = planePoint - linePoint;
 
  double distance, denominator; // d , cos(theta).
 
  denominator = planeDir.dot( lineDir );
  if ( denominator != 0 )
  {
    distance = ( planeDir.dot( basePoint ) ) / denominator;
    cross    = linePoint + lineDir * distance;
 
    // calculation of the uncertainty in intercept
 
    double x0 = pLine.x();
    double y0 = pLine.y();
    double z0 = pLine.z();
 
    double vx = vectLine.x();
    double vy = vectLine.y();
    double vz = vectLine.z();
 
    double dx0 = pLineSigma.x();
    double dy0 = pLineSigma.y();
    double dz0 = pLineSigma.z();
 
    double dvx = vectLineSigma.x();
    double dvy = vectLineSigma.y();
    double dvz = vectLineSigma.z();
 
    double pa = plane.a();
    double pb = plane.b();
    double pc = plane.c();
 
    double Const1 = planeDir.dot( planePoint );
    double Const2 = pa * vx + pb * vy + pc * vz;
    double Const3 = pa * x0 + pb * y0 + pc * z0;
 
    double sigmaX_1 = dx0 * dx0 * ( pb * vy + pc * vz ) / Const2;
    double sigmaX_2 = ( -1 ) * dy0 * dy0 * pb * vx / Const2;
    double sigmaX_3 = ( -1 ) * dz0 * dz0 * pc * vx / Const2;
    double sigmaX_4 =
        dvx * dvx * ( Const1 - Const3 ) * ( pb * vy + pc * vz ) / Const2 / Const2;
    double sigmaX_5 = ( -1 ) * dvy * dvy * vx * ( Const1 - Const3 ) * pb / Const2 / Const2;
    double sigmaX_6 = ( -1 ) * dvz * dvz * vx * ( Const1 - Const3 ) * pc / Const2 / Const2;
 
    double sigmaX = sigmaX_1 + sigmaX_2 + sigmaX_3 + sigmaX_4 + sigmaX_5 + sigmaX_6;
 
    double sigmaY_1 = ( -1 ) * dx0 * dx0 * pa * vy / Const2;
    double sigmaY_2 = dy0 * dy0 * ( 1 - pb * vy / Const2 );
    double sigmaY_3 = ( -1 ) * dz0 * dz0 * pc * vy / Const2;
    double sigmaY_4 = ( -1 ) * dvx * dvx * ( Const1 - Const3 ) * pa * vy / Const2 / Const2;
    double sigmaY_5 =
        dvy * dvy * ( Const1 - Const3 ) * ( pa * vx + pc * vz ) / Const2 / Const2;
    double sigmaY_6 = ( -1 ) * dvz * dvz * ( Const1 - Const3 ) * pc * vy / Const2 / Const2;
 
    double sigmaY = sigmaY_1 + sigmaY_2 + sigmaY_3 + sigmaY_4 + sigmaY_5 + sigmaY_6;
 
    double sigmaZ_1 = ( -1 ) * dx0 * dx0 * pa * vz / Const2;
    double sigmaZ_2 = ( -1 ) * dy0 * dy0 * pb * vz / Const2;
    double sigmaZ_3 = dz0 * dz0 * ( 1 - pc * vz / Const2 );
    double sigmaZ_4 = ( -1 ) * dvx * dvx * ( Const1 - Const3 ) * pa * vz / Const2 / Const2;
    double sigmaZ_5 = ( -1 ) * dvy * dvy * ( Const1 - Const3 ) * pb * vz / Const2 / Const2;
    double sigmaZ_6 =
        dvz * dvz * ( Const1 - Const3 ) * ( pa * vx + pb * vy ) / Const2 / Const2;
 
    double sigmaZ = sigmaZ_1 + sigmaZ_2 + sigmaZ_3 + sigmaZ_4 + sigmaZ_5 + sigmaZ_6;
 
    HepPoint3D c1( sqrt( abs( sigmaX ) ), sqrt( abs( sigmaY ) ), sqrt( abs( sigmaZ ) ) );
    crossSigma = c1;
    return 1;
  }
  else { return 0; }
}
 
bool MucGeometron::GetIntersectionQuadPlaneLocal( const int   part, // liangyt 2009.3.12
                                                  const int   orient,
                                                  const float a, // y = a * x * x + b * x + c;
                                                  const float b, const float c,
                                                  const HepPlane3D plane, HepPoint3D& cross1,
                                                  HepPoint3D& cross2 ) {
 
  if ( a == 0 )
  {
    // not quad
    return 0;
  }
 
  HepPoint3D planePoint = plane.point();
  Hep3Vector planeDir( plane.a(), plane.b(), plane.c() );
 
  float A, B, C, D;
  A = plane.a();
  B = plane.b();
  C = plane.c();
  D = planeDir.dot( planePoint );
 
  // B*a*x^2 + (Bb + A)*x + Bc - D = 0
 
  // cout<<"in geomtron ABCD = "<<A<<" "<<B<<" "<<C<<" "<<D<<endl;
 
  float a1 = B * a;
  float b1 = B * b + A;
  float c1 = B * c - D;
 
  float b2sub4ac = b1 * b1 - 4 * a1 * c1;
 
  // cout<<"in geomtron abc delta = "<<a1<<" "<<b1<<" "<<c1<<" "<<b2sub4ac<<endl;
 
  float x1 = 0.0, x2 = 0.0, y1 = 0.0, y2 = 0.0, z1 = 0.0, z2 = 0.0;
  if ( abs( a1 ) > 10E-10 )
  {
    if ( b2sub4ac >= 0 )
    {
      x1 = ( -b1 + sqrt( b2sub4ac ) ) / ( 2 * a1 );
      x2 = ( -b1 - sqrt( b2sub4ac ) ) / ( 2 * a1 );
      y1 = a * x1 * x1 + b * x1 + c;
      y2 = a * x2 * x2 + b * x2 + c;
      // cout<<"in MucGeometron------ x1 y1 z1 = "<<x1<<" "<<y1<<" "<<z1<<endl;
      // cout<<"in MucGeometron------ x2 y2 z2 = "<<x2<<" "<<y2<<" "<<z2<<endl;
      cross1.set( x1, y1, z1 );
      cross2.set( x2, y2, z2 );
      return 1;
    }
    else
    {
      // cout<<"MucGeometron: no intersection!"<<endl;
      cross1.set( -9999, -9999, -9999 );
      cross2.set( -9999, -9999, -9999 );
      return 0;
    }
  }
  else
  {
    x1 = -c1 / b1;
    y1 = a * x1 * x1 + b * x1 + c;
    // cout<<"in MucGeometron------ x1 y1 z1 = "<<x1<<" "<<y1<<" "<<z1<<endl;
    cross1.set( x1, y1, z1 );
    cross2.set( x1, y1, z1 );
    // cout<<"in MucGeometron: abs(a1)<10E-10"<<endl;
    return 1;
  }
}
 
bool MucGeometron::GetIntersectionQuadPlane( const HepPoint3D pLine, // liangyt 2007.4.9
                                             const float vy, const float y0,
                                             const float a, // y = a * x * x + b * x + c;
                                             const float b, const float c,
                                             const HepPlane3D plane, HepPoint3D& cross1,
                                             HepPoint3D& cross2 ) {
  // Given a parabola  by a b c,
  // find the intersection with the plane.
 
  // y = vy * z + y0;
  // y = a * x * x + b * x + c;
  // plane.a() * x + plane.b() * y + plane.c() * z = planeDir.dot(planePoint)
 
  HepPoint3D planePoint = plane.point();
  Hep3Vector planeDir( plane.a(), plane.b(), plane.c() );
 
  float A, B, C, D;
  A = plane.a();
  B = plane.b();
  C = plane.c();
  D = planeDir.dot( planePoint );
 
  //(B+Cv)*a*x^2 + ((B+Cv)b + A)*x + (B+Cv)c + Cy0 - D = 0
 
  float a1 = ( B + C / vy ) * a;
  float b1 = ( B + C / vy ) * b + A;
  float c1 = ( B + C / vy ) * c - C / y0 - D;
 
  float b2sub4ac = b1 * b1 - 4 * a1 * c1;
 
  float x1, x2, y1, y2, z1, z2;
  if ( abs( a1 ) > 10E-10 )
  {
    if ( b2sub4ac >= 0 )
    {
      x1 = ( -b1 + sqrt( b2sub4ac ) ) / ( 2 * a1 );
      x2 = ( -b1 - sqrt( b2sub4ac ) ) / ( 2 * a1 );
      y1 = a * x1 * x1 + b * x1 + c;
      y2 = a * x2 * x2 + b * x2 + c;
      z1 = ( y1 - y0 ) / vy;
      z2 = ( y2 - y0 ) / vy;
      // cout<<"in MucGeometron------ x1 y1 z1 = "<<x1<<" "<<y1<<" "<<z1<<endl;
      // cout<<"in MucGeometron------ x2 y2 z2 = "<<x2<<" "<<y2<<" "<<z2<<endl;
      cross1.set( x1, y1, z1 );
      cross2.set( x2, y2, z2 );
      return 1;
    }
    else
    {
      // cout<<"MucGeometron: no intersection!"<<endl;
      cross1.set( -9999, -9999, -9999 );
      cross2.set( -9999, -9999, -9999 );
      return 0;
    }
  }
  else
  {
    x1 = -c1 / b1;
    y1 = a * x1 * x1 + b * x1 + c;
    z1 = ( y1 - y0 ) / vy;
    // cout<<"in MucGeometron------ x1 y1 z1 = "<<x1<<" "<<y1<<" "<<z1<<endl;
    cross1.set( x1, y1, z1 );
    cross2.set( x1, y1, z1 );
    // cout<<"in MucGeometron: abs(a1)<10E-10"<<endl;
    return 1;
  }
}