#include<iostream>
#include"Riostream.h"
#include"TFile.h"
#include"TH1.h"
#include"TNtuple.h"
#include"TMatrixD.h"
#include"TMatrixT.h"
#include"TStyle.h"
#include"TCanvas.h"

using namespace std;

void FSA2()
{
	TFile *f1 = new TFile("E:/Data/file_root/ketqua/fileketqua.root","recreate");
	Double_t phomau[2048],phochuan[4][2048];
	Double_t anpha[4][4],beta[4],a[4],xichma[4];
	Int_t i;
	const int nlines=2048; // so dong trong file TKA

	//pho betong
	TFile file("E:/Data/file_root/mauthuong/hopdami20_truphong.root","h3");
	TH1F *hpho = (TH1F*)file.Get("pho");
	for (i=0;i<nlines;i++)
	{
		Float_t value = h3->GetBinContent(i);
		phomau[i]=value;
	}
	//file.close();

// write phochuan

//cat
    TFile file1("E:/Data/file_root/mauthanhphan/hopcat20_truphong.root","h3");
	TH1F *hpho1 = (TH1F*)file1.Get("pho1");
	for (i=0;i<nlines;i++)
	{
		Float_t value1 = h3->GetBinContent(i);
		phochuan[0][i]=value1;
	}
	//file1.close();

	//da
	TFile file2("E:/Data/file_root/mauthanhphan/hopdami20_truphong.root","h3");
	TH1F *hpho2 = (TH1F*)file2.Get("pho2");
	for (i=0;i<nlines;i++)
	{
		Float_t value2 = h3->GetBinContent(i);
		phochuan[1][i]=value2;
	}
	//file2.close();

	//nuoc
	TFile file3("E:/Data/file_root/mauthanhphan/hopnuoc20_truphong.root","h3");
	TH1F *hpho3 = (TH1F*)file3.Get("pho3");
	for (i=0;i<nlines;i++)
	{
		Float_t value3 = h3->GetBinContent(i);
		phochuan[2][i]=value3;
	}
	//file3.close();

	//xi mang
	TFile file4("E:/Data/file_root/mauthanhphan/hopximang20_truphong.root","h3");
	TH1F *hpho4 = (TH1F*)file4.Get("pho4");
	for (i = 0; i < nlines; i++)
	{
		Float_t value4 = h3->GetBinContent(i);
		phochuan[3][i]=value4;
	}
	//file4.close();

	//Tinh ma tran
	Int_t j,x;
	for (j = 0; j < 4; j++)
	{
		beta[j]=0;
		a[j]=0;
		for (x=0;x<4;x++)
			anpha[j][x]=0;
	}
	/*for (j=0;j<4;j++){
		beta[j]=0;
		a[j]=0;
		for (x=0;x<4,x++)
		{
			anpha[j][x]=0;
		}
	}*/
	Int_t c;
	for (c=0;c<nlines;c++)
	{
		beta[0] += phomau[c]*phochuan[0][c];
		beta[1] += phomau[c]*phochuan[1][c];
		beta[2] += phomau[c]*phochuan[2][c];
		beta[3] += phomau[c]*phochuan[3][c];
	}
	Int_t m,n,b;
	for (m=0; m < 4; m++)
	{
		for (n=0; n < 4; n++)
		{
			for (b=0; b < nlines; b++)
			{
				anpha[m][n] += phochuan[m][b]*phochuan[n][b];
	}}}
	// Tim ma tran nghich dao cua anpha la epsilon
	/*Int_t a,q;
	for (a=0;a<4;a++)
	{
		for (q=0;q<4;q++)
			epsilon = anpha.Invert();
	}*/
	//TMatrixD epsilon(kInvert,anpha);
	TMatrixD epsilon = anpha.InvertFast();
	//Tinh ket qua ma tran a
	Int_t r,t;
	for (r=0;r<4;r++)
	{
		for (t=0;t<4;t++)
			a[r] += beta[r]*epsilon[r][t];
	}
	//Tìm sai so cua a
	Int_t z;
	for (z=0;z<4;z++)
	{
		xichma[z] = sqrt(epsilon[z][z]);
	}
	f1->Write();
    f1->Close();
    f1->cd();
	cout<<a =<<"\t"<<a[0]<<"\t"<<a[1]<<"\t"<<a[2]<<"\t"<<a[3]<<endl;
	cout<<ss a =<<"\t"<<xichma[0]<<"\t"<<xichma[1]<<"\t"<<xichma[2]<<"\t"<<xichma[3]<<endl;
}