#include "iostream"
#include "fstream"
#include "cstdio"
#include "TROOT.h"
#include "TFile.h"
#include "TRFIOFile.h"
#include "TTree.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "string.h"
#include "TChain.h"
#include "TLine.h"
#include "TStopwatch.h"
#include "TGraphErrors.h"
#include "TGraph.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TMath.h"

//#ifdef __CINT__
using namespace std; 
#include "corrhm.C"
void corrhp(const Int_t year, const Int_t iweek, string *week, const Int_t iconfig, const Int_t nbinxbj, Double_t *x, Double_t *err_x, Double_t *A1p, Double_t *err_A1p, Double_t *CorrAV, Double_t *err_CorrAV, Int_t *HADP, Int_t *HADN, Int_t *pRun, Int_t *nRun, Int_t &RunMin, Int_t &RunMax);
//#endif 
void corr_() {


  gROOT->Reset();
  gStyle->SetPalette(1,0);
  gStyle->SetOptStat(1111111);
  gStyle->SetOptFit(1);

  /// Inputs
  const int nbinxbj = 12; // number of Xbj ranges
  /// End of inputs

  ofstream xout("x.out");

//   FILE * pFile, * nFile;
//   char pLine [10000], nLine [10000]; 
//   string stp, stn;
  Int_t pRun[100], nRun[100];
  //   size_t posp, posn;
  ifstream in; 

  Int_t RunMin, RunMax;
  Double_t x[nbinxbj], A1p[nbinxbj], CorrAV[nbinxbj], err_x[nbinxbj], err_A1p[nbinxbj], err_CorrAV[nbinxbj];
  Double_t A1av[nbinxbj][2], err_A1av[nbinxbj][2];

  Double_t A1num[nbinxbj][2], A1den[nbinxbj][2];
  Double_t errsum[nbinxbj][2];

  Double_t A1avper[40][nbinxbj][2], err_A1avper[40][nbinxbj][2];
  Double_t A1matrix[40][12][nbinxbj][2], err_A1matrix[40][12][nbinxbj][2];   /// [12] max number of configurations
 
  Int_t HADP[nbinxbj], HADN[nbinxbj], HADPF[nbinxbj], HADNF[nbinxbj];
  for (Int_t i = 0; i < nbinxbj; i++){ // loop to fill A1 matrices
    HADP[i] = 0;
    HADN[i] = 0;
    HADPF[i] = 0;
    HADNF[i] = 0;
  }

  FILE *hadfout;
  hadfout = fopen("hadron_final.out","w"); // output file for configuration i

  Double_t xbinsXBj[101];
  for(Int_t i=0; i<101; i++ ) xbinsXBj[i] = 1e-3 * pow(10,0.03*i);

  for (int year=2004; year<=2004; year++){// loop overall years
    Int_t nbweek = 13;     string week[nbweek]; // start by the largest number of nbweeks. Otherwise will not work. To be clarified, why. Related to Valid index range of array 'a' is 0 to 19. You may wonder why a[20] is not
    //an error.  The reason why is that ANSI-C standard requires C language 
    //processor to allow access of an array element one beyond the boundary.  
    //In this case, a[20] must be allowed. 

    FILE *a1out;
    a1out= fopen(Form("A1_%i.out",year),"w");
    Int_t totalperiods = 0;
    for (int iweek=0; iweek <1; iweek++){ // loop overall weeks
      //first for p
//       gSystem->Exec(Form("awk '(/^%i /) && (/ W%i /) {print substr($0,12,300)}' /comp02/comp/helena/DB/config_p.list > $PWD/my_%i_p%i.txt",year,iweek,year,iweek)); // extracts all first runs of ALL configurations
//       gSystem->Exec(Form("awk '(/^%i /) && (/ W%i /) {print substr($0,6,3)}' /comp02/comp/helena/DB/config_p.list > $PWD/tmpfile",year,iweek)); // extracts week name
//       in.open("tmpfile");
//       in>>week[iweek]; // fills the string array of weeks
//       in.close();
//       //      gSystem->Exec("rm -f tmpfile"); 
//       pFile = fopen (Form("my_%i_p%i.txt",year,iweek), "r");  
//       Int_t p = 0;
//       while (!feof(pFile)) { // in order to reach the end of line
// 	fgetc (pFile);
// 	p++; 
//       }
//       rewind (pFile); // to put the cursor at the beginning of the line
//       stp=fgets (pLine , 10000 , pFile);
//       posp = stp.find(" "); // finds the first space

//       //then for n
//       gSystem->Exec(Form("awk '(/^%i /) && (/ W%i /) {print substr($0,12,300)}' /comp02/comp/helena/DB/config_n.list > $PWD/my_%i_n%i.txt",year,iweek,year,iweek));
//       nFile = fopen (Form("my_%i_n%i.txt",year,iweek), "r");  
//       Int_t n = 0;
//       while (!feof(nFile)) {
// 	fgetc (nFile);
// 	n++;
//       }

//       rewind (nFile);
//       stn=fgets (nLine , 10000 , nFile);
//       posn = stn.find(" ");
      Int_t config = 0;
      for (Int_t iconfig = 0; iconfig<1; iconfig++){
      //      for (Int_t iconfig = 0; iconfig<1; iconfig++){
	pRun[iconfig]=0; nRun[iconfig]=0;
// 	sscanf (stp.substr (posp,6).c_str(), "%d", &pRun[iconfig]); // extracts the first run of the FIRST configuration and fills the array pRun
// 	sscanf (stn.substr (posn,6).c_str(), "%d", &nRun[iconfig]);
  
// 	if ((week[iweek].c_str()=="P1C" && year==2002 && iconfig==0)){ //empty configuration
// 	  posp = posp+6; // advances to the first run of the NEXT configuration 
// 	  posn = posn+6;
// 	  gSystem->Exec("rm -f pFile"); 
// 	  gSystem->Exec("rm -f nFile"); 
// 	  continue;}

	fprintf (a1out, "\n");
	fprintf (a1out, "%i data, Period %s \n",year, week[iweek].c_str());
	fprintf (a1out, "\n");
	fprintf (a1out, "\n");
	fprintf (hadfout, "\n");

	corrhp(year, iweek, week,iconfig, nbinxbj, x, err_x, A1p, err_A1p, CorrAV, err_CorrAV, HADP, HADN, pRun, nRun, RunMin, RunMax); 
 
	fprintf (a1out, "Period %s - configuration %i, %i::%i \n", week[iweek].c_str(),iconfig, RunMin, RunMax); // week[iweek].c_str() transforms the string week[iweek] in an array of characters

	fprintf (a1out, "\n");

	for (Int_t j = 0; j < nbinxbj; j++){ // loop to fill A1 matrices
	  A1matrix[iweek][iconfig][j][0] = A1p[j];	 
	  A1matrix[iweek][iconfig][j][1] = CorrAV[j];
	  err_A1matrix[iweek][iconfig][j][0] = err_A1p[j];	 
	  err_A1matrix[iweek][iconfig][j][1] = err_CorrAV[j];
	  fprintf (a1out, "A1[h+][%s][%i][%i] =\t % 8.6f \t err_A1[h+][%s][%i][%i] = \t % 8.6f A1[h-][%s][%i][%i] =\t % 8.6f \t err_A1[h-][%s][%i][%i] = \t % 8.6f \n",week[iweek].c_str(),iconfig,j,A1matrix[iweek][iconfig][j][0],week[iweek].c_str(),iconfig,j,err_A1matrix[iweek][iconfig][j][0],week[iweek].c_str(),iconfig,j,A1matrix[iweek][iconfig][j][1],week[iweek].c_str(),iconfig,j,err_A1matrix[iweek][iconfig][j][1]);
	  xout<<x[j]<<" "<<err_x[j]<<endl;
	  HADPF[j] += HADP[j];
	  HADNF[j] += HADN[j];
	  fprintf(hadfout, "HADPF[%i] = %i \t HADNF[%i] = %i  \t\t year = %i \t period = %s \t config = %i \n", j, HADPF[j],j, HADNF[j], year, week[iweek].c_str(), iconfig );
	}
	fprintf (a1out, "\n");
	fprintf (a1out, "\n");

// 	posp = posp+6; // advances to the first run of the NEXT configuration 
// 	posn = posn+6;
	config++;

      } ///      for (int iconfig=0; iconfig<n/6;iconfig++){
      gSystem->Exec("rm -f pFile"); 
      gSystem->Exec("rm -f nFile"); 
      
  
      /////////////////////  Calculates the average of the several asymmetries for 1 week /////////////////////


      fprintf (a1out, "\n");
      fprintf (a1out, "%i data, Period %s - Average overall configurations %i \n", year, week[iweek].c_str(), config);
      fprintf (a1out, "\n");

      for (Int_t j = 0; j < nbinxbj; j++){  
	for (Int_t k = 0; k < 2; k++){  
	  A1num[j][k] = 0.0;
	  A1den[j][k] = 0.0;
	  errsum[j][k] = 0.0;
	  A1av[j][k] = 0.0;
	  err_A1av[j][k] = 0.0;
	
	  for (Int_t i = 0; i < config; i++){
	    if (A1matrix[iweek][i][j][k] ==0 || err_A1matrix[iweek][i][j][k] ==0) continue;
	    //	  if (year==2002&&(j==10||j==11)) continue;
	    A1num[j][k] += A1matrix[iweek][i][j][k]*(1/(err_A1matrix[iweek][i][j][k]*err_A1matrix[iweek][i][j][k])); // num = numerador
	    A1den[j][k] += (1/(err_A1matrix[iweek][i][j][k]*err_A1matrix[iweek][i][j][k])); //peso = 1/sigma2 den = denominador

	    errsum[j][k] += 1/(err_A1matrix[iweek][i][j][k]*err_A1matrix[iweek][i][j][k]);
      
	    A1av[j][k] = A1num[j][k] / A1den[j][k];
	    err_A1av[j][k] = sqrt( 1 / errsum[j][k]);
	  }
	  A1avper[iweek][j][k]=A1av[j][k];
	  err_A1avper[iweek][j][k]=err_A1av[j][k];
	}
      }
      for (Int_t j = 0; j < nbinxbj; j++){  
	fprintf (a1out, "A1[%i][h+] = \t\t % 8.6f \t err_A1[%i][h+] = \t\t % 8.6f A1[%i][h-] = \t\t % 8.6f \t err_A1[%i][h-] = \t\t % 8.6f \n",j,A1av[j][0],j,err_A1av[j][0],j,A1av[j][1],j,err_A1av[j][1]);
      }
      fprintf (a1out, "\n");
    
      ///////////////////////////////////////// period plot ////////////////////////////////////////////////////// 
      //create the histogram 

      TLine line;
      line.SetLineStyle(1);
      //      line.SetLineWidth(0.3);
      line.SetLineColor(1);

      Double_t A1avp[nbinxbj], err_A1avp[nbinxbj], A1avn[nbinxbj], err_A1avn[nbinxbj];
      for (Int_t j = 0; j < nbinxbj; j++){ 
	A1avp[j] = A1av[0][j];
	err_A1avp[j] = err_A1av[0][j];
	A1avn[j] = A1av[1][j];
	err_A1avn[j] = err_A1av[1][j];
      }

      TCanvas *c1001 = new TCanvas("c1001","A_{1}",200,10,700,500);
      c1001->Divide(1,2);
      //      c1001->SetGrid();

      // create a 2-d histogram to define the range
      c1001->cd(1);
      TH2F *a1p = new TH2F("a1p",Form("A_{1}^{h-} %s %i",week[iweek].c_str(),year), 100, xbinsXBj, 100, -0.5 , 1 );
      a1p->SetXTitle("x_{Bj}");
      a1p->SetYTitle("A_{1}^{h+}");
      a1p->Draw();
      line.DrawLine(0.001, 0, 1, 0); 
      gStyle->SetOptStat(0);
      gPad->SetLogx();
 
      TGraphErrors *gr = new TGraphErrors(nbinxbj,x,A1avp,err_x,err_A1avp);
      gr->SetMarkerColor(kRed);
      gr->SetMarkerStyle(21);
      gr->SetMarkerSize(0.45);
      gr->Draw("P");
      gPad->SetLogx();
      c1001->Update();
 
      c1001->cd(2);
      TH2F *a1n = new TH2F("a1n",Form("A_{1}^{h-} %s %i",week[iweek].c_str(),year), 100, xbinsXBj, 100, -0.5 , 1 );
      a1n->SetXTitle("x_{Bj}");
      a1n->SetYTitle("A_{1}^{h-}");
      a1n->Draw();
      line.DrawLine(0.001, 0, 1, 0); 
      gStyle->SetOptStat(0);
      gPad->SetLogx();
 
      TGraphErrors *gr1 = new TGraphErrors(nbinxbj,x,A1avn,err_x,err_A1avn);
      gr1->SetMarkerColor(kRed);
      gr1->SetMarkerStyle(21);
      gr1->SetMarkerSize(0.45);
      gr1->Draw("P");
      gPad->SetLogx();
      c1001->Update();
      c1001->Print(Form("%i_%s.ps",year,week[iweek].c_str()));
      a1p->Delete();
      a1n->Delete();

      totalperiods++;
    }//for (iweek=0; iweek < nbweek; iweek++){// loop overall periods

    /////////////////////  Calculates the average of the asymmetries overall periods  /////////////////////
    Double_t A1_year[nbinxbj][2], err_A1_year[nbinxbj][2];
    Double_t A1numper[nbinxbj][2], A1denper[nbinxbj][2];
    Double_t errsumper[nbinxbj][2];
    Double_t A1_year_p[nbinxbj], err_A1_year_p[nbinxbj], A1_year_n[nbinxbj], err_A1_year_n[nbinxbj];

    fprintf (a1out, "%i data - Average overall periods (%i) \n",year, totalperiods);
    fprintf (a1out, "\n");

    for (Int_t j = 0; j < nbinxbj; j++){  
      for (Int_t k = 0; k < 2; k++){  
	A1numper[j][k] = 0.0;
	A1denper[j][k] = 0.0;
	errsumper[j][k] = 0.0;
	A1_year[j][k] = 0.0;
	err_A1_year[j][k] = 0.0;

	for (Int_t i = 0; i < nbweek; i++){ 
	  if (A1avper[i][j][k] ==0 || err_A1avper[i][j][k] ==0) continue;
	  //	  if (year==2002&&(j==10||j==11)) continue;

	  A1numper[j][k] += A1avper[i][j][k]*(1/(err_A1avper[i][j][k]*err_A1avper[i][j][k])); // num = numerador
	  A1denper[j][k] += (1/(err_A1avper[i][j][k]*err_A1avper[i][j][k])); //peso = 1/sigma2 den = denominador

	  errsumper[j][k] += 1/(err_A1avper[i][j][k]*err_A1avper[i][j][k]);

	  A1_year[j][k] = A1numper[j][k] / A1denper[j][k];
	  err_A1_year[j][k] = sqrt( 1 / errsumper[j][k]);
	}
      }
      fprintf (a1out, "A1[h+][%i][%i] =\t\t % 8.6f \t err_A1[h+][%i][%i] =\t % 8.6f A1[h-][%i][%i] =\t\t % 8.6f \t err_A1[h-][%i][%i] =\t % 8.6f \n", year, j, A1_year[j][0], year, j, err_A1_year[j][0], year, j, A1_year[j][1], year, j, err_A1_year[j][1]);  
      A1_year_p[j] = A1_year[j][0];
      err_A1_year_p[j] = err_A1_year[j][0];
      A1_year_n[j] = A1_year[j][1];
      err_A1_year_n[j] = err_A1_year[j][1];
     
    }
    fclose(a1out);
    fclose(hadfout);
    ///////////////////////////////////////// year plot ////////////////////////////////////////////////////// 
    //create the histogram 
    TLine line;
    TCanvas *c1002 = new TCanvas("c1002","A_{1} - %i data",200,10,700,500);

    TH2F *ya1n = new TH2F("ya1n","", 100, xbinsXBj, 100, 0. , 0.65 );
    ya1n->SetXTitle("x_{Bj}");
    ya1n->SetYTitle("corr(A_{1d}^{h+}, A_{1d}}");
    ya1n->Draw();
    line.DrawLine(0.001, 0, 1, 0); 
    gStyle->SetOptStat(0);
    gPad->SetLogx();
    gPad->SetGridy();
    TGraph *lr1 = new TGraphErrors(nbinxbj,x,A1_year_n);
    lr1->SetMarkerColor(kBlue);
    lr1->SetMarkerStyle(21);
    lr1->SetMarkerSize(0.45);
    lr1->Draw("LP");
    gPad->SetLogx();
    TLegend *leg2 = new TLegend(0.1,0.75,0.6,0.89);
    leg2-> AddEntry(lr1,Form("W27 - %i data",year),"p");
    leg2->Draw();
    leg2-> SetFillColor(0);
    c1002->Update();
    //    c1002->Print(Form("A1_%i.ps",year));
    c1002->Print(Form("A1_inc_2004_27_hp.ps",year));
  } //for (yr...// loop overall years
}