#include <stdio.h>
#include <iostream>
#include <fstream>
#include <iomanip>

#include <TMatrixDSym.h>
#include <TVectorD.h>


#include <TH2.h>

void RetrPionFFData();

const Int_t nel = 195;

const Int_t nk08 = 60;

const Int_t nk10 = 75;

const Int_t nk12 = 60;

TMatrixDSym K08Acc(nk08), K08FSR(nk08), K08Lmnsty(nk08), K08Trg(nk08);
TMatrixDSym K08Unf(nk08), K08radH(nk08), K08Bkgd(nk08), K08L3(nk08);
TMatrixDSym K08Mtrk(nk08), K08Trk(nk08), K08sqrtsH(nk08);

TMatrixDSym K10Acc(nk10), K10L3(nk10), K10Mtrk(nk10), K10RcnFil(nk10);
TMatrixDSym K10Trk(nk10), K10pffFSR(nk10), K10radH(nk10);
TMatrixDSym K10Bkgd(nk10), K10Lmnsty(nk10), K10OmCut(nk10), K10Trg(nk10);
TMatrixDSym K10Unf(nk10), K10mdFSR(nk10), K10pieID(nk10), K10rhopiBg(nk10);

TMatrixDSym K12AccRatio(nk12), K12FSR(nk12), K12MTrkTail(nk12);
TMatrixDSym K12TrgRatio(nk12), K12Unf(nk12), K12BgRatio(nk12);
TMatrixDSym K12L3Ratio(nk12), K12MtrkRatio(nk12), K12TrkRatio(nk12), K12WghtdBg(nk12), K12VP(nk12);

TMatrixD K0810MTrk(nk08,nk10),K0810Unf(nk08,nk10),K0810Trk(nk08,nk10);
TMatrixD K0810Acc(nk08,nk10),K0810L3(nk08,nk10),K0810Lmnsty(nk08,nk10);
TMatrixD K0810mdFSR(nk08,nk10),K0810radH(nk08,nk10);

TMatrixD K0812Bkgd(nk08,nk12),K0812MTrk(nk08,nk12),K0812Trk(nk08,nk12);
TMatrixD K0812Trg(nk08,nk12),K0812Unf(nk08,nk12),K0812Acc(nk08,nk12);
TMatrixD K0812L3(nk08,nk12),K0812FSR(nk08,nk12);

TMatrixD K1012MTrk(nk10,nk12),K1012Trk(nk10,nk12),K1012Unf(nk10,nk12);
TMatrixD K1012Acc(nk10,nk12),K1012L3(nk10,nk12),K1012FSR(nk10,nk12);

TMatrixDSym Msyst(195);

TH2D *hSyst = new TH2D("hSyst","KLOE systematic covariance matrix",195,0,195,195,0,195);

void PionFF() {

  RetrPionFFData();

  //Outfile:
  ofstream systfile("KLOE_syst_pionFF.dat");
  
    //Start creation of covariance matrix
  for (UInt_t i=0; i<195; i++)
  {
    for (UInt_t j=0; j<195; j++)
      {
//==== "diagonal" contributions of the 3 covariance matrices:
	if (i < 60 && j < 60 )
	  {
	    // syst cov KLOE08:
	    Msyst(i,j)  = K08Acc(i,j); //Acceptance
	    Msyst(i,j) += K08FSR(i,j); //FSR
	    Msyst(i,j) += K08Lmnsty(i,j); //Luminosity
	    Msyst(i,j) += K08Trg(i,j); //Trigger
	    Msyst(i,j) += K08Unf(i,j); //Unfolding
	    Msyst(i,j) += K08radH(i,j); //radiator funtion
	    Msyst(i,j) += K08Bkgd(i,j); //Background
	    Msyst(i,j) += K08L3(i,j); //Level 3 Trigger
	    Msyst(i,j) += K08Mtrk(i,j); //Trackmass
	    Msyst(i,j) += K08Trk(i,j); //Tracking
	    Msyst(i,j) += K08sqrtsH(i,j); //sqrts dependence of radiator
	    //	    Msyst(i,j) += K08VP(i,j); //Vacuum Polarization
	    
	  }
	else if ((i >= 60 && i < 135) && (j >= 60 && j < 135))
	  {
	    // syst cov KLOE10:
	    Msyst(i,j)  = K10Acc(i-60,j-60); //Acceptance
	    Msyst(i,j) += K10L3(i-60,j-60); //Level 3 Trigger
	    Msyst(i,j) += K10Mtrk(i-60,j-60); //Trackmass
	    Msyst(i,j) += K10RcnFil(i-60,j-60); //Reconstruction Filter
	    Msyst(i,j) += K10Trk(i-60,j-60); //Tracking
	    Msyst(i,j) += K10pffFSR(i-60,j-60); //FSR dep. F_pi(1GeV)
	    Msyst(i,j) += K10radH(i-60,j-60); //Radiator Function
	    Msyst(i,j) += K10Bkgd(i-60,j-60); //Background
	    Msyst(i,j) += K10Lmnsty(i-60,j-60); //Luminosity
	    Msyst(i,j) += K10OmCut(i-60,j-60); //Omega cut
	    Msyst(i,j) += K10Trg(i-60,j-60); //Trigger
	    Msyst(i,j) += K10Unf(i-60,j-60); //Unfolding
	    Msyst(i,j) += K10mdFSR(i-60,j-60); //FSR model dep.
	    Msyst(i,j) += K10pieID(i-60,j-60); //pi-e separation
	    Msyst(i,j) += K10rhopiBg(i-60,j-60); //Background rho-pi
	    //	    Msyst(i,j) += K10VP(i-60,j-60); //Vacuum Polarization
	    
	  }
	else if ((i >= 135 && i < 195) && (j >= 135 && j < 195)) 
	  {
	    // syst cov KLOE12:
	    Msyst(i,j)  = K12AccRatio(i-135,j-135); //Acceptance
	    Msyst(i,j) += K12FSR(i-135,j-135); //FSR
	    Msyst(i,j) += K12MTrkTail(i-135,j-135); //Trackmass Tail
	    Msyst(i,j) += K12TrgRatio(i-135,j-135); //Trigger
	    Msyst(i,j) += K12Unf(i-135,j-135); //Unfolding
	    Msyst(i,j) += K12BgRatio(i-135,j-135); //Background
	    Msyst(i,j) += K12L3Ratio(i-135,j-135); //Level 3 Trigger
	    Msyst(i,j) += K12MtrkRatio(i-135,j-135); //Trackmass
	    Msyst(i,j) += K12TrkRatio(i-135,j-135); //Tracking
	    Msyst(i,j) += K12WghtdBg(i-135,j-135); //Background weights
	    Msyst(i,j) += K12VP(i-135,j-135); //Vacuum Polarization	    
	  }
//==== cross correlations start:
//====KLOE08 vs KLOE10:
        if (i < 60 && (j >= 60 && j <  135))
	  {
	    Msyst(i,j)  = K0810MTrk(i,j-60); //Trackmass
	    Msyst(i,j) += K0810Unf(i,j-60); //Unfolding
	    Msyst(i,j) += K0810Trk(i,j-60); //Tracking
	    Msyst(i,j) += K0810Acc(i,j-60); //Acceptance
	    Msyst(i,j) += K0810L3(i,j-60); //Level 3 Trigger
	    Msyst(i,j) += K0810Lmnsty(i,j-60); //Lumminosity
	    Msyst(i,j) += K0810mdFSR(i,j-60); //FSR model dep.
	    Msyst(i,j) += K0810radH(i,j-60); //Radiator Function
	    //	    Msyst(i,j) += K0810VP(i,j-60); //Vacuum Polarization
	    
	    Msyst(j,i) = Msyst(i,j);
	  }
//====KLOE08 vs KLOE12:
        if ( i < 60 && j >= 135)	
	  {
	    Msyst(i,j)  = K0812Bkgd(i,j-135); //Background
	    Msyst(i,j) += K0812MTrk(i,j-135); //Trackmass
	    Msyst(i,j) += K0812Trk(i,j-135); //Tracking
	    Msyst(i,j) += K0812Trg(i,j-135); //Trigger
	    Msyst(i,j) += K0812Unf(i,j-135); //Unfolding
	    Msyst(i,j) += K0812Acc(i,j-135); //Acceptance
	    Msyst(i,j) += K0812L3(i,j-135); //Level 3 Trigger
	    Msyst(i,j) += K0812FSR(i,j-135); //FSR

	    Msyst(j,i) = Msyst(i,j);
	  }
//====KLOE10 vs KLOE12:
	if (((i >= 60 && i < 135)) && (j >= 135 && j < 195))
	  {
	    Msyst(i,j)  = K1012MTrk(i-60,j-135); //Trackmass
	    Msyst(i,j) += K1012Trk(i-60,j-135); //Tracking
	    Msyst(i,j) += K1012Unf(i-60,j-135); //Unfolding
	    Msyst(i,j) += K1012Acc(i-60,j-135); //Acceptance
	    Msyst(i,j) += K1012L3(i-60,j-135); //Level 3 Trigger
	    Msyst(i,j) += K1012FSR(i-60,j-135); //FSR

	    Msyst(j,i) = Msyst(i,j); 
	  }
	
            hSyst->SetBinContent(i+1,j+1, Msyst(i,j));
	    systfile << setw(20) << setiosflags(ios::scientific) << setprecision(16) << Msyst(i,j) << endl;	
      }
  }

  systfile.close();

}

void RetrPionFFData()
{
  Int_t nn=0;
  Double_t val;

  ifstream ifK08Acc("KLOE08/K08Acc_FF.mat");
  nn=0;
    while(ifK08Acc>>val){K08Acc(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08Acc.close();
  
  ifstream ifK08FSR("KLOE08/K08FSR_FF.mat");
  nn=0;
    while(ifK08FSR>>val){K08FSR(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08FSR.close();
  
  ifstream ifK08Lmnsty("KLOE08/K08Lmnsty_FF.mat");
  nn=0;
    while(ifK08Lmnsty>>val){K08Lmnsty(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08Lmnsty.close();
  
  ifstream ifK08Trg("KLOE08/K08Trg_FF.mat");
  nn=0;
    while(ifK08Trg>>val){K08Trg(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08Trg.close();
  
  ifstream ifK08Unf("KLOE08/K08Unf_FF.mat");
  nn=0;
    while(ifK08Unf>>val){K08Unf(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08Unf.close();
  
  ifstream ifK08radH("KLOE08/K08radH_FF.mat");
  nn=0;
    while(ifK08radH>>val){K08radH(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08radH.close();
  
  ifstream ifK08Bkgd("KLOE08/K08Bkgd_FF.mat");
  nn=0;
    while(ifK08Bkgd>>val){K08Bkgd(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08Bkgd.close();
  
  ifstream ifK08L3("KLOE08/K08L3_FF.mat");  
  nn=0;
    while(ifK08L3>>val){K08L3(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08L3.close();
  
  ifstream ifK08Mtrk("KLOE08/K08Mtrk_FF.mat");
  nn=0;
    while(ifK08Mtrk>>val){K08Mtrk(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08Mtrk.close();
  
  ifstream ifK08Trk("KLOE08/K08Trk_FF.mat");
  nn=0;
    while(ifK08Trk>>val){K08Trk(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08Trk.close();
  
  ifstream ifK08sqrtsH("KLOE08/K08sqrtsH_FF.mat");
  nn=0;
    while(ifK08sqrtsH>>val){K08sqrtsH(nn%nk08,nn/nk08)=val; nn++;}  
  ifK08sqrtsH.close();
  
  ifstream ifK10Acc("KLOE10/K10Acc_FF.mat");
  nn=0;
    while(ifK10Acc>>val){K10Acc(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10Acc.close();
  
  ifstream ifK10L3("KLOE10/K10L3_FF.mat");
  nn=0;
    while(ifK10L3>>val){K10L3(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10L3.close();
  
  ifstream ifK10Mtrk("KLOE10/K10Mtrk_FF.mat");
  nn=0;
    while(ifK10Mtrk>>val){K10Mtrk(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10Mtrk.close();
  
  ifstream ifK10RcnFil("KLOE10/K10RcnFil_FF.mat");
  nn=0;
    while(ifK10RcnFil>>val){K10RcnFil(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10RcnFil.close();
  
  ifstream ifK10Trk("KLOE10/K10Trk_FF.mat");
  nn=0;
    while(ifK10Trk>>val){K10Trk(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10Trk.close();
  
  ifstream ifK10pffFSR("KLOE10/K10pffFSR_FF.mat");
  nn=0;
    while(ifK10pffFSR>>val){K10pffFSR(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10pffFSR.close();
  
  ifstream ifK10radH("KLOE10/K10radH_FF.mat");
  nn=0;
    while(ifK10radH>>val){K10radH(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10radH.close();
  
  ifstream ifK10Bkgd("KLOE10/K10Bkgd_FF.mat");
  nn=0;
    while(ifK10Bkgd>>val){K10Bkgd(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10Bkgd.close();
  
  ifstream ifK10Lmnsty("KLOE10/K10Lmnsty_FF.mat");
  nn=0;
    while(ifK10Lmnsty>>val){K10Lmnsty(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10Lmnsty.close();
  
  ifstream ifK10OmCut("KLOE10/K10OmCut_FF.mat");
  nn=0;
    while(ifK10OmCut>>val){K10OmCut(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10OmCut.close();
  
  ifstream ifK10Trg("KLOE10/K10Trg_FF.mat");
  nn=0;
    while(ifK10Trg>>val){K10Trg(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10Trg.close();
  
  ifstream ifK10Unf("KLOE10/K10Unf_FF.mat");
  nn=0;
    while(ifK10Unf>>val){K10Unf(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10Unf.close();
  
  ifstream ifK10mdFSR("KLOE10/K10mdFSR_FF.mat");
  nn=0;
    while(ifK10mdFSR>>val){K10mdFSR(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10mdFSR.close();
  
  ifstream ifK10pieID("KLOE10/K10pieID_FF.mat");
  nn=0;
    while(ifK10pieID>>val){K10pieID(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10pieID.close();
  
  ifstream ifK10rhopiBg("KLOE10/K10rhopiBg_FF.mat");
  nn=0;
    while(ifK10rhopiBg>>val){K10rhopiBg(nn%nk10,nn/nk10)=val; nn++;}  
  ifK10rhopiBg.close();
  
  ifstream ifK12AccRatio("KLOE12/K12AccRatio_FF.mat");
  nn=0;
    while(ifK12AccRatio>>val){K12AccRatio(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12AccRatio.close();
  
  ifstream ifK12FSR("KLOE12/K12FSR_FF.mat");
  nn=0;
    while(ifK12FSR>>val){K12FSR(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12FSR.close();
  
  ifstream ifK12MTrkTail("KLOE12/K12MTrkTail_FF.mat");
  nn=0;
    while(ifK12MTrkTail>>val){K12MTrkTail(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12MTrkTail.close();
  
  ifstream ifK12TrgRatio("KLOE12/K12TrgRatio_FF.mat");
  nn=0;
    while(ifK12TrgRatio>>val){K12TrgRatio(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12TrgRatio.close();
  
  ifstream ifK12Unf("KLOE12/K12Unf_FF.mat");
  nn=0;
    while(ifK12Unf>>val){K12Unf(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12Unf.close();
  
  ifstream ifK12BgRatio("KLOE12/K12BgRatio_FF.mat");
  nn=0;
    while(ifK12BgRatio>>val){K12BgRatio(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12BgRatio.close();
  
  ifstream ifK12L3Ratio("KLOE12/K12L3Ratio_FF.mat");
  nn=0;
    while(ifK12L3Ratio>>val){K12L3Ratio(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12L3Ratio.close();
  
  ifstream ifK12MtrkRatio("KLOE12/K12MtrkRatio_FF.mat");
  nn=0;
    while(ifK12MtrkRatio>>val){K12MtrkRatio(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12MtrkRatio.close();  
  ifstream ifK12TrkRatio("KLOE12/K12TrkRatio_FF.mat");
  nn=0;
    while(ifK12TrkRatio>>val){K12TrkRatio(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12TrkRatio.close();  
  ifstream ifK12WghtdBg("KLOE12/K12WghtdBg_FF.mat");
  nn=0;
    while(ifK12WghtdBg>>val){K12WghtdBg(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12WghtdBg.close();  

  ifstream ifK12VP("KLOE12/K12VP_FF.mat");
  nn=0;
    while(ifK12VP>>val){K12VP(nn%nk12,nn/nk12)=val; nn++;}  
  ifK12VP.close();  
  
  ifstream ifK0810MTrk("KLOE0810/K0810MTrk_FF.mat");
  ifstream ifK0810Unf("KLOE0810/K0810Unf_FF.mat");
  ifstream ifK0810Trk("KLOE0810/K0810Trk_FF.mat");
  ifstream ifK0810Acc("KLOE0810/K0810Acc_FF.mat");
  ifstream ifK0810L3("KLOE0810/K0810L3_FF.mat");
  ifstream ifK0810Lmnsty("KLOE0810/K0810Lmnsty_FF.mat");
  ifstream ifK0810mdFSR("KLOE0810/K0810FSR_FF.mat");  
  ifstream ifK0810radH("KLOE0810/K0810radH_FF.mat"); 

  for (UInt_t j=0; j<nk08; j++)
  {
    for (UInt_t k=0; k<nk10; k++)
      {
	ifK0810MTrk >> K0810MTrk(j,k);
	ifK0810Unf >> K0810Unf(j,k);	
       	ifK0810Trk >> K0810Trk(j,k);
	ifK0810Acc >> K0810Acc(j,k);
	ifK0810L3 >> K0810L3(j,k);
	ifK0810Lmnsty >> K0810Lmnsty(j,k);
	ifK0810mdFSR >> K0810mdFSR(j,k);
	ifK0810radH >> K0810radH(j,k);
      }
  }

  ifK0810MTrk.close();
  ifK0810Unf.close();
  ifK0810Trk.close();
  ifK0810Acc.close();
  ifK0810L3.close();
  ifK0810Lmnsty.close();
  ifK0810mdFSR.close();
  ifK0810radH.close();
  
  ifstream ifK0812Bkgd("KLOE0812/K0812Bkgd_FF.mat");  
  ifstream ifK0812MTrk("KLOE0812/K0812MTrk_FF.mat"); 
  ifstream ifK0812Trk("KLOE0812/K0812Trk_FF.mat"); 
  ifstream ifK0812Trg("KLOE0812/K0812Trg_FF.mat"); 
  ifstream ifK0812Unf("KLOE0812/K0812Unf_FF.mat"); 
  ifstream ifK0812Acc("KLOE0812/K0812Acc_FF.mat"); 
  ifstream ifK0812L3("KLOE0812/K0812L3_FF.mat"); 
  ifstream ifK0812FSR("KLOE0812/K0812FSR_FF.mat"); 

  //KLOE08 - KLOE12 correlation matrices:  
  for (UInt_t j=0; j<nk08; j++)
  {
    for (UInt_t k=0; k<nk12; k++)
      {
	ifK0812Bkgd >> K0812Bkgd(j,k);
	ifK0812MTrk >> K0812MTrk(j,k);
	ifK0812Trk >> K0812Trk(j,k);
	ifK0812Trg >> K0812Trg(j,k);
	ifK0812Unf >> K0812Unf(j,k);
	ifK0812Acc >> K0812Acc(j,k);
	ifK0812L3 >> K0812L3(j,k);
	ifK0812FSR >> K0812FSR(j,k);
      }
  }
  
  ifK0812Bkgd.close();
  ifK0812MTrk.close();
  ifK0812Trk.close();
  ifK0812Trg.close();
  ifK0812Unf.close();
  ifK0812Acc.close();
  ifK0812L3.close();
  ifK0812FSR.close();
  
  ifstream ifK1012MTrk("KLOE1012/K1012MTrk_FF.mat");    
  ifstream ifK1012Trk("KLOE1012/K1012Trk_FF.mat");    
  ifstream ifK1012Unf("KLOE1012/K1012Unf_FF.mat");    
  ifstream ifK1012Acc("KLOE1012/K1012Acc_FF.mat");    
  ifstream ifK1012L3("KLOE1012/K1012L3_FF.mat");    
  ifstream ifK1012FSR("KLOE1012/K1012FSR_FF.mat");    

  //KLOE10 - KLOE12 correlation matrices:  
  for (UInt_t j=0; j<nk10; j++)
  {
    for (UInt_t k=0; k<nk12; k++)
      {
	ifK1012MTrk >> K1012MTrk(j,k);
	ifK1012Trk >> K1012Trk(j,k);
	ifK1012Unf >> K1012Unf(j,k);
	ifK1012Acc >> K1012Acc(j,k);
	ifK1012L3 >> K1012L3(j,k);
	ifK1012FSR >> K1012FSR(j,k);
      }
  }
  
  ifK1012MTrk.close();   
  ifK1012Trk.close();  
  ifK1012Unf.close();  
  ifK1012Acc.close();  
  ifK1012L3.close();   
  ifK1012FSR.close();  

}