#include "HistoCollector3D.h"
#include "HbtFitter.h"
#include <TCanvas.h>
#include <TLegend.h>
#include <TMath.h>
#include <TROOT.h>
#include <TGraphErrors.h>
#include <vector>
#include <iostream>
#include <fstream>

using namespace std;

int PartType = 4;
int qInvBins = 40;
double qInvLo = 0.;
double qInvHi = 0.4;
int NKtBins = 11;
int SourceFunc = 0;   //Gauss
int NonFemtoFunc = 0; //Constant
double FitRangeLo = 0.;
double FitRangeHi = 0.2;
double CoulombSR = 5.;

HbtFitter *DoFit(TH3F* hNum, TH3F* hDen, TH3F* hQinv, 
		 int aPartType=3, int aSourceFunc=0, 
		 int aNonFemtoFunc=0);
void SetTGraphStyle(TGraphErrors *gr, int charge);
void SetProjectionStyle(TH1* h, int proj, Bool_t isFit=false);

//_________________
void ProcessMesons3D(int PID=3, int cent=8, int energy=200) {

  double PionMassSqr = 0.019479835;
  double KaonMassSqr = 0.24371698;
  HistoCollector3D mHistoCollector;
  PartType = PID;
  int Centrality = cent;
  TString inFileName, oFileName, oTxtFileName;
  inFileName = "/home/nigmatkulov/work/star/data/bes/cent";
  inFileName += cent;
  inFileName += "/auau";
  oFileName = "oAuAu";
  oTxtFileName = "oAuAu";
  inFileName += energy;
  oFileName += energy;
  oTxtFileName += energy;
  if(PartType == 3) {
    inFileName += "_pions_";
    oFileName += "_pions_";
    oTxtFileName += "_pions_";
  }
  else if(PartType == 4) {
    inFileName += "_kaons_";
    oFileName += "_kaons_";
    oTxtFileName += "_kaons_";
  }
  inFileName += "3d_cent";
  oFileName += "3d_cent";
  oTxtFileName += "3d_cent";
  inFileName += cent;
  oFileName += cent;
  oTxtFileName += cent;
  inFileName += ".root";
  oFileName += ".root";
  oTxtFileName += ".txt";
  //Here... add centralitRoutKtGrErr_0y and make output
  
  //Pions
  const char* fname;
  if(PartType == 3) {
    fname = "auau200_pions_3d_cent8.root";
  }
  //Kaons
  if(PartType == 4) {
    fname = "auau200_kaons_3d_cent8.root";
  }
  mHistoCollector.SetFileName(inFileName.Data());
  mHistoCollector.SetQinvBins(qInvBins,qInvLo,qInvHi);
  mHistoCollector.SetPartType(PartType);
  mHistoCollector.SetKtBins(NKtBins);
  mHistoCollector.LoadData();
  int nKtSumBins;
  if(PartType == 3) {
    nKtSumBins = 7;
  }
  else if(PartType == 4) {
    nKtSumBins = 5;
  }
  //AuAu200 pions
  //int FitKtArray[7][2] = {{2,2}, {3,3}, {4,4}, {5,5}, {6,6}, {7,7}, {8,10} };
  //AuAu200 kaons
  int FitKtArray[5][2] = {{1,3}, {4,5}, {5,6}, {7,8}, {8,10} };
  std::vector<TH3F *> hNumer;
  std::vector<TH3F *> hDenom;
  std::vector<TH3F *> hQinv;
  for(int iCharge=0; iCharge<3; iCharge++) {
    for(int iKtBin=0; iKtBin<nKtSumBins; iKtBin++) {
      TH3F* hA = mHistoCollector.BuildKt3DNum(iCharge,FitKtArray[iKtBin][0],FitKtArray[iKtBin][1]);
      TH3F* hB = mHistoCollector.BuildKt3DDen(iCharge,FitKtArray[iKtBin][0],FitKtArray[iKtBin][1]);
      TH3F* hC = mHistoCollector.BuildKt3DQinv(iCharge,FitKtArray[iKtBin][0],FitKtArray[iKtBin][1]);
      hNumer.push_back(hA);
      hDenom.push_back(hB);
      hQinv.push_back(hC);
    }
  }
  
  vector<double> LamVect, RoutVect, RsideVect, RlongVect, Chi2Vect;
  vector<double> LamErrVect, RoutErrVect, RsideErrVect, RlongErrVect;
  vector<double> retParVect;
  
  //HbtFitter *myFitter = NULL;
  double LamArr[nKtSumBins];
  double RoutArr[nKtSumBins];
  double RsideArr[nKtSumBins];
  double RlongArr[nKtSumBins];
  double LamErrArr[nKtSumBins];
  double RoutErrArr[nKtSumBins];
  double RsideErrArr[nKtSumBins];
  double RlongErrArr[nKtSumBins];
  double RosRatArr[nKtSumBins];
  double RosRatErrArr[nKtSumBins];
  double RosDiffArr[nKtSumBins];
  double RosDiffErrArr[nKtSumBins];
  double KtArr[nKtSumBins];
  double KtErrArr[nKtSumBins];
  double MtArr[nKtSumBins];
  double KtShift = 0.1;
  double KtStep = 0.1;
  double KtLo,KtHi,KtVal,MtVal;
  double MassSqr;
  if(PartType == 3) {
    MassSqr = PionMassSqr;
  }
  else if(PartType == 4) {
    MassSqr = KaonMassSqr;
  }

  //Kt dependence
  TGraphErrors *LamKtGrErr[3];
  TGraphErrors *RoutKtGrErr[3];
  TGraphErrors *RsideKtGrErr[3];
  TGraphErrors *RlongKtGrErr[3];
  TGraphErrors *RosRatKtGrErr[3];
  TGraphErrors *RosDiffKtGrErr[3];

  //Mt dependence
  TGraphErrors *LamMtGrErr[3];
  TGraphErrors *RoutMtGrErr[3];
  TGraphErrors *RsideMtGrErr[3];
  TGraphErrors *RlongMtGrErr[3];
  TGraphErrors *RosRatMtGrErr[3];
  TGraphErrors *RosDiffMtGrErr[3];

  vector<TH1F *> hOutProjVect,hSideProjVect,hLongProjVect;
  vector<TH1F *> hFitOutProjVect,hFitSideProjVect,hFitLongProjVect;

  //For all charges: +, -, sum
  for(int iCharge=0; iCharge<3; iCharge++) {
    for(int iKtBin=0; iKtBin<nKtSumBins; iKtBin++) {
      
      int iter = iCharge*nKtSumBins + iKtBin;
      HbtFitter *myFitter = DoFit(hNumer.at(iter),hDenom.at(iter),hQinv.at(iter),
				  PartType,SourceFunc,NonFemtoFunc);

      TH1F *hExpOutProj = myFitter->Get3DOutProjection();
      TH1F *hExpSideProj = myFitter->Get3DSideProjection();
      TH1F *hExpLongProj = myFitter->Get3DLongProjection();

      TH1F *hFitOutProj = myFitter->Get3DFitOutProjection();
      TH1F *hFitSideProj = myFitter->Get3DFitSideProjection();
      TH1F *hFitLongProj = myFitter->Get3DFitLongProjection();

      hExpOutProj->SetName(Form("hOutProj_%d_%d",iCharge,iKtBin));
      hExpSideProj->SetName(Form("hSideProj_%d_%d",iCharge,iKtBin));
      hExpLongProj->SetName(Form("hLongProj_%d_%d",iCharge,iKtBin));

      hFitOutProj->SetName(Form("hFitOutProj_%d_%d",iCharge,iKtBin));
      hFitSideProj->SetName(Form("hFitSideProj_%d_%d",iCharge,iKtBin));
      hFitLongProj->SetName(Form("hFitLongProj_%d_%d",iCharge,iKtBin));

      hOutProjVect.push_back(hExpOutProj);
      hSideProjVect.push_back(hExpSideProj);
      hLongProjVect.push_back(hExpLongProj);

      hFitOutProjVect.push_back(hFitOutProj);
      hFitSideProjVect.push_back(hFitSideProj);
      hFitLongProjVect.push_back(hFitLongProj);
      
      //Calculate kT and mT binning
      KtLo=KtShift + FitKtArray[iKtBin][0]*KtStep;
      KtHi=KtShift + FitKtArray[iKtBin][1]*KtStep + KtStep;
      KtVal = (KtHi + KtLo)*0.5;
      MtVal = TMath::Sqrt(KtVal*KtVal + MassSqr);

      KtArr[iKtBin] = KtVal;
      KtErrArr[iKtBin] = 0.01;
      MtArr[iKtBin] = MtVal;

      //Obtain fit parameters
      retParVect.clear();
      retParVect = myFitter->GetFemto3DParameters();

      Chi2Vect.push_back(myFitter->GetChi2PerNDF());
      LamVect.push_back(retParVect.at(0));
      RoutVect.push_back(retParVect.at(1));
      RsideVect.push_back(retParVect.at(2));
      RlongVect.push_back(retParVect.at(3));

      /*
      ///////// DEBUG ///////////
      std::cout << "-------------------------------------------------------------" << std::endl;
      std::cout << "Rout: " << RoutVect.back() << " Rside: " << RsideVect.back() 
		<< " Rlong: " << RlongVect.back() << std::endl;
      std::cout << "-------------------------------------------------------------" << std::endl;
      */

      //Obtain fit error parameters
      retParVect.clear();
      retParVect = myFitter->GetFemto3DParError();

      LamErrVect.push_back(retParVect.at(0));
      RoutErrVect.push_back(retParVect.at(1));
      RsideErrVect.push_back(retParVect.at(2));
      RlongErrVect.push_back(retParVect.at(3));

      /*
      ///////// DEBUG ///////////
      std::cout << "-------------------------------------------------------------" << std::endl;
      std::cout << "RoutErr: " << RoutErrVect.back() << " RsideErr: " << RsideErrVect.back() 
		<< " RlongErr: " << RlongErrVect.back() << std::endl;
      std::cout << "-------------------------------------------------------------" << std::endl;
      */

      //Parameters for the given charge and each kT bin
      LamArr[iKtBin] = LamVect.back();
      LamErrArr[iKtBin] = LamErrVect.back();
      RoutArr[iKtBin] = RoutVect.back();
      RoutErrArr[iKtBin] = RoutErrVect.back();
      RsideArr[iKtBin] = RsideVect.back();
      RsideErrArr[iKtBin] = RsideErrVect.back();
      RlongArr[iKtBin] = RlongVect.back();
      RlongErrArr[iKtBin] = RlongErrVect.back();
      RosRatArr[iKtBin] = RoutArr[iKtBin]/RsideArr[iKtBin];
      RosRatErrArr[iKtBin] = ( (RoutArr[iKtBin] * RsideErrArr[iKtBin] +
				RsideArr[iKtBin] * RoutErrArr[iKtBin]) /
			       (RsideArr[iKtBin] * RsideArr[iKtBin]) );
      RosDiffArr[iKtBin] = (RoutArr[iKtBin] * RoutArr[iKtBin] - 
			    RsideArr[iKtBin] * RsideArr[iKtBin]);
      RosDiffErrArr[iKtBin] = (2 * RoutArr[iKtBin] * RoutErrArr[iKtBin] -
			       2 * RsideArr[iKtBin] * RsideErrArr[iKtBin]);
      
    } //for(int iKtBin=0; iKtBin<nKtSumBins; iKtBin++)


    /*
    ///////// DEBUG ///////////
    std::cout << "Rout array" << std::endl;
    for(int iKt=0; iKt<nKtSumBins; iKt++) {
      std::cout << RoutArr[iKt] << "\t";
    }
    std::cout << std::endl;
    std::cout << "Rout error array" << std::endl;
    for(int iKt=0; iKt<nKtSumBins; iKt++) {
      std::cout << RoutErrArr[iKt] << "\t";
    }
    std::cout << std::endl;
    */

    //Lambda
    LamKtGrErr[iCharge] = new TGraphErrors(nKtSumBins,KtArr,LamArr,KtErrArr,LamErrArr);
    LamKtGrErr[iCharge]->SetName(Form("LamKtGrErr_%d",iCharge));
    SetTGraphStyle(LamKtGrErr[iCharge], iCharge);
    LamKtGrErr[iCharge]->GetXaxis()->SetTitle("k_{T} (GeV/c)");
    LamKtGrErr[iCharge]->GetYaxis()->SetTitle("#lambda");
    LamMtGrErr[iCharge] = new TGraphErrors(nKtSumBins,MtArr,LamArr,KtErrArr,LamErrArr);
    LamMtGrErr[iCharge]->SetName(Form("LamMtGrErr_%d",iCharge));
    SetTGraphStyle(LamMtGrErr[iCharge], iCharge);
    LamMtGrErr[iCharge]->GetXaxis()->SetTitle("m_{T} (GeV/c^{2})");
    LamMtGrErr[iCharge]->GetYaxis()->SetTitle("#lambda");

    //Rout
    RoutKtGrErr[iCharge] = new TGraphErrors(nKtSumBins,KtArr,RoutArr,KtErrArr,RoutErrArr);
    RoutKtGrErr[iCharge]->SetName(Form("RoutKtGrErr_%d",iCharge));
    SetTGraphStyle(RoutKtGrErr[iCharge], iCharge);
    RoutKtGrErr[iCharge]->GetXaxis()->SetTitle("k_{T} (GeV/c)");
    RoutKtGrErr[iCharge]->GetYaxis()->SetTitle("R_{out} (fm)");
    RoutMtGrErr[iCharge] = new TGraphErrors(nKtSumBins,MtArr,RoutArr,KtErrArr,RoutErrArr);
    RoutMtGrErr[iCharge]->SetName(Form("RoutMtGrErr_%d",iCharge));
    SetTGraphStyle(RoutMtGrErr[iCharge], iCharge);
    RoutMtGrErr[iCharge]->GetXaxis()->SetTitle("m_{T} (GeV/c^{2})");
    RoutMtGrErr[iCharge]->GetYaxis()->SetTitle("R_{out} (fm)");

    //Rside
    RsideKtGrErr[iCharge] = new TGraphErrors(nKtSumBins,KtArr,RsideArr,KtErrArr,RsideErrArr);
    RsideKtGrErr[iCharge]->SetName(Form("RsideKtGrErr_%d",iCharge));
    SetTGraphStyle(RsideKtGrErr[iCharge], iCharge);
    RsideKtGrErr[iCharge]->GetXaxis()->SetTitle("k_{T} (GeV/c)");
    RsideKtGrErr[iCharge]->GetYaxis()->SetTitle("R_{side} (fm)");
    RsideMtGrErr[iCharge] = new TGraphErrors(nKtSumBins,MtArr,RsideArr,KtErrArr,RsideErrArr);
    RsideMtGrErr[iCharge]->SetName(Form("RsideMtGrErr_%d",iCharge));
    SetTGraphStyle(RsideMtGrErr[iCharge], iCharge);
    RsideMtGrErr[iCharge]->GetXaxis()->SetTitle("m_{T} (GeV/c^{2})");
    RsideMtGrErr[iCharge]->GetYaxis()->SetTitle("R_{side} (fm)");

    //Rlong
    RlongKtGrErr[iCharge] = new TGraphErrors(nKtSumBins,KtArr,RlongArr,KtErrArr,RlongErrArr);
    RlongKtGrErr[iCharge]->SetName(Form("RlongKtGrErr_%d",iCharge));
    SetTGraphStyle(RlongKtGrErr[iCharge], iCharge);
    RlongKtGrErr[iCharge]->GetXaxis()->SetTitle("k_{T} (GeV/c)");
    RlongKtGrErr[iCharge]->GetYaxis()->SetTitle("R_{long} (fm)");
    RlongMtGrErr[iCharge] = new TGraphErrors(nKtSumBins,MtArr,RlongArr,KtErrArr,RlongErrArr);
    RlongMtGrErr[iCharge]->SetName(Form("RlongMtGrErr_%d",iCharge));
    SetTGraphStyle(RlongMtGrErr[iCharge], iCharge);
    RlongMtGrErr[iCharge]->GetXaxis()->SetTitle("m_{T} (GeV/c^{2})");
    RlongMtGrErr[iCharge]->GetYaxis()->SetTitle("R_{long} (fm)");

    //RosRat
    RosRatKtGrErr[iCharge] = new TGraphErrors(nKtSumBins,KtArr,RosRatArr,KtErrArr,RosRatErrArr);
    RosRatKtGrErr[iCharge]->SetName(Form("RosRatKtGrErr_%d",iCharge));
    SetTGraphStyle(RosRatKtGrErr[iCharge], iCharge);
    RosRatKtGrErr[iCharge]->GetXaxis()->SetTitle("k_{T} (GeV/c)");
    RosRatKtGrErr[iCharge]->GetYaxis()->SetTitle("R_{o}/R_{s}");
    RosRatMtGrErr[iCharge] = new TGraphErrors(nKtSumBins,MtArr,RosRatArr,KtErrArr,RosRatErrArr);
    RosRatMtGrErr[iCharge]->SetName(Form("RosRatMtGrErr_%d",iCharge));
    SetTGraphStyle(RosRatMtGrErr[iCharge], iCharge);
    RosRatMtGrErr[iCharge]->GetXaxis()->SetTitle("m_{T} (GeV/c^{2})");
    RosRatMtGrErr[iCharge]->GetYaxis()->SetTitle("R_{o}/R_{s}");

    //RosDiff
    RosDiffKtGrErr[iCharge] = new TGraphErrors(nKtSumBins,KtArr,RosDiffArr,KtErrArr,RosDiffErrArr);
    RosDiffKtGrErr[iCharge]->SetName(Form("RosDiffKtGrErr_%d",iCharge));
    SetTGraphStyle(RosDiffKtGrErr[iCharge], iCharge);
    RosDiffKtGrErr[iCharge]->GetXaxis()->SetTitle("k_{T} (GeV/c)");
    RosDiffKtGrErr[iCharge]->GetYaxis()->SetTitle("R_{o}^{2}-R_{s}^{2} (fm)");
    RosDiffMtGrErr[iCharge] = new TGraphErrors(nKtSumBins,MtArr,RosDiffArr,KtErrArr,RosDiffErrArr);
    RosDiffMtGrErr[iCharge]->SetName(Form("RosDiffMtGrErr_%d",iCharge));
    SetTGraphStyle(RosDiffMtGrErr[iCharge], iCharge);
    RosDiffMtGrErr[iCharge]->GetXaxis()->SetTitle("m_{T} (GeV/c^{2})");
    RosDiffMtGrErr[iCharge]->GetYaxis()->SetTitle("R_{o}^{2}-R_{s}^{2} (fm)");

  } //for(unsigned int iCharge=0; iCharge<3; iCharge++)

  TFile *oFile = new TFile(oFileName.Data(),"recreate");

  for(unsigned int iHisto = 0; iHisto<hOutProjVect.size(); iHisto++) {
    SetProjectionStyle(hOutProjVect.at(iHisto),0);
    hOutProjVect.at(iHisto)->Write();
    SetProjectionStyle(hSideProjVect.at(iHisto),1);
    hSideProjVect.at(iHisto)->Write();
    SetProjectionStyle(hLongProjVect.at(iHisto),2);
    hLongProjVect.at(iHisto)->Write();

    SetProjectionStyle(hFitOutProjVect.at(iHisto),0,true);
    hFitOutProjVect.at(iHisto)->Write();
    SetProjectionStyle(hFitSideProjVect.at(iHisto),1,true);
    hFitSideProjVect.at(iHisto)->Write();
    SetProjectionStyle(hFitLongProjVect.at(iHisto),2,true);
    hFitLongProjVect.at(iHisto)->Write();
  }
  
  for(int iCharge=0; iCharge<3; iCharge++) {
    LamKtGrErr[iCharge]->Write();
    LamMtGrErr[iCharge]->Write();
    RoutKtGrErr[iCharge]->Write();
    RoutMtGrErr[iCharge]->Write();
    RsideKtGrErr[iCharge]->Write();
    RsideMtGrErr[iCharge]->Write();
    RlongKtGrErr[iCharge]->Write();
    RlongMtGrErr[iCharge]->Write();
    RosRatKtGrErr[iCharge]->Write();
    RosRatMtGrErr[iCharge]->Write();
    RosDiffKtGrErr[iCharge]->Write();
    RosDiffMtGrErr[iCharge]->Write();
  }

  ofstream oTxtFile;
  oTxtFile.open(oTxtFileName.Data());
  oTxtFile << Form("Collision energy: %3i \t Centrality: %1i \t Coulomb R: %2.1f \n",
		   energy, cent, CoulombSR);
  oTxtFile << "-------------------------------\n";
  int point;
  Double_t x,y,ey;

  //Write to file
  for(int iCharge=0; iCharge<3; iCharge++) {

    switch(iCharge) {
    case 0:
      oTxtFile << "Positive charge: \n";
      break;
    case 1:
      oTxtFile << "Negative charge: \n";
      break;
    case 2:
      oTxtFile << "Sum of charges: \n";
      break;
    default:
      oTxtFile << "Error in the charge cycle \n";
      break;
    };
    
    oTxtFile << "-------------------------------\n";
    for(int iKtBin=0; iKtBin<nKtSumBins; iKtBin++) {
      if(iKtBin==0) {
	oTxtFile << "lambda: ";
      }
      point = LamKtGrErr[iCharge]->GetPoint(iKtBin,x,y);
      ey =  LamKtGrErr[iCharge]->GetErrorY(iKtBin);
      oTxtFile << Form(" %5.4f +- %5.4f ", y, ey);
      if(iKtBin == (nKtSumBins-1)) {
	oTxtFile << "\n";
      }
      else {
	oTxtFile << " | ";
      }
    }
    for(int iKtBin=0; iKtBin<nKtSumBins; iKtBin++) {
      if(iKtBin==0) {
	oTxtFile << "Rout  : ";
      }
      point = RoutKtGrErr[iCharge]->GetPoint(iKtBin,x,y);
      ey =  RoutKtGrErr[iCharge]->GetErrorY(iKtBin);
      oTxtFile << Form(" %5.4f +- %5.4f ", y, ey);
      if(iKtBin == (nKtSumBins-1)) {
	oTxtFile << "\n";
      }
      else {
	oTxtFile << " | ";
      }
    }
    for(int iKtBin=0; iKtBin<nKtSumBins; iKtBin++) {
      if(iKtBin==0) {
	oTxtFile << "Rside : ";
      }
      point = RsideKtGrErr[iCharge]->GetPoint(iKtBin,x,y);
      ey =  RsideKtGrErr[iCharge]->GetErrorY(iKtBin);
      oTxtFile << Form(" %5.4f +- %5.4f ", y, ey);
      if(iKtBin == (nKtSumBins-1)) {
	oTxtFile << "\n";
      }
      else {
	oTxtFile << " | ";
      }
    }
    for(int iKtBin=0; iKtBin<nKtSumBins; iKtBin++) {
      if(iKtBin==0) {
	oTxtFile << "Rlong : ";
      }
      point = RlongKtGrErr[iCharge]->GetPoint(iKtBin,x,y);
      ey =  RlongKtGrErr[iCharge]->GetErrorY(iKtBin);
      oTxtFile << Form(" %5.4f +- %5.4f ", y, ey);
      if(iKtBin == (nKtSumBins-1)) {
	oTxtFile << "\n";
      }
      else {
	oTxtFile << " | ";
      }
    }
    for(int iKtBin=0; iKtBin<nKtSumBins; iKtBin++) {
      if(iKtBin==0) {
	oTxtFile << "RosRat : ";
      }
      point = RosRatKtGrErr[iCharge]->GetPoint(iKtBin,x,y);
      ey =  RosRatKtGrErr[iCharge]->GetErrorY(iKtBin);
      oTxtFile << Form(" %5.4f +- %5.4f ", y, ey);
      if(iKtBin == (nKtSumBins-1)) {
	oTxtFile << "\n";
      }
      else {
	oTxtFile << " | ";
      }
    }
    for(int iKtBin=0; iKtBin<nKtSumBins; iKtBin++) {
      if(iKtBin==0) {
	oTxtFile << "RosDiff : ";
      }
      point = RosDiffKtGrErr[iCharge]->GetPoint(iKtBin,x,y);
      ey =  RosDiffKtGrErr[iCharge]->GetErrorY(iKtBin);
      oTxtFile << Form(" %5.4f +- %5.4f ", y, ey);
      if(iKtBin == (nKtSumBins-1)) {
	oTxtFile << "\n";
      }
      else {
	oTxtFile << " | ";
      }
    }
    oTxtFile << "-------------------------------\n";
  } //for(int iCharge=0; iCharge<3; iCharge++)

  oTxtFile << "\n";
  oTxtFile.close();

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

//_________________
HbtFitter *DoFit(TH3F* hNum, TH3F* hDen, TH3F* hQinv, 
		 int aPartType, int aSourceFunc, int aNonFemtoFunc) {
  
  HbtFitter *fitter = new HbtFitter(aPartType,aSourceFunc,aNonFemtoFunc,false);
  fitter->SetFitMethod(1);  //Log-likelihood
  fitter->SetCorrectOnMc(false);
  fitter->SetProjectionType(false); //Project by vals
  fitter->Set3DBinsNum(qInvBins,qInvBins,qInvBins);
  fitter->Set3DHistoRange(qInvLo,qInvHi,
			  qInvLo,qInvHi,
			  qInvLo,qInvHi);
  fitter->UseCoulombCorrection(true);  //on/off Coulomb correction
  fitter->SetCoulombSourceRadius(CoulombSR);

  hQinv->Divide(hDen);      //!!!!!! Hey dude, do not forget to divide !!!!!!!
  fitter->Set3DNumer(hNum);
  fitter->Set3DDenom(hDen);
  fitter->Set3DCoulomb(hQinv);

  //Set fit range
  fitter->Set3DFitRange(FitRangeLo, FitRangeHi, 
			FitRangeLo, FitRangeHi, 
			FitRangeLo, FitRangeHi);
  
  fitter->SetInit3DFitParameters(0.3, 3., 3., 3., 0., 0., 0., 0.1,
				 0., 0., 0., 0., 0., 0.);
  fitter->SetInit3DFitParErrors(0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.001,
				0.1, 0.1, 0.1, 0.1, 0.1, 0.1);
  fitter->SetInit3DFitParLimitLo(0.05, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.,
				 -2., -2., -2., -2., -2., -2.);
  fitter->SetInit3DFitParLimitHi(1.05, 10., 10., 10., 10., 10., 10., 10.,
				 +2., +2., +2., +2., +2., +2.);
  //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;
  
  return fitter;
}

//_________________
void SetTGraphStyle(TGraphErrors *gr, int charge) {
  
  Int_t mColor;
  Int_t mMarkerStyle;
  Double_t mMarkerSize = 1.5;
  Int_t mLineWidth = 2;
  
  if(charge == 0) {
    mColor = 2;       //Red
    mLineWidth = 2;
    if(PartType == 3) {
      mMarkerStyle = 20;
    }
    else if(PartType == 4) {
      mMarkerStyle = 29;
    }
  }
  else if(charge == 1) {
    mColor = 4;       //Blue
    mLineWidth = 2;
    if(PartType == 3) {
      mMarkerStyle = 20;
    }
    else if(PartType == 4) {
      mMarkerStyle = 29;
    }
  }
  else {
    mColor = 1;       //Black
    mLineWidth = 2;
    if(PartType == 3) {
      mMarkerStyle = 20;
    }
    else if(PartType == 4) {
      mMarkerStyle = 29;
    }
  }

  gr->SetLineColor(mColor);
  gr->SetMarkerColor(mColor);
  gr->SetLineWidth(mLineWidth);
  gr->SetMarkerStyle(mMarkerStyle);
  gr->SetMarkerSize(mMarkerSize);
}

//_________________
void SetProjectionStyle(TH1* h, int proj, Bool_t isFit) {
  TString Xstr,Ystr;
  if(proj == 0) {
    Xstr = "q_{out} (GeV/c)";
    Ystr =  "C(q_{out})";
  }
  else if(proj == 1) {
    Xstr = "q_{side} (GeV/c)";
    Ystr =  "C(q_{side})";
  }
  else if(proj == 2) {
    Xstr = "q_{long} (GeV/c)";
    Ystr =  "C(q_{long})";
  }
  int color;
  if(isFit == false) {
    color = 1;
  }
  else {
    color = 2;
  }
  
  h->SetLineWidth(2);
  h->SetLineColor(color);
  h->SetMarkerColor(color);
  h->SetAxisRange(0.9, 1.6,"Y");
  h->SetAxisRange(FitRangeLo,FitRangeHi,"X");
  h->GetXaxis()->SetTitle(Xstr.Data());
  h->GetYaxis()->SetTitle(Ystr.Data());

  h->SetNdivisions(505,"X");
  h->SetNdivisions(505,"Y");
}