#include <MdcLine.h>

Public Member Functions
	MdcLine (int n)

	~MdcLine ()

int	fit (int nUse)

Public Attributes
double *	x

double *	y

double *	sigma

double *	resid

int	nPoint

double	slope

double	intercept

double	chisq

double	errmat [3]

Detailed Description

Definition at line 13 of file MdcLine.h.

Constructor & Destructor Documentation

◆ MdcLine()

MdcLine::MdcLine ( int n )

inline

Definition at line 25 of file MdcLine.h.

                 { 
    nPoint = n;
    x = new double[n];
    y = new double[n];
    sigma = new double[n];
    resid = new double[n];
  };

◆ ~MdcLine()

MdcLine::~MdcLine ( )

inline

Definition at line 32 of file MdcLine.h.

             { 
    delete [] x;
    delete [] y;
    delete [] sigma;
    delete [] resid;
  };

Member Function Documentation

◆ fit()

int MdcLine::fit ( int nUse )

Definition at line 10 of file MdcLine.cxx.

                         {
 
  double s1, sx, sy, sxx, sxy;  // The various sums. 
  double delinv, denom;
  double weight = 0.013;     // weight for hits, calculated from sigmas. 
  int ihit;
  if (nUse > nPoint) nUse = nPoint;
  s1 = sx = sy = sxx = sxy = 0.0;
  chisq = 0.0;
  for (ihit = 0; ihit < nUse; ihit++) {
    
    if (sigma[ihit] < 0.0) continue;
    weight = 1. / (sigma[ihit] * sigma[ihit]);//NEED sigma of MdcHits
    s1 += weight;    
    sx += x[ihit] * weight;
    sy += y[ihit] * weight;
    sxx += x[ihit] * (x[ihit] * weight);
    sxy += y[ihit] * (x[ihit] * weight);
 
    //std::cout<<__FILE__<<" "<<__LINE__<< " "<<ihit<<" sigma "<<
      //sigma[ihit] <<" weight "
      //<< weight<<" x "
      //<< x[ihit]<<" y "
      //<< y[ihit]<< std::endl;
  }
  
  /* Calculate parameters. */
  denom = s1 * sxx - sx * sx;
  delinv = (denom == 0.0) ? 1.e20 : 1. / denom;
  intercept = (sy * sxx - sx * sxy) * delinv;
  slope = (s1 * sxy - sx * sy) * delinv;
  errmat[0] = sxx * delinv;
  errmat[1] = -sx * delinv;
  errmat[2] = s1 * delinv;
  //std::cout << " After Fit:" << std::endl;//yzhang debug
  /* Calculate residuals. */
  for (ihit = 0; ihit < nUse; ihit++) {
    if (sigma[ihit] < 0.0) continue;
    resid[ihit] = ( y[ihit] - intercept - slope * x[ihit] );
    chisq += resid[ihit] * resid[ihit] / (sigma[ihit] * sigma[ihit]);
    //std::cout << "resid "<<ihit<<" "<<resid[ihit]<<" sigma "<<sigma[ihit]<<" chisq add "<<resid[ihit] * resid[ihit] / (sigma[ihit] * sigma[ihit]) << std::endl;//yzhang debug 
  }
 
  //std::cout<<__FILE__<<" "<<__LINE__
  //<< " intercept "<<intercept<<" slope "<<slope 
      //<< " chisq "<<chisq<<std::endl;//yzhang debug 
  return 0;
}