using System;
using System.Collections.Generic;
using System.Text;

namespace BrazilTimestampReconstruction
{
    public class LlseFit
    {
        double a, b, siga, sigb, chi2, q;       // b: Alpha, a: Beta; y = a + bx
        Int16 numPoints = -1;
        Int16 dof = -1;
        public LlseFit()
        {
            a = -1;
            b = -1;
        }
        public LlseFit(double alpha, double beta)
        {
            this.b = alpha;
            this.a = beta;
        }
        public Int16 NumPoints
        {
            set { numPoints = value; }
            get { return numPoints; }
        }
        public Int16 DegreesOfFreedom
        {
            set { dof = value; }
            get { return dof; }
        }
        public double Alpha
		{
			get{return b;}
		}
		public double Beta
		{
			get{return a;}
		}
		public double UncertaintyOfA
		{
			get{return siga;}
		}
		public double UncertaintyOfB
		{
			get{return sigb;}
		}
		public double ChiSquare
		{
            set { chi2 = value; }
			get{return chi2;}
		}
		public double GoodnessOfFitProbability
		{
			get{return q;}
		}
		public void FitData(double[] x, double[] y, int ndata, double[] sig, int mwt)
		{
			Gammq gam = new Gammq();
	        int i;
	        double wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat;
            numPoints = (Int16) ndata;
            dof = (Int16) (numPoints - 1);

	        b=0.0;
	        if (mwt!=0) 
			{
		       ss=0.0;
		       for (i=0;i<ndata;i++) 
			   {
			       wt=1.0/sig[i]/sig[i];
			       ss += wt;
			       sx += x[i]*wt;
			       sy += y[i]*wt;
		       }
	        } 
			else 
			{
		       for (i=0;i<ndata;i++) 
			   {
			       sx += x[i];
			       sy += y[i];
		       }
		       ss=ndata;
	        }
	        sxoss=sx/ss;
	        if (mwt!=0) 
			{
		       for (i=0;i<ndata;i++) 
			   {
			      t=(x[i]-sxoss)/sig[i];
			      st2 += t*t;
			      b += t*y[i]/sig[i];
               }
	        } 
			else 
			{
		       for (i=0;i<ndata;i++) 
			   {
			      t=x[i]-sxoss;
			      st2 += t*t;
			      b += t*y[i];
		       }
	        }
	        b /= st2;
	        a=(sy-sx*b)/ss;
	        siga=Math.Sqrt((1.0+sx*sx/(ss*st2))/ss);
	        sigb=Math.Sqrt(1.0/st2);
	        chi2=0.0;
	        if (mwt == 0) 
			{
		       for (i=0;i<ndata;i++)
		          chi2 += (y[i]-a-b*x[i])*(y[i]-a-b*x[i]);
		       q=1.0;
		       sigdat=Math.Sqrt(chi2/(ndata-2));
		       siga *= sigdat;
		       sigb *= sigdat;
	        } 
			else 
			{
		       for (i=0;i<ndata;i++)
			      chi2 += ((y[i]-a-b*x[i])/sig[i])*((y[i]-a-b*x[i])/sig[i]);
		       q=gam.gammq(0.5*(ndata-2),0.5*chi2);
	        }

            // Test this part
            // compute reduced chi square
            chi2 /= numPoints;
        }
    }
}