Segment.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
///
/// \file     Segment.cxx
/// \brief    Segments are a basic unit for finding vertices
/// \author   messier@indiana.edu
///
////////////////////////////////////////////////////////////////////////
#include "VertexDT/func/Segment.h"

#include <cmath>
#include <iostream>

namespace vdt
{

  Segment::Segment(bool flip) :
    fFlipped(flip),
    fX1(0),
    fY1(0),
    fX2(0),
    fY2(0),
    fChi2Line0(-1),
    fSw(0),
    fSwx(0),
    fSwx2(0),
    fSwy(0),
    fSwy2(0),
    fSwxy(0)
  { }

  //......................................................................

  void Segment::Flip(bool flip) 
  {
    if (flip!=fFlipped) {
      fFlipped = flip;
      if (fSw!=0) {
        std::cerr << __FILE__ << ":" << __LINE__ 
                  << " Flip called after segment initialization. "
                  << "Will lead to incorrect results!" << std::endl;
        abort();
      }
    }
  }

  //......................................................................

  void Segment::SetEndPoint1(double x, double y, double resin, double resex) 
  {
    if (!fFlipped) { fX1 = x; fY1 = y; }
    else           { fX1 = y; fY1 = x; }
  
    fdX = fX2-fX1;
    fdY = fY2-fY1;

    fResIN1sqr = resin*resin;
    fResEX1    = resex;
  
    double mag = sqrt(fdX*fdX+fdY*fdY);
    fdX /= mag;
    fdY /= mag;
  }

  //......................................................................

  void Segment::SetEndPoint2(double x, double y, double resin, double resex) 
  {
    if (!fFlipped) { fX2 = x; fY2 = y; }
    else           { fX2 = y; fY2 = x; }
  
    fdX = fX2-fX1;
    fdY = fY2-fY1;

    fResIN2sqr = resin*resin;
    fResEX2    = resex;

    double mag = sqrt(fdX*fdX+fdY*fdY);
    fdX /= mag;
    fdY /= mag;
  }

  //......................................................................

  void Segment::AddPoint(double xx, double sigx, double yy, double sigy) 
  {
    double x, y, w;
    if (!fFlipped) { x = xx; y = yy; w = 1.0/(sigy*sigy); }
    else           { x = yy; y = xx; w = 1.0/(sigx*sigx); }
  
    double wx = w*x;
    double wy = w*y;
  
    fSw   += w;
    fSwx  += wx;
    fSwx2 += wx*x;
    fSwy  += wy;
    fSwy2 += wy*y;
    fSwxy += wx*y;

    //
    // Changing the sums invalidates this chi^2 min
    //
    fChi2Line0 = -1;
  }

  //......................................................................

  double Segment::BestM(double x, double y) 
  {
    return ((fSw*y-fSwy)*x-fSwx*y+fSwxy)/((fSw*x-2*fSwx)*x+fSwx2);
  }

  //......................................................................

  double Segment::Chi2(double x, double y)
  {
    return this->Chi2Line(x,y)+this->Chi2EndPoint(x,y);
  }

  //......................................................................

  void Segment::CalcChi2Line0() 
  {
    double d = fSw*fSwx2-fSwx*fSwx;
    if (d==0.0) d = 1.0E-9;

    double m = (fSw*fSwxy-fSwx*fSwy)/d;
    double b = (fSwx2*fSwy-fSwx*fSwxy)/d;

    fChi2Line0 = 
      m*m*fSwx2 + 2*b*m*fSwx - 2*m*fSwxy + b*b*fSw - 2*b*fSwy + fSwy2;
  }

  //......................................................................

  double Segment::Chi2Line(double xx, double yy) 
  {
    double x, y;
    if (!fFlipped) { x = xx; y = yy; }
    else           { x = yy; y = xx; }

    if (fChi2Line0<0) this->CalcChi2Line0();

    double m = this->BestM(x,y);
    double c =
      ( ( (x*fSw - 2*fSwx)*x + fSwx2 )*m + 
        ( (-2*y*fSw + 2*fSwy)*x + 2*y*fSwx - 2*fSwxy ) )*m +
      (y*fSw - 2*fSwy)*y + fSwy2;

    return c-fChi2Line0;
  }

  //......................................................................

  double Segment::Chi2EndPoint(double xx, double yy) 
  {
    double x, y;
    if (!fFlipped) { x = xx; y = yy; }
    else           { x = yy; y = xx; }

    double chi1;
    double d1 = this->GetD1(x,y);
    if (d1>0.0) chi1 = 2.0*d1/fResEX1;
    else        chi1 = d1*d1/fResIN1sqr;

    double chi2;
    double d2 = this->GetD2(x,y);
    if (d2>0.0) chi2 = 2.0*d2/fResEX2;
    else        chi2 = d2*d2/fResIN2sqr;
  
    if (chi1<chi2) return chi1;
    return chi2;
  }

  //......................................................................

  double Segment::GetD1(double x, double y) 
  {
    double dx10 = fX1-x;
    double dy10 = fY1-y;
    return (dx10*fdX)+(dy10*fdY);
  }

  //......................................................................

  double Segment::GetD2(double x, double y) 
  {
    double dx20 = fX2-x;
    double dy20 = fY2-y;
    return -(dx20*fdX)-(dy20*fdY);
  }

} // end namespace vdt

////////////////////////////////////////////////////////////////////////