---

LSQFit.cxx

Go to the documentation of this file.

//$Id: LSQFit.cxx,v 1.2 2000/12/14 23:06:13 burnett Exp $

// Author: Bill Atwood
// Creation: Monday, June  28, 1993

// Contents ----------------------------------------------------------
//
//
// Description
//
//    Implementation of class LSQFit member functions.
//
// End ---------------------------------------------------------------


// Interface Dependencies --------------------------------------------


#include "reconstruction/LSQFit.h"

#include "geometry/Ray.h"


#include <cmath>

inline double sqr(double x){return x*x;}

LSQFit::LSQFit ()
{
}

LSQFit::LSQFit (const LSQFit& fit)
{
        unsigned numData = fit.dataList.size();
        for(unsigned i=0; i<numData; i++) {
          dataList.push_back(FitData( fit.dataList[i]) );
        }
        chisq = fit.chisq;
        rmsResid = fit.rmsResid;
        x0 = fit.x0;
        slopeX = fit.slopeX;
}


LSQFit::~LSQFit ()
{
   clear();
}



void LSQFit::clear ()
{
  //  Remove all  hits.
  dataList.clear();
}

float LSQFit::maxResidual (int* index)
{
  float test, maxRes = 0.;
  int idx = -1;

  for(unsigned i=0; i<dataList.size(); i++) {
     test = dataList[i].res*dataList[i].res;
      if(maxRes < test) {
         idx = (int)dataList[i].index;
         maxRes = test;
      }
  }

   if(idx >= 0) maxRes = sqrt(maxRes);
   *index = idx;
   return maxRes;
  //## end LSQFit::maxResidual%-1294063902.body
}

unsigned int LSQFit::numDataPoints () const
{
  //## begin LSQFit::numDataPoints%770251369.body preserve=yes
    return (int)dataList.size();
}

int LSQFit::compareFits (const LSQFit& fit)
{
  //## begin LSQFit::compareFits%-248894077.body preserve=yes
   int numComData = 0;
   unsigned int fit_numHits = fit.dataList.size();
   unsigned int numHits     = dataList.size(); 
   for(unsigned i=0; i<numHits; i++) {
       for(unsigned j=0; j<fit_numHits; j++) {
          if(dataList[i].index == fit.dataList[j].index) {
             numComData++;
             break;
          }
       }
    }
    return numComData;
  //## end LSQFit::compareFits%-248894077.body
}

Point LSQFit::getHit (unsigned int i)
{
 if(i < dataList.size())  return Point(dataList[i].x, 0, dataList[i].z);
 else                       return Point(FLT_MAX, FLT_MAX, FLT_MAX);
}

unsigned int LSQFit::getHitIndex (unsigned int i)
{
 if(i < dataList.size())  return dataList[i].index;
 else                        return 0;
}

float LSQFit::doFit ()
{
//  Check on number of points  with non-zero weights
 chisq = -.5;    // a small negative number to allow for zero in special cases
 if(dataList.size() < 3) {chisq = -1; rmsResid = -1; return chisq;}

 double s1x,  sx,  szwx, szzwx,  sxx, sxz,  detx,  fit;
 s1x  = sx = szwx = szzwx = sxx = sxz =  0.0;
 unsigned i;

 for(i=0; i<dataList.size(); i++) {
       float x = dataList[i].x;
       float z = dataList[i].z;
       s1x += dataList[i].wt;

       szwx += dataList[i].wt*z;
       szzwx += dataList[i].wt*z*z;

       sx  += dataList[i].wt*x;
       sxx += dataList[i].wt*x*x;
       sxz += dataList[i].wt*z*x;
 }

 detx = s1x*szzwx - szwx*szwx;
 if(detx  == 0) {chisq = -2; rmsResid = 0.; return chisq;}

 x0      = (sx*szzwx - szwx*sxz)/detx;
 slopeX  = (s1x*sxz - szwx*sx) /detx;

 chisq = 0; rmsResid = 0;
  float lst_slope = 0., this_slope, x_lst = 0., z_lst = 0.;
 for(i=0; i<dataList.size(); i++) {

       float x = dataList[i].x;
       float z = dataList[i].z; 
       fit = x0 + slopeX * z;
       dataList[i].res = sqrt(dataList[i].wt) * (x-fit);
       chisq += sqr( dataList[i].res );

       if(i>0) {
           this_slope = fabs(atan((x-x_lst)/(z-z_lst)));
           if(i> 1) rmsResid += fabs(this_slope - lst_slope);
           lst_slope  = this_slope;
       }
       x_lst = x;
       z_lst = z;
  }
  
  chisq =  chisq/(dataList.size() - 2);
  rmsResid = rmsResid/(dataList.size() - 2); 

  return chisq;
  
}

int LSQFit::addData (float x, float z, float wx)
{
  
    unsigned numData = dataList.size();
    dataList.push_back(  FitData( numData, x, z, wx) );
    return dataList.size();
  
}

int LSQFit::addData (unsigned int idx, float x, float z, float wx)
{
  
    dataList.push_back( FitData( idx, x, z, wx) );
    return (int)dataList.size();
  
}

int LSQFit::removeData (unsigned int index)
{
  
    for(unsigned i=0; i <dataList.size(); i++) {
        if ( index == dataList[i].index ) {
            dataList.erase(&dataList[i]);
            break;
        }
     }
    return (int)dataList.size();
  
}

float LSQFit::truncateFit ()
{
  
   if(dataList.size() > 3) {
     int   maxIndex;
     maxResidual(&maxIndex);
     removeData(maxIndex);
     return doFit();
   }
   return chisq;
  
}

void LSQFit::printOn (std::ostream& os) const
{
}

---