#include <TFile.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TH3.h>
#include <iostream>
#include <vector>
#include <TString.h>
#include "HbtFitter.h"

//_________________
void Test3DFitter(const char* inFileName = "auau200_pions_3d_cent8.root",
		  const char* oFileName = "oTest3DFitter.root") {

  Int_t mPartType = 3;
  Int_t mSourceForm = 0;
  Int_t mNonFemtoForm = 0;
  Double_t mCoulombSR = 6.;
  Bool_t mExcludePoints = false;
  Double_t mProjValLo = 0.;
  Double_t mProjValHi = 0.05;

  TFile *inFile = TFile::Open(inFileName);
  //TH3F* hNumer = (TH3F*)inFile->Get("hQinvNum");
  //TH3F* hDenom = (TH3F*)inFile->Get("hQinvDen");
  TH3F* hNumer = (TH3F*)inFile->Get("hPionBPLCMSFrameKtCF_1_Num_bin_2");
  TH3F* hDenom = (TH3F*)inFile->Get("hPionBPLCMSFrameKtCF_1_Den_bin_2");
  TH3F* hQinv = (TH3F*)inFile->Get("hPionBPLCMSFrameKtCF_1_Qinv_bin_2");
  //Normalazing Qinv histogram in order to get pure Qinv weights
  hQinv->Divide(hDenom);

  Int_t mQoutNbins = hNumer->GetNbinsX();
  Double_t mQoutLo = hNumer->GetXaxis()->GetBinLowEdge(1);
  Double_t mQoutHi = hNumer->GetXaxis()->GetBinUpEdge(mQoutNbins);
  Int_t mQsideNbins = hNumer->GetNbinsY();
  Double_t mQsideLo = hNumer->GetYaxis()->GetBinLowEdge(1);
  Double_t mQsideHi = hNumer->GetYaxis()->GetBinUpEdge(mQsideNbins);
  Int_t mQlongNbins = hNumer->GetNbinsZ();
  Double_t mQlongLo = hNumer->GetZaxis()->GetBinLowEdge(1);
  Double_t mQlongHi = hNumer->GetZaxis()->GetBinUpEdge(mQlongNbins);
  HbtFitter *fitter = new HbtFitter(mPartType,mSourceForm,mNonFemtoForm,mExcludePoints);
  fitter->SetFitMethod(1);          //Log-Likelihood
  fitter->SetCorrectOnMc(false);
  fitter->SetProjectionType(false); //Project by vals
  fitter->Set3DBinsNum(mQoutNbins,mQsideNbins,mQlongNbins);
  fitter->Set3DHistoRange(mQoutLo,mQoutHi,
			  mQsideLo,mQsideHi,
			  mQlongLo,mQlongHi);
  fitter->UseCoulombCorrection(true);  //on/off Coulomb correction

  fitter->Set3DNumer(hNumer);
  fitter->Set3DDenom(hDenom);
  fitter->Set3DCoulomb(hQinv);
  fitter->SetCoulombSourceRadius(mCoulombSR);

  //Set fit range
  fitter->Set3DFitRange(0., 0.25, 0., 0.25, 0., 0.25);

  //Fix cross-term parameters
  fitter->FixFitParameter3D(4); //out-side
  fitter->FixFitParameter3D(5); //out-long
  fitter->FixFitParameter3D(6); //side-long

  //Main magic
  fitter->Fit3D();
  fitter->MakeProjections();
  std::cout << "Chi2: " << fitter->GetChi2() << " NDF: " << fitter->GetNDF()
	    << " Chi2/NDF: " << fitter->GetChi2PerNDF() << std::endl;

  TFile *outFile = new TFile(oFileName,"recreate");
  
  //Obtain projections
  TH1F* hCorrFctnQoutProj = fitter->Get3DOutProjection();
  TH1F* hCorrFctnQsideProj = fitter->Get3DSideProjection();
  TH1F* hCorrFctnQlongProj = fitter->Get3DLongProjection();

  TH1F* hFitQoutProj = fitter->Get3DFitOutProjection();
  TH1F* hFitQsideProj = fitter->Get3DFitSideProjection();
  TH1F* hFitQlongProj = fitter->Get3DFitLongProjection();

  TCanvas *canv = new TCanvas("canv","canv",800,600);
  canv->Divide(3,1);
  canv->cd(1);
  hCorrFctnQoutProj->Draw();
  hCorrFctnQoutProj->GetXaxis()->SetRange(0.,0.3);
  hFitQoutProj->SetLineColor(kRed);
  hFitQoutProj->Draw("sameL");
  canv->cd(2);
  hCorrFctnQsideProj->Draw();
  hCorrFctnQsideProj->GetXaxis()->SetRange(0.,0.3);
  hFitQsideProj->SetLineColor(kRed);
  hFitQsideProj->Draw("sameL");
  canv->cd(3);
  hCorrFctnQlongProj->Draw();
  hCorrFctnQlongProj->GetXaxis()->SetRange(0.,0.3);
  hFitQlongProj->SetLineColor(kRed);
  hFitQlongProj->Draw("sameL");
  
  //outFile->Write();
  hCorrFctnQoutProj->Write();
  hCorrFctnQsideProj->Write();
  hCorrFctnQlongProj->Write();
  hFitQoutProj->Write();
  hFitQsideProj->Write();
  hFitQlongProj->Write();
  outFile->Close();
  std::vector<double> mRetVal = fitter->GetFemto3DParameters();
  std::vector<double> mRetValErr = fitter->GetFemto3DParError();
  std::vector<double> mRetNonFemtoVal = fitter->GetNonFemto3DParameters();
  std::vector<double> mRetNonFemtoValErr = fitter->GetNonFemto3DParError();
  double chi2PerNDF = fitter->GetChi2PerNDF();
  std::cout << "Fit parameters:" << std::endl;
  std::cout << "Chi2/NDF: " << chi2PerNDF << std::endl;
  for(unsigned int iVal = 0; iVal < mRetVal.size(); iVal++) {
    std::cout << Form("Par: %d %4.3f+/-%4.3f",iVal,mRetVal.at(iVal),
		      mRetValErr.at(iVal)) << std::endl;
  }
  std::cout << Form("Norm: %4.3f+/-%4.3f",mRetNonFemtoVal.at(0),
		      mRetNonFemtoValErr.at(0)) << std::endl;
}