#include <TFile.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TString.h>
#include <TMath.h>
#include <TLegend.h>

//_________________
void BuildDependencies() {

  const Char_t *mInFileName[2];
  TFile *mInFile[2];
  mInFileName[0] = "oKaonHBT_auau200_pid2.root";
  mInFileName[1] = "oPionHBT_auau200_pid2.root";
  for(Int_t i=0; i<2; i++) {
    mInFile[i] = TFile::Open(mInFileName[i]);
  }
  
  TFile *mOutFile = new TFile("oHBTdep_auau200_pid2.root","recreate");

  ///////////////////////
  //       Au+Au       //
  ///////////////////////
  TGraphErrors *mMultDepGrErr[2][3];   //Type, Charge
  TGraphErrors *mMultKtGrErr[2][3][7]; //Type, Charge, Cenrality
  TGraphErrors *mMultMtGrErr[2][3][7]; //Type, Charge, Cenrality
  
  //Loop over files and graphs
  for(Int_t iFile=0; iFile<2; iFile++) {
    for(Int_t iCharge=0; iCharge<3; iCharge++) {

      mMultDepGrErr[iFile][iCharge] = (TGraphErrors*)mInFile[iFile]
	->Get(Form("gMultDepRadiiGrErr_%d",iCharge));
      if(iFile == 0) {
	mMultDepGrErr[iFile][iCharge]->SetName(Form("gKaonMultDepRadiiGrErr_%d",iCharge));
      }
      else {
	mMultDepGrErr[iFile][iCharge]->SetName(Form("gPionMultDepRadiiGrErr_%d",iCharge));
      }

      for(Int_t iCent=0; iCent<7; iCent++) {
	mMultKtGrErr[iFile][iCharge][iCent] = (TGraphErrors*)mInFile[iFile]
	  ->Get(Form("gMultKtDepRadiiGrErr_%d_%d",iCharge,(iCent+2)));
	mMultMtGrErr[iFile][iCharge][iCent] = (TGraphErrors*)mInFile[iFile]
	  ->Get(Form("gMultKtDepRadiiMtGrErr_%d_%d",iCharge,(iCent+2)));
	if(iFile == 0) {
	  mMultKtGrErr[iFile][iCharge][iCent]->
	    SetName(Form("gKaonMultKtDepRadiiGrErr_%d_%d",iCharge,(iCent+2)));
	  mMultMtGrErr[iFile][iCharge][iCent]->
	    SetName(Form("gKaonMultMtDepRadiiGrErr_%d_%d",iCharge,(iCent+2)));
	}
	else {
	  mMultKtGrErr[iFile][iCharge][iCent]->
	    SetName(Form("gPionMultKtDepRadiiGrErr_%d_%d",iCharge,(iCent+2)));
	  mMultMtGrErr[iFile][iCharge][iCent]->
	    SetName(Form("gPionMultMtDepRadiiGrErr_%d_%d",iCharge,(iCent+2)));
	}
      } //for(Int_t iCent=2; iCent<9; iCent++)
    } //for(Int_t iCharge=0; iCharge<3; iCharge++)
  } //for(Int_t iFile=0; iFile<2; iFile++)

  TMultiGraph *mgMultDep = new TMultiGraph;
  for(Int_t iFile=0; iFile<2; iFile++) {
    mgMultDep->Add(mMultDepGrErr[iFile][2]);
  }

  TMultiGraph *mgMultKtDep = new TMultiGraph;
  TMultiGraph *mgMultMtDep = new TMultiGraph;
  TMultiGraph *mgMultMtSumDep = new TMultiGraph;
  for(Int_t iFile=0; iFile<2; iFile++) {

    //Separate charges
    for(Int_t iCharge=0; iCharge<2; iCharge++) {
      for(Int_t iCent=0; iCent<7; iCent++) {
	if(iFile == 0) {
	  mgMultKtDep->Add(mMultKtGrErr[iFile][iCharge][iCent]);
	}
	mgMultMtDep->Add(mMultMtGrErr[iFile][iCharge][iCent]);
      }
    }

    //Sum of charges
    for(Int_t iCent=0; iCent<7; iCent++) {
      mgMultMtSumDep->Add(mMultMtGrErr[iFile][2][iCent]);
    }
  }

  Double_t mCentrality[9] = {15,31,58,98,157,234,335,430,482};
  Double_t mCentralityErr[9] = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1};
  Double_t mCentralityRad[9] = {1.78, 2.13, 2.45, 2.77, 3.13, 3.45, 3.78, 3.93, 4.04};
  Double_t mCentralityRadErr[9] = {0.03, 0.04, 0.03, 0.05, 0.04, 0.04, 0.05, 0.06, 0.05};
  TGraphErrors *mCentrSysErrTGr = new TGraphErrors(9, mCentrality, mCentralityRad, 
						   mCentralityErr, mCentralityRadErr);
  mCentrSysErrTGr->SetMarkerStyle(20);
  mCentrSysErrTGr->SetMarkerColor(1);
  mCentrSysErrTGr->SetLineWidth(2);
  //mgMultDep->Add(mCentrSysErrTGr);

  TCanvas *cCentDep = new TCanvas("cCentDep","cCentDep",800,600);
  mgMultDep->Draw("AP");
  mCentrSysErrTGr->Draw("[]");

  TCanvas *cCentKtDep = new TCanvas("cCentKtDep","cCentKtDep",800, 600);
  mgMultKtDep->Draw("AP");

  TCanvas *cCentMtDep = new TCanvas("cCentMtDep","cCentMtDep",800,600);
  mgMultMtDep->Draw("AP");

  TCanvas *cCentMtSumDep = new TCanvas("cCentMtSumDep","cCentMtSumDep",800.600);
  mgMultMtSumDep->Draw("APL");

  /*
  ///////////////////////
  //         pp        //
  ///////////////////////
  Double_t mKaonRadiiKtVal[4] = { 1.24, 1.14, 1.12, 1.16};
  Double_t mKaonRadiiKtValErr[4] = {0.11, 0.10, 0.11, 0.08};
  Double_t mKaonKtVal[4] = {0.21, 0.35, 0.45, 0.66};
  Double_t mKaonKtValErr[4] = {0.01, 0.01, 0.01, 0.01};
  Double_t mKaonMtVal[4];
  Double_t KAON_MASS_SQR = 0.24371698;
  for(Int_t i=0; i<4; i++) {
    mKaonMtVal[i] = TMath::Sqrt(KAON_MASS_SQR + mKaonKtVal[i]);
  } 
  TGraphErrors *mKaonMtRadii = new TGraphErrors(4, mKaonMtVal, mKaonRadiiKtVal, 
						mKaonKtValErr, mKaonRadiiKtValErr);
  mKaonMtRadii->SetName(Form("mKaonMtRadii"));
  Double_t mPionRadiiKtVal[4] = {1.32, 1.26, 1.18, 1.05};
  Double_t mPionRadiiKtValErr[4] = {0.02, 0.02, 0.02, 0.03};
  Double_t mPionKtVal[4] = {0.25, 0.35, 0.45, 0.6};
  Double_t mPionKtValErr[4] = {0.01, 0.01, 0.01, 0.01};
  Double_t mPionMtVal[4];
  Double_t PION_MASS_SQR = 0.0194797;
  for(Int_t i=0; i<4; i++) {
    mPionMtVal[i] = TMath::Sqrt(PION_MASS_SQR + mPionKtVal[i]);
  } 
  TGraphErrors *mPionMtRadii = new TGraphErrors(4, mPionMtVal, mPionRadiiKtVal, 
						mPionKtValErr, mPionRadiiKtValErr);
  mPionMtRadii->SetName(Form("mPionMtRadii"));

  TMultiGraph *mgSmallMtDep = new TMultiGraph;
  mgSmallMtDep->Add(mKaonMtRadii);
  mgSmallMtDep->Add(mPionMtRadii);
  TCanvas *cSmallMtDep = new TCanvas("cSmallMtDep","cSmallMtDep",800,600);
  mgSmallMtDep->Draw("AP");
  */

  cCentDep->Write();
  cCentKtDep->Write();
  cCentMtDep->Write();
  cCentMtSumDep->Draw();

  mOutFile->Write();
  mOutFile->Close();
}