#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 = "Crab3D.root",
		  const char* oFileName = "oTest3DFitter.root") {

  Int_t mPartType = 3;
  Int_t mSourceForm = 0;
  Int_t mNonFemtoForm = 0;
  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_0_Num_bin_2");
  //TH3F* hDenom = (TH3F*)inFile->Get("hPionBPLCMSFrameKtCF_0_Den_bin_2");
  //TH3F* hQinv = (TH3F*)inFile->Get("hPionBPLCMSFrameKtCF_0_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(false);  //on/off Coulomb correction

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

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

  //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();
  hFitQoutProj->SetLineColor(kRed);
  hFitQoutProj->Draw("sameL");
  canv->cd(2);
  hCorrFctnQsideProj->Draw();
  hFitQsideProj->SetLineColor(kRed);
  hFitQsideProj->Draw("sameL");
  canv->cd(3);
  hCorrFctnQlongProj->Draw();
  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();
  double chi2PerNDF = fitter->GetChi2PerNDF();
  std::cout << "Fit parameters:" << std::endl;
  std::cout << "Chi2/NDF: " << chi2PerNDF << std::endl;
  char*  mParNames[14] = {"Lambda",       "R_{out}",      "R_{side}",       "R_{long}", 
			  "R_{out-side}", "R_{out-long}", "R_{side-long}",
			  "Norm",         "A_{out}",      "A_{side}",       "A_{long}", 
			  "B_{out}",      "B_{side}",     "B_{long}"};
  for(unsigned int iVal = 0; iVal < mRetVal.size(); iVal++) {
    std::cout << Form("%s: %4.3f+/-%4.3f",mParNames[iVal],mRetVal.at(iVal),
		      mRetValErr.at(iVal)) << std::endl;
  }
}