my_star/source/HbtFitter/HbtFitter.cxx

#include "HbtFitter.h"
#include <TMath.h>
#include <TString.h>
#include <iostream>

Double_t R_GAUS_TO_EXPO = 1./TMath::Sqrt( TMath::Pi() );
Double_t hBar = 0.1973269718;

//_________________
HbtFitter* HbtFitter::objHbtFitter=0;

//_________________
HbtFitter::HbtFitter(Int_t PartType, Int_t SourceForm,
                     Int_t NonFemtoForm, Bool_t ExcludePoints) {

  //Set basic parameters
  SetPartType(PartType);
  SetSourceForm(SourceForm);
  SetNonFemtoForm(NonFemtoForm);
  SetExcludeFromFit(ExcludePoints);
  SetFitMethod(0);

  //Correct on MC or not
  mCorrectOnMc = false;

  //Project by bins or by vals
  mProjectionType = false;

  //Number of femtoscopic parameters: lambda, R
  mFemto1DParNum = 2;
  mFemto3DParNum = 7;

  //Fixed parameter vector initialization
  if(!mFixParameters1D.empty()) {
    mFixParameters1D.clear();
  }
  if(!mFixParameters3D.empty()) {
    mFixParameters3D.clear();
  }

  //Initialize fit parameters
  mQinvVal[0] = 0.; mQinvVal[1] = 1.;
  mQoutVal[0] = 0.; mQoutVal[1] = 1.;
  mQsideVal[0] = 0.; mQsideVal[1] = 1.;
  mQlongVal[0] = 0.; mQlongVal[1] = 1.;

  //Number of Q bins
  mQinvNbins = 50;
  mQoutNbins = 50;
  mQsideNbins = 50;
  mQlongNbins = 50;

  //Projection ranges (in bins and in vals)
  mOutProjBins[0] = 1; mOutProjBins[1] = 6;
  mSideProjBins[0] = 1; mSideProjBins[1] = 6;
  mLongProjBins[0] = 1; mLongProjBins[1] = 6;

  mOutProjVals[0] = 0.; mOutProjVals[1] = 0.05;
  mSideProjVals[0] = 0.; mSideProjVals[1] = 0.05;
  mLongProjVals[0] = 0.; mLongProjVals[1] = 0.05;

  //Clean fit parameters
  if(!mFemto1DVal.empty()) {
    mFemto1DVal.clear();
  }
  if(!mFemto1DValErr.empty()) {
    mFemto1DValErr.clear();
  }
  if(!mNonFemto1DVal.empty()) {
    mNonFemto1DVal.clear();
  }
  if(!mNonFemto1DValErr.empty()) {
    mNonFemto1DValErr.clear();
  }

  if(!mFemto3DVal.empty()) {
    mFemto3DVal.clear();
  }
  if(!mFemto3DValErr.empty()) {
    mFemto3DValErr.clear();
  }
  if(!mNonFemto3DVal.empty()) {
    mNonFemto3DVal.clear();
  }
  if(!mNonFemto3DValErr.empty()) {
    mNonFemto3DValErr.clear();
  }

  //Use Coulomb correction if it is not set
  mUseCoulombCorrection = true;

  //Initialize 1D non-femto parameters
  for(Int_t iIter=0; iIter<3; iIter++) {
    mNonFemto1DParArr[iIter] = 0;
    mNonFemto1DParErrArr[iIter] = 0;
  }

  //Initialize 3D non-femto parameters
  for(Int_t iIter1=0; iIter1<3; iIter1++) {  //Out,Side,Long
    //iIter1 corresponds to out,side,long
    for(Int_t iIter2=0; iIter2<3; iIter2++) {
      //iIter2 corresponds to the mNonFemtoParNum
      mNonFemto3DParArr[iIter1][iIter2] = 0;
      mNonFemto3DParErrArr[iIter1][iIter2] = 0;
    }
  }

  //Set null pointers to the histograms
  h1dNumer = NULL;
  h1dDenom = NULL;
  h1dCorrFctn = NULL;
  h1dMcNumer = NULL;
  h1dMcDenom = NULL;
  h1dMcCorrFctn = NULL;

  h3dNumer = NULL;
  h3dDenom = NULL;
  h3dCorrFctn = NULL;
  h3dCoulombWeight = NULL;
  h3dMcNumer = NULL;
  h3dMcDenom = NULL;
  h3dMcCorrFctn = NULL;

  //3D fit numerator
  h3dFitNumerator = NULL;
  h3dFitQoutProj = NULL;
  h3dFitQsideProj = NULL;
  h3dFitQlongProj = NULL;

  //3D correlation function and 3D fit projections
  h3dQoutProjNum = NULL;
  h3dQoutProjDen = NULL;
  h3dMcQoutProjNum = NULL;
  h3dMcQoutProjDen = NULL;

  h3dQsideProjNum = NULL;
  h3dQsideProjDen = NULL;
  h3dMcQsideProjNum = NULL;
  h3dMcQsideProjDen = NULL;
  
  h3dQlongProjNum = NULL;
  h3dQlongProjDen = NULL;
  h3dMcQlongProjNum = NULL;
  h3dMcQlongProjDen = NULL;
  
  h3dCorrFctnOutProj = NULL;
  h3dCorrFctnFitOutProj = NULL;
  h3dCorrFctnSideProj = NULL;
  h3dCorrFctnFitSideProj = NULL;
  h3dCorrFctnLongProj = NULL;
  h3dCorrFctnFitLongProj = NULL;
  
  //Pointer to the current object
  SetObject(this);
}

//_________________
HbtFitter::~HbtFitter() {
  /* empty */
}

//_________________
void HbtFitter::SetCorrectOnMc(Bool_t corr) {
  mCorrectOnMc = corr;
  std::cout << "HbtFitter::SetCorrectOnMc : " << mCorrectOnMc << std::endl;
}

//_________________
void HbtFitter::SetSourceForm(Int_t SourceForm) {

  mSourceForm = SourceForm;
  std::cout << "HbtFitter::SetSourceForm : " << mSourceForm << std::endl;
}

//_________________
void HbtFitter::SetNonFemtoForm(Int_t form) {

  /*
    Non-femtoscopic forms:
    0 - constant
    1 - linear
    2 - 2nd oreder polynomial
    3 - sqrt(1 + a x^2 + b x^4)  (sqrt tail)
    4 - a ( 1 + e^{b x^2})       (Gaussian tail)
    5 - a / (1 + b x^2)          (hyperbolic tail)
   */

  mNonFemtoForm = form;

  switch(mNonFemtoForm) {
  case 0 :
    mNonFemto1DParNum = 1;
    mNonFemto3DParNum = 1;
    break;
  case 1 :
    mNonFemto1DParNum = 2;
    mNonFemto3DParNum = 4;
    break;
  case 2 :
    mNonFemto1DParNum = 3;
    mNonFemto3DParNum = 7;
    break;
  case 3 :
    mNonFemto1DParNum = 3;
    mNonFemto3DParNum = 7;
    break;
  case 4 :
    mNonFemto1DParNum = 2;
    mNonFemto3DParNum =4;
    break;
  case 5 :
    mNonFemto1DParNum = 2;
    mNonFemto3DParNum = 4;
    break;
  default:
    std::cout << "SetNonFemtoForm: Unknown non-femtoscopic form" << std::cout;
    mNonFemto1DParNum = 1;
    mNonFemto3DParNum = 1;
    break;
  }
  std::cout << "HbtFitter::SetNonFemtoForm : " << mNonFemtoForm << std::endl;
}

//_________________
void HbtFitter::SetFitMethod(Int_t method) {

  switch(method) {
  case 0: //ChiSquare
    mFitMethod = 0;
    std::cout << "HbtFitter::FitMethod : ChiSquare" << std::endl;
    break;
  case 1: //Log-likelihood
    mFitMethod = 1;
    std::cout << "HbtFitter::FitMethod : Log-Likelihood" << std::endl;
    break;
  default: //ChiSquare
    mFitMethod = 0;
    std::cout << "HbtFitter::FitMethod : [WARNING] - unknown fit method" 
	      << "\nReset to ChiSquare" << std::endl;
    break;
  }
}

//_________________
void HbtFitter::Set1DFitRange(Double_t qLo, Double_t qHi) {
  mQinvFitRange[0] = qLo;
  mQinvFitRange[1] = qHi;
  if(h1dNumer) {
    mQinvFitRangeBin[0] = h1dNumer->GetXaxis()->FindBin(mQinvFitRange[0]);
    mQinvFitRangeBin[1] = h1dNumer->GetXaxis()->FindBin(mQinvFitRange[1]);
  }
  else {
    std::cout << "HbtFitter::Set1DFitRange [WARNING] Fit range is set to default values" << std::endl;
    mQinvFitRangeBin[0] = 1;
    mQinvFitRangeBin[1] = mQinvNbins;
  }
}

//_________________
void HbtFitter::Set3DFitRange(Double_t qOutLo, Double_t qOutHi,
			      Double_t qSideLo, Double_t qSideHi,
			      Double_t qLongLo, Double_t qLongHi) {
  mQoutFitRange[0] = qOutLo; mQoutFitRange[1] = qOutHi;
  mQsideFitRange[0] = qSideLo; mQsideFitRange[1] = qSideHi;
  mQlongFitRange[0] = qLongLo; mQlongFitRange[1] = qLongHi;
  if(h3dNumer) {
    mQoutFitRangeBin[0] = h3dNumer->GetXaxis()->FindBin(mQoutFitRange[0]);
    mQoutFitRangeBin[1] = h3dNumer->GetXaxis()->FindBin(mQoutFitRange[1]);
    mQsideFitRangeBin[0] = h3dNumer->GetXaxis()->FindBin(mQsideFitRange[0]);
    mQsideFitRangeBin[1] = h3dNumer->GetXaxis()->FindBin(mQsideFitRange[1]);
    mQlongFitRangeBin[0] = h3dNumer->GetXaxis()->FindBin(mQlongFitRange[0]);
    mQlongFitRangeBin[1] = h3dNumer->GetXaxis()->FindBin(mQlongFitRange[1]);
  }
  else {
    std::cout << "HbtFitter::Set3DFitRange [WARNING] Fit range is set to default values" << std::endl;
    mQoutFitRangeBin[0] = 1; mQoutFitRangeBin[1] = mQoutNbins;
    mQsideFitRangeBin[0] = 1; mQsideFitRangeBin[1] = mQsideNbins;
    mQlongFitRangeBin[0] = 1; mQlongFitRangeBin[1] = mQlongNbins;
  }

  std::cout << "Fit ranges: " << std::endl
	    << "Qout : " << mQoutFitRange[0] << "<=Qout<=" << mQoutFitRange[1] 
	    << " " << mQoutFitRangeBin[0] << "<=QoutBin<=" << mQoutFitRangeBin[1] << std::endl
	    << "Qside : " << mQsideFitRange[0] << "<=Qside<=" << mQsideFitRange[1] 
	    << " " << mQsideFitRangeBin[0] << "<=QsideBin<=" << mQsideFitRangeBin[1] << std::endl
	    << "Qlong : " << mQlongFitRange[0] << "<=Qlong<=" << mQlongFitRange[1] 
	    << " " << mQlongFitRangeBin[0] << "<=QlongBin<=" << mQlongFitRangeBin[1] << std::endl;
}

//_________________
void HbtFitter::Fit1D() {


  //Double check if the points were really exluded
  /*
  if(mExcludeRange == 1) {
    if(x[0] > mExcludedVal[0] && x[0]<mExcludedVal[1]) {
      TF1::RejectPoint();
      return 0;
    }
  }
  */
  Double_t mAutoParVal[5] = {0.5, 5., 1., 0., 0.};
  Double_t mAutoParValErr[5] = {0.01, 0.01, 0.01, 0.01, 0.01};
  TString  mParNames[5] = {"#lambda", "R_{inv}", "N", "a", "b"};

  Int_t mFitParNum = mFemto1DParNum + mNonFemto1DParNum;
  //From Victor M. Perevozchikov
  SetObject(this);
  TMinuit mFit(mFitParNum);
  switch(mFitMethod) {
  case 0:
    std::cout << "Fit1D: ChiSquare method has been initialized" << std::endl;
    mFit.SetFCN(HbtFitter::FitFunc1DChiSquare);
    break;
  case 1:
    std::cout << "Fit1D: Log-likelihood method has been initialized" << std::endl;
    mFit.SetFCN(HbtFitter::FitFunc1DLogLike);
    break;
  default:
    std::cout << "Fit1D: Default method has been initialized" << std::endl;
    mFit.SetFCN(HbtFitter::FitFunc1DChiSquare);
    break;
  }

  Double_t mInitParVal[mFitParNum];
  Double_t mInitParValErr[mFitParNum];
  TString  mInitParName[mFitParNum];
  for(Int_t iIter=0; iIter<mFitParNum; iIter++) {
    mInitParVal[iIter] = mAutoParVal[iIter];
    mInitParValErr[iIter] = mAutoParValErr[iIter];
    mInitParName[iIter] = mParNames[iIter];
    mFit.DefineParameter(iIter,mInitParName[iIter].Data(),mInitParVal[iIter],mInitParValErr[iIter], 0, 0);
  }
  //Fix fit parameters
  if(!mFixParameters1D.empty()) {
    for(UInt_t iIter=0; iIter<mFixParameters1D.size(); iIter++) {
      mFit.FixParameter(mFixParameters1D.at(iIter));
    }
  }

  Double_t mArgList[mFitParNum];
  Int_t ierflag = 0;
  Int_t mNFitCycles = 0;
  do {
    mNFitCycles++;
    std::cout << "********** Fit loop " <<mNFitCycles << std::endl;
    //mFit.mnexcm("SIMPLEX", mArgList, 0, ierflag);
    mFit.mnexcm("MIGRAD", mArgList, 0, ierflag);
  } while( (ierflag!=0) && (mNFitCycles!=50) );

  Double_t mFitVal, mFitValErr;
  //Return femtoscopic fit parameters
  for(Int_t iIter=0; iIter<mFemto1DParNum; iIter++) {
    mFitVal=0.; mFitValErr=0.;
    mFit.GetParameter(iIter, mFitVal, mFitValErr);
    mFemto1DVal.push_back(mFitVal);
    mFemto1DValErr.push_back(mFitValErr);
    std::cout << Form("%s = %3.2f +- %3.2f", mParNames[iIter].Data(), 
		      mFitVal, mFitValErr) << std::endl;
  }
  //Return non-femtoscopic fit parameters
  for(Int_t iIter=mFemto1DParNum; iIter<mFitParNum; iIter++) {
    mFitVal=0.; mFitValErr=0.;
    mFit.GetParameter(iIter, mFitVal, mFitValErr);
    mNonFemto1DVal.push_back(mFitVal);
    mNonFemto1DValErr.push_back(mFitValErr);
    std::cout << Form("%s = %3.2f +- %3.2f", mParNames[iIter].Data(), 
		      mFitVal, mFitValErr) << std::endl;
  }
}

//_________________
void HbtFitter::Fit3D() {
 
  //Initial values:              L      Ro     Rs     Rl     Ros    Rol    Rsl 
  //                             N      ao     as     al     bo     bs     bl
  Double_t mAutoParVal[14] =    {0.5,   4.,    4.,    4.,    0.,    0.,    0., 
				 2.0,  0.,    0.,    0.,    0.,    0.,    0.}; 

  //                             L      Ro     Rs     Rl     Ros    Rol    Rsl 
  //                             N      ao     as     al     bo     bs     bl
  Double_t mAutoParValErr[14] = {0.01,  0.01,  0.01,  0.01,  0.01,  0.01,  0.01,
				 0.001, 0.01,  0.01,  0.01,  0.01,  0.01,  0.01};

  //                             L      Ro     Rs     Rl     Ros    Rol    Rsl 
  //                             N      ao     as     al     bo     bs     bl
  Double_t mAutoParValLo[14] =  {0.,    0.,    0.,    0.,    0.,    0.,    0., 
				 0.,    -5.,   -5.,   -5.,  -2.,  -2.,   -2.};

  //                             L      Ro     Rs     Rl     Ros    Rol    Rsl 
  //                             N      ao     as     al     bo     bs     bl
  Double_t mAutoParValHi[14] =  {1.05,  60.,   60.,   60.,   60.,   60.,   60.,
				 5.,    5.,    5.,    5.,    2.,    2.,    2.}; 
  TString  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}"};
  
  Int_t mFitParNum = mFemto3DParNum + mNonFemto3DParNum;
  TMinuit mFit(mFitParNum);
  switch(mFitMethod) {
  case 0:
    std::cout << "Fit3D: ChiSquare method has been initialized" << std::endl;
    mFit.SetFCN(HbtFitter::FitFunc3DChiSquare);
    break;
  case 1:
    std::cout << "Fit3D: Log-likelihood method has been initialized" << std::endl;
    mFit.SetFCN(HbtFitter::FitFunc3DLogLike);
    break;
  default:
    std::cout << "Fit3D: Default method has been initialized" << std::endl;
    mFit.SetFCN(HbtFitter::FitFunc3DLogLike);
    break;
  }
  Double_t mInitParVal[mFitParNum];
  Double_t mInitParValErr[mFitParNum];
  Double_t mInitParValLo[mFitParNum];
  Double_t mInitParValHi[mFitParNum];
  TString  mInitParName[mFitParNum];
  for(Int_t iIter=0; iIter<mFitParNum; iIter++) {
    mInitParVal[iIter] = mAutoParVal[iIter];
    mInitParValErr[iIter] = mAutoParValErr[iIter];
    mInitParValLo[iIter] = mAutoParValLo[iIter];
    mInitParValHi[iIter] = mAutoParValHi[iIter];
    mInitParName[iIter] = mParNames[iIter];
    mFit.DefineParameter(iIter, mInitParName[iIter].Data(),
			 mInitParVal[iIter], mInitParValErr[iIter],
			 mInitParValLo[iIter], mInitParValHi[iIter]);
  }
  //Fix fit parameters
  if(!mFixParameters3D.empty()) {
    for(UInt_t iIter=0; iIter<mFixParameters3D.size(); iIter++) {
      mFit.FixParameter(mFixParameters3D.at(iIter));
    }
  }
  mFit.SetErrorDef(0.5); //Log-likelihood minimization

  Double_t mArgList[mFitParNum];
  Int_t ierflag = 0;
  Int_t mNFitCycles = 0;
  std::cout << "Start 3D fit" << std::endl;
  do {
    mNFitCycles++;
    std::cout << "********** Fit loop " <<mNFitCycles << std::endl;
    mFit.mnexcm("MIGRAD", mArgList, 0, ierflag);
  } while( (ierflag!=0) && (mNFitCycles!=50) );
  std::cout << "The 3D fit has been finished" << std::endl;

  Double_t mFitVal, mFitValErr;
  //Return femtoscopic fit parameters
  for(Int_t iIter=0; iIter<mFemto3DParNum; iIter++) {
    mFitVal=0.; mFitValErr=0.;
    mFit.GetParameter(iIter, mFitVal, mFitValErr);
    mFemto3DVal.push_back(mFitVal);
    mFemto3DValErr.push_back(mFitValErr);
  }
  //Return non-femtoscopic fit parameters
  for(Int_t iIter=mFemto3DParNum; iIter<mFitParNum; iIter++) {
    mFitVal=0.; mFitValErr=0.;
    mFit.GetParameter(iIter, mFitVal, mFitValErr);
    mNonFemto3DVal.push_back(mFitVal);
    mNonFemto3DValErr.push_back(mFitValErr);
  }

  if(!h3dCorrFctn) {
    h3dCorrFctn =  new TH3F("h3dCorrFctn","h3dCorrFctn",
			    mQoutNbins,mQoutVal[0],mQoutVal[1],
			    mQsideNbins,mQsideVal[0],mQsideVal[1],
			    mQlongNbins,mQlongVal[0],mQlongVal[1]);
    h3dCorrFctn->Divide(h3dNumer,h3dDenom);
  }
}

//_________________
void HbtFitter::Set3DBinsNum(Int_t OutNBins, Int_t SideNBins, Int_t LongNBins) {
  mQoutNbins = OutNBins;
  mQsideNbins = SideNBins;
  mQlongNbins = LongNBins;
}

//_________________
void HbtFitter::Set3DHistoRange(Double_t outLo, Double_t outHi,
				Double_t sideLo, Double_t sideHi,
				Double_t longLo, Double_t longHi) {
  mQoutVal[0] = outLo; mQoutVal[1] = outHi;
  mQsideVal[0] = sideLo; mQsideVal[1] = sideHi;
  mQlongVal[0] = longLo; mQlongVal[1] = longHi;
}

//_________________
void HbtFitter::Set1DNonFemtoParameters(vector<double> par) {
  for(UInt_t i=0; i<par.size(); i++) {
    mNonFemto1DParArr[i] = par.at(i);
  }
}

//_________________
void HbtFitter::Set3DNonFemtoParameters(vector<double> parOut,
					  vector<double> parSide,
					  vector<double> parLong) {
  //Assuming the same amount of parameters 
  //for out, side and long components
  for(UInt_t i=0; i<parOut.size(); i++) {
    mNonFemto3DParArr[0][i] = parOut.at(i);
    mNonFemto3DParArr[1][i] = parSide.at(i);
    mNonFemto3DParArr[2][i] = parLong.at(i);
  }
}

//_________________
void HbtFitter::SetExcludeRange(Double_t qInvLo, Double_t qInvHi) {
  //Exlude range from the fit
  mExcludedVal[0] = qInvLo;
  mExcludedVal[1] = qInvHi;
}

//_________________
Double_t HbtFitter::Coulomb1D(Int_t qInvBin) {
  Double_t coulombVal;
  switch(mPartType) {
  case 3:
    coulombVal = PI_LONG[qInvBin-1];
    break;
  case 4:
    coulombVal = KK_LONG[qInvBin-1];
    break;
  default:
    std::cout << "Coulomb1D: Unknown particle type" << std::endl;
    coulombVal = 1.;
    break;
  }

  return coulombVal;
}

//_________________
Double_t HbtFitter::Coulomb3D(Int_t qOutBin, Int_t qSideBin, Int_t qLongBin) {

  Double_t coulombVal;

  if(!h3dCoulombWeight) { //If coulomb qInv weight function is not set

    std::cout << "Coulomb3D: The histogram with Coulomb weights is not set" << std::endl;
    coulombVal = 1.;
  }
  else { //If 3D qInv weight function is set

    //Obtain qInv weight from 3D
    Double_t qInvWeight = h3dCoulombWeight->GetBinContent(qOutBin, qSideBin, qLongBin);

    //Calculation of the Coulomb bin correction
    Int_t qBins = mQoutNbins;
    Double_t qHi = mQoutVal[1];
    Double_t qLow = h3dCoulombWeight->GetXaxis()->GetBinLowEdge(1);
    Double_t qStep = (qHi-qLow)/qBins;
    Int_t qCoulBin = (Int_t)(qInvWeight/qStep);

    //Obtain 1D Coulomb correction
  
    switch(mPartType) {
    case 3:
      coulombVal = PI_LONG[qCoulBin];
      break;
    case 4:
      coulombVal = KK_LONG[qCoulBin];
      break;
    default:
      std::cout << "Coulomb3D: Unknown particle type" << std::endl;
      coulombVal = 1.;
      break;
    }

    if(coulombVal<0.0001) {
      coulombVal=1.;
    }
    /*
      std::cout << Form("Coulomb3D: oBin=%2d sBin=%2d lBin=%2d cVal=%4.3f",
      qOutBin, qSideBin, qLongBin, coulombVal) << std::endl;
    */
  }
  return coulombVal;
}

//_________________
Double_t HbtFitter::TheorCorrFctn1D(Double_t *par, Int_t qInvBin) {

  Double_t mArg = 0;
  Double_t mTheory = 0;
  Double_t mNonFemto = 0;
  Double_t qInv = h1dNumer->GetBinCenter(qInvBin);
  Double_t qCoul;
  if(mUseCoulombCorrection==0) {
    qCoul = 1.;
  }
  else {
    qCoul = Coulomb1D(qInvBin);
  }

  switch(mSourceForm) {
  case 0: //Gaussian source
    mArg = qInv*qInv*par[1]*par[1];
    mArg /= -1.*hBar*hBar;
    break;
  case 1: //Exponential source
    mArg = qInv*par[1];
    mArg /= -1.*hBar*hBar;
    mArg /= R_GAUS_TO_EXPO;          //Scale from expo to gauss
    break;
  default:
    std::cout << "TheorCorrFctn1D: Unknown source form" << std::endl;
    mArg = qInv*qInv*par[1]*par[1];
    break;
  }
  
  mTheory = (1 - par[0]) + par[0]*qCoul*(1+TMath::Exp(mArg)) * 
    NonFemto1D(&par[2], qInv); //Ideal * NonFemto

  return mTheory;
}

//_________________
Double_t HbtFitter::TheorCorrFctn3D(Double_t *par, Int_t qOutBin, Int_t qSideBin, Int_t qLongBin) {

  Double_t mArg = 0.;
  Double_t mTheory = 0.;
  Double_t mNonFemto = 0.;
  Double_t qOut = h3dNumer->GetXaxis()->GetBinCenter(qOutBin);
  Double_t qSide = h3dNumer->GetYaxis()->GetBinCenter(qSideBin);
  Double_t qLong = h3dNumer->GetZaxis()->GetBinCenter(qLongBin);
  
  Double_t qCoul;
  if(mUseCoulombCorrection == 0) {
    qCoul = 1.;
  }
  else {
    qCoul = Coulomb3D(qOutBin, qSideBin, qLongBin);
  }

  switch(mSourceForm) {
  case 0: //Gaussian source
    mArg =  par[1]*par[1]*qOut*qOut;             //Out^2
    mArg += par[2]*par[2]*qSide*qSide;           //Side^2
    mArg += par[3]*par[3]*qLong*qLong;           //Long^2
    mArg += 2*par[4]*par[4]*qOut*qSide;            //Out*Side
    mArg += 2*par[5]*par[5]*qOut*qLong;            //Out*Long
    mArg += 2*par[6]*par[6]*qSide*qLong;           //Side*Long
    mArg /= -1.*hBar*hBar;
    break;
  case 1: //Exponential source. Note that the cross-terms should not present in fit
    mArg =  TMath::Abs(par[1]*par[1]*qOut*qOut);             //Out
    mArg += TMath::Abs(par[2]*par[2]*qSide*qSide);            //Side
    mArg += TMath::Abs(par[3]*par[3]*qLong*qLong);            //Long
    mArg += 2*par[4]*qOut*qSide;       //Out*Side
    mArg += 2*par[5]*qOut*qLong;       //Out*Long
    mArg += 2*par[6]*qSide*qLong;      //Side*Long
    mArg /= -1.*hBar*hBar;
    mArg /= R_GAUS_TO_EXPO;          //Scale from expo to gaus
    break;
  }

  mTheory = (1 - par[0] + par[0]*qCoul*(1+TMath::Exp(mArg))) * 
    NonFemto3D(&par[7], qOut, qSide, qLong);  //Ideal*NonFemto
  return mTheory;
}

//_________________
Double_t HbtFitter::NonFemto1D(Double_t *par, Double_t qInv) {

  Double_t nonFemtoVal = 0;
  switch(mNonFemtoForm) {
  case 0: //Constant
    nonFemtoVal = par[0];
    break;
  case 1: //Linear
    nonFemtoVal = par[0]*(1 + par[1]*qInv);
    break;
  case 2: //Second oreder polynomial
    nonFemtoVal = par[0]*(1 + par[1]*qInv + par[2]*qInv*qInv);
    break;
  case 3: //Sqrt(pol4)
    nonFemtoVal = par[0]*TMath::Sqrt(1 + par[1]*qInv*qInv + par[2]*qInv*qInv*qInv*qInv);
    break;
  case 4: //Gaussian-like
    nonFemtoVal = par[0]*(1 + TMath::Exp(-par[1]*qInv*qInv));
    break;
  case 5: //Hyperbolic
    nonFemtoVal = par[0] / (1 + par[1]*qInv*qInv);
    break;
  default:
    nonFemtoVal = par[0];
    break;
  }

  return nonFemtoVal;
}

//_________________
Double_t HbtFitter::NonFemto3D(Double_t *par, Double_t qOut, Double_t qSide, Double_t qLong) {

  Double_t nonFemtoVal = 0;
  Double_t mArg;
  switch(mNonFemtoForm) {
  case 0: //Constant
    nonFemtoVal = par[0];
    break;
  case 1: //Linear
    nonFemtoVal  = 1;
    nonFemtoVal += par[1]*qOut;
    nonFemtoVal += par[2]*qSide;
    nonFemtoVal += par[3]*qLong;
    nonFemtoVal *= par[0];
    break;
  case 2: //Second order polynomial
    nonFemtoVal  = 1;
    nonFemtoVal += par[1]*qOut;
    nonFemtoVal += par[2]*qSide;
    nonFemtoVal += par[3]*qLong;
    nonFemtoVal += par[4]*qOut*qOut;
    nonFemtoVal += par[5]*qSide*qSide;
    nonFemtoVal += par[6]*qLong*qLong;
    nonFemtoVal *= par[0];
    break;
  case 3: //Sqrt(Pol4)
    mArg = 1;
    mArg += par[1]*qOut*qOut;
    mArg += par[2]*qSide*qSide;
    mArg += par[3]*qLong*qLong;
    mArg += par[4]*qOut*qOut*qOut*qOut;
    mArg += par[5]*qSide*qSide*qSide*qSide;
    mArg += par[6]*qLong*qLong*qLong*qLong;
    nonFemtoVal = TMath::Sqrt(mArg);
    nonFemtoVal *= par[0];
    break;
  case 4: //Gaussian-like
    nonFemtoVal  = 1;
    nonFemtoVal += TMath::Exp(-1.*par[1]*qOut*qOut);
    nonFemtoVal += TMath::Exp(-1.*par[2]*qSide*qSide);
    nonFemtoVal += TMath::Exp(-1.*par[3]*qLong*qLong);
    nonFemtoVal *= par[0];
    break;
  case 5: //Hyperbolic
    nonFemtoVal = par[0];
    mArg  = 1;
    mArg += par[1]*qOut*qOut;
    mArg += par[2]*qSide*qSide;
    mArg += par[3]*qLong*qLong;
    nonFemtoVal /= mArg;
    break;
  default:
    nonFemtoVal = par[0];
    break;
  }

  return nonFemtoVal;
}

//_________________
Double_t HbtFitter::CorrFctn1DChiSqr(Int_t &npar, Double_t *gin, Double_t &f, 
				     Double_t *par, Int_t iflag) {

  f = 0.;

  Double_t mChisq = 0;
  Double_t mDelta;
  Double_t mQinvErr, A, B, C;
  Double_t mcA, mcB;
  Int_t mNQinvBins = mQinvNbins;

  Double_t mNonFemtoCorrFactor = 1.;
  TH1F hCorrectedOnMc("hCorrectedOnMc","hCorrectedOnMc",mQinvNbins, mQinvVal[0], mQinvVal[1]);
  if(mCorrectOnMc) {
    hCorrectedOnMc.Divide(h1dCorrFctn,h1dMcCorrFctn);
  }

  //Calculating overall chi2
  for(Int_t iQbin = mQinvFitRangeBin[1]; iQbin>=mQinvFitRangeBin[0]; iQbin--) {

    mQinvErr = h1dCorrFctn->GetBinError(iQbin);
    A = h1dNumer->GetBinContent(iQbin);
    B = h1dDenom->GetBinContent(iQbin);
    C = A/B;
    //If divide on MC CorrFctn
    if(mCorrectOnMc) {
      mNonFemtoCorrFactor = h1dMcCorrFctn->GetBinContent(iQbin);
      C /= mNonFemtoCorrFactor;
      mQinvErr = hCorrectedOnMc.GetBinError(iQbin);
    }

    if( (A>=0.0001) && (B>=0.0001)) {
      mDelta = (C - TheorCorrFctn1D(par, iQbin))/mQinvErr;
      mChisq += mDelta*mDelta;
    }
  }
  f = mChisq;

  return f;
}

//_________________
Double_t HbtFitter::CorrFctn1DLogLike(Int_t &npar, Double_t *gin, Double_t &f, 
				      Double_t *par, Int_t iflag) {

  f = 0.;
  Double_t mChisq = 0;
  Double_t mDelta;
  Double_t A, B, C, mcA, mcB;
  Int_t mNQinvBins = mQinvNbins;

  Double_t mNonFemtoCorrFactor;

  //Calculating of the log-likelihood chi2
  for(Int_t iQbin = mQinvFitRangeBin[1]; iQbin>=mQinvFitRangeBin[0]; iQbin--) {
    
    A = h1dNumer->GetBinContent(iQbin);
    B = h1dDenom->GetBinContent(iQbin);
    if(mCorrectOnMc) {
      mcA = h1dMcNumer->GetBinContent(iQbin);
      mcB = h1dMcDenom->GetBinContent(iQbin);
      A *= mcB;
      B *= mcA;
    }

    if( (A>=0.0001) && (B>=0.0001) ) {
      
      C = TheorCorrFctn1D(par, iQbin);
      mDelta = -2* ( A * TMath::Log( (C/A)*((A+B)/(C+1)) ) +
		     B * TMath::Log( (A+B)/(B*(C+1)) ) );
      mChisq += mDelta;
    }
  } //for(Int_t iQbin = 1; iQbin<=mNQinvBins; iQbin++)
  f = mChisq;

  return f;
}

//_________________
Double_t HbtFitter::CorrFctn3DChiSqr(Int_t &npar, Double_t *gin, Double_t &f, 
				     Double_t *par, Int_t iflag) {

  f = 0.;
  Double_t mChisq = 0;
  Double_t mDelta = 0;
  Double_t A, B, C, mQcorrErr, mcA, mcB, mNonFemtoCorrFactor;
  Int_t mNqOutBins = mQoutNbins;
  Int_t mNqSideBins = mQsideNbins;
  Int_t mNqLongBins = mQlongNbins;
  mNonFemtoCorrFactor = 1.;
  TH3F hCorrectedOnMc("hCorrectedOnMc","hCorrectedOnMc",
		      mQoutNbins,mQoutVal[0],mQoutVal[1],
		      mQsideNbins,mQsideVal[0],mQsideVal[1],
		      mQlongNbins,mQlongVal[0],mQlongVal[1]);
  if(mCorrectOnMc) {
    hCorrectedOnMc.Divide(h3dCorrFctn,h3dMcCorrFctn);
  }

  for(Int_t iOutBin=mQoutFitRangeBin[1]; iOutBin>=mQoutFitRangeBin[0]; iOutBin--) {
    for(Int_t iSideBin=mQsideFitRangeBin[1]; iSideBin>=mQsideFitRangeBin[0]; iSideBin--) {
      for(Int_t iLongBin=mQlongFitRangeBin[1]; iLongBin>=mQlongFitRangeBin[0]; iLongBin--) {
	
	mQcorrErr = h3dCorrFctn->GetBinError(iOutBin, iSideBin, iLongBin);
	if(mQcorrErr<=0.0) continue;
	A = h3dNumer->GetBinContent(iOutBin, iSideBin, iLongBin);
	B = h3dDenom->GetBinContent(iOutBin, iSideBin, iLongBin);
	C = A/B;
	//If divide on MC CorrFctn
	if(mCorrectOnMc) {
	  mNonFemtoCorrFactor = h3dMcCorrFctn->GetBinContent(iOutBin, iSideBin, iLongBin);
	  C /= mNonFemtoCorrFactor;
	  mQcorrErr = hCorrectedOnMc.GetBinError(iOutBin, iSideBin, iLongBin);
	}

	if( (A>=0.0001) && (B>=0.0001)) {
	  mDelta = (C - TheorCorrFctn3D(par, iOutBin, iSideBin, iLongBin))/mQcorrErr;
	  mChisq += mDelta*mDelta;
	}

      } //for(Int_t iLongBin=1; iLongBin<=mNqLongBins; iLongBin++)
    } //for(Int_t iSideBin=1; iSideBin<=mNqSideBins; iSideBin++)
  } //for(Int_t iOutBin=1; iOutBin<=mNqOutBins; iOutBin++)

  f = mChisq;
  return f;
}

//_________________
Double_t HbtFitter::CorrFctn3DLogLike(Int_t &npar, Double_t *gin, Double_t &f, 
				      Double_t *par, Int_t iflag) {

  f = 0.;
  Double_t mChisq = 0;
  Double_t mDelta = 0;
  Double_t A, B, C, mcA, mcB;
  Int_t mNqOutBins = mQoutNbins;
  Int_t mNqSideBins = mQsideNbins;
  Int_t mNqLongBins = mQlongNbins;

  //Calculation of the log-likelihood chi2
  for(Int_t iOutBin=mQoutFitRangeBin[1]; iOutBin>=mQoutFitRangeBin[0]; iOutBin--) {
    for(Int_t iSideBin=mQsideFitRangeBin[1]; iSideBin>=mQsideFitRangeBin[0]; iSideBin--) {
      for(Int_t iLongBin=mQlongFitRangeBin[1]; iLongBin>=mQlongFitRangeBin[0]; iLongBin--) {

	A = h3dNumer->GetBinContent(iOutBin, iSideBin, iLongBin);
	B = h3dDenom->GetBinContent(iOutBin, iSideBin, iLongBin);

	if( (A>=0.0001) && (B>=0.0001) ) {
	  
	  C = TheorCorrFctn3D(par, iOutBin, iSideBin, iLongBin);
	  if(mCorrectOnMc) {
	    mcA = h3dMcNumer->GetBinContent(iOutBin, iSideBin, iLongBin);
	    mcB = h3dMcDenom->GetBinContent(iOutBin, iSideBin, iLongBin);
	    C *= mcA;
	    C /= mcB;
	  }
	  mDelta = -2. * ( A * TMath::Log( (C/A) * ((A+B)/(C+1)) ) +
			   B * TMath::Log( (1.0/B) * ((A+B)/(C+1))) );
	  mChisq += mDelta;	  
	}

      } //for(Int_t iLongBin=1; iLongBin<=mNqLongBins; iLongBin++)
    } //for(Int_t iSideBin=1; iSideBin<=mNqSideBins; iSideBin++)
  } //for(Int_t iOutBin=1; iOutBin<=mNqOutBins; iOutBin++)

  f = mChisq;
  return f;
}

//_________________
void HbtFitter::FitFunc1DChiSquare(Int_t &npar, Double_t *gin, Double_t &f, 
				   Double_t *par, Int_t iflag) {

  f = objHbtFitter->CorrFctn1DChiSqr(npar, gin, f, par, iflag);
}

//_________________
void HbtFitter::FitFunc1DLogLike(Int_t &npar, Double_t *gin, Double_t &f, 
				 Double_t *par, Int_t iflag) {

  f = objHbtFitter->CorrFctn1DLogLike(npar, gin, f, par, iflag);
}

//_________________
void HbtFitter::FitFunc3DChiSquare(Int_t &npar, Double_t *gin, Double_t &f, 
				   Double_t *par, Int_t iflag) {
  
  f = objHbtFitter->CorrFctn3DChiSqr(npar, gin, f, par, iflag);
}

//_________________
void HbtFitter::FitFunc3DLogLike(Int_t &npar, Double_t *gin, Double_t &f, 
				 Double_t *par, Int_t iflag) {
  
  f = objHbtFitter->CorrFctn3DLogLike(npar, gin, f, par, iflag);
}

//_________________
void HbtFitter::MakeProjections() {
  
  std::cout << "Start making projections";
  //Experimental Q distribution projections
  MakeExpProjections();

  if(mCorrectOnMc) {
    //MC Q distribution projections
    MakeMcProjections();

    //Correct projections
    h3dCorrFctnOutProj->Divide(h3dMcCorrFctnOutProj);
    h3dCorrFctnSideProj->Divide(h3dMcCorrFctnSideProj);
    h3dCorrFctnLongProj->Divide(h3dMcCorrFctnLongProj);
  }

  //Calculate 3D fit projections
  MakeFitProjections();

  //Scale distributions
  Double_t mScale = 1./mNonFemto3DVal.at(0);
  h3dCorrFctn->Scale(mScale);
  h3dFitNumerator->Divide(h3dDenom);
  h3dFitNumerator->Scale(mScale);

  Calculate3DChi2();

  //Correlation function scale
  h3dCorrFctnOutProj->Scale(mScale);
  h3dCorrFctnSideProj->Scale(mScale);
  h3dCorrFctnLongProj->Scale(mScale);
  //Fit function scale
  h3dCorrFctnFitOutProj->Scale(mScale);
  h3dCorrFctnFitSideProj->Scale(mScale);
  h3dCorrFctnFitLongProj->Scale(mScale);
  
  std::cout << "\t[DONE]" << std::endl;
}

//_________________
void HbtFitter::Calculate3DChi2() {

  mChi2PerNDF = 0.;
  mChi2 = 0;
  double mNBinsNDF = 0;
  double mError;
  double mCorrFctnVal, mFitFctnVal, mDelta;
  for(Int_t iOutBin=mQoutFitRangeBin[1]; iOutBin>=mQoutFitRangeBin[0]; iOutBin--) {
    for(Int_t iSideBin=mQsideFitRangeBin[1]; iSideBin>=mQsideFitRangeBin[0]; iSideBin--) {
      for(Int_t iLongBin=mQlongFitRangeBin[1]; iLongBin>=mQlongFitRangeBin[0]; iLongBin--) {

	mError = h3dCorrFctn->GetBinError(iOutBin,iSideBin,iLongBin);
	if(mError<=0) continue;
	mNBinsNDF++;
	mCorrFctnVal = h3dCorrFctn->GetBinContent(iOutBin,iSideBin,iLongBin);
	mFitFctnVal = h3dFitNumerator->GetBinContent(iOutBin,iSideBin,iLongBin);
	mDelta = (mCorrFctnVal - mFitFctnVal)/mError;
	mChi2 += (mDelta*mDelta);
      }
    }
  }

  mNDF = mNBinsNDF - (mFemto3DParNum+mNonFemto3DParNum) + mFixParameters3D.size();
  mChi2PerNDF = mChi2/mNDF;
}

//_________________
void HbtFitter::MakeFitProjections() {
  
  //Make fit numerator
  MakeFitNumerator();
  
  if(!h3dFitQoutProj || !h3dQoutProjDen) {
    std::cout << "MakeFitProjections: Out projections are empty" << std::endl;
    std::cout << "Fit out proj numerator address: " << h3dFitQoutProj 
	      << " Fit out proj denominator address: " << h3dQoutProjDen << std::endl;
    return;
  }

  if(!h3dFitQsideProj || !h3dQsideProjDen) {
    std::cout << "MakeFitProjections: Side projections are empty" << std::endl;
    std::cout << "Fit side proj numerator address: " << h3dFitQsideProj 
	      << " Fit side proj denominator address: " << h3dQsideProjDen << std::endl;
    return;
  }

  if(!h3dFitQlongProj || !h3dQlongProjDen) {
    std::cout << "MakeFitProjections: Long projections are empty" << std::endl;
    std::cout << "Fit long proj numerator address: " << h3dFitQlongProj 
	      << " Fit long proj denominator address: " << h3dQlongProjDen << std::endl;
    return;
  }

  //TString name, title;

  //Out

  TH1F* hCorrOut = new TH1F(*h3dFitQoutProj);
  /*
  name=""; title="";
  name = hCorrOut->GetName();
  title = hCorrOut->GetTitle();
  name += "_fit";
  title += " fit";
  */
  hCorrOut->Divide(h3dQoutProjDen);
  for(Int_t iBin=1; iBin<=(hCorrOut->GetNbinsX()); iBin++) {
    if(hCorrOut->GetBinContent(iBin) != 0) {
      hCorrOut->SetBinError(iBin, 0.);
    }
  }
  h3dCorrFctnFitOutProj = hCorrOut;
  //h3dCorrFctnFitOutProj->Scale(1./mNonFemto3DVal.at(0));
  //Side
  TH1F* hCorrSide = new TH1F(*h3dFitQsideProj);
  /*
  name=""; title="";
  name = hCorrSide->GetName();
  title = hCorrSide->GetTitle();
  name += "_fit";
  title += " fit";
  */
  hCorrSide->Divide(h3dQsideProjDen);
  for(Int_t iBin=1; iBin<=(hCorrSide->GetNbinsX()); iBin++) {
    if(hCorrSide->GetBinContent(iBin) != 0) {
      hCorrSide->SetBinError(iBin, 0.);
    }
  }
  h3dCorrFctnFitSideProj = hCorrSide;
  //h3dCorrFctnFitSideProj->Scale(1./mNonFemto3DVal.at(0));
  //Long
  TH1F* hCorrLong = new TH1F(*h3dFitQlongProj);
  /*
  name=""; title="";
  name = hCorrLong->GetName();
  title = hCorrLong->GetTitle();
  name += "_fit";
  title += " fit";
  */
  hCorrLong->Divide(h3dQlongProjDen);
  for(Int_t iBin=1; iBin<=(hCorrLong->GetNbinsX()); iBin++) {
    if(hCorrLong->GetBinContent(iBin) != 0) {
      hCorrLong->SetBinError(iBin, 0.);
    }
  }
  h3dCorrFctnFitLongProj = hCorrLong;
  //h3dCorrFctnFitLongProj->Scale(1./mNonFemto3DVal.at(0));
}

//_________________
void HbtFitter::MakeFitNumerator() {
			     
  //Obtain fit parameters
  Double_t fitPar[mFemto3DParNum+mNonFemto3DParNum];
  for(Int_t iIter=0; iIter<mFemto3DParNum; iIter++) {
    fitPar[iIter] = mFemto3DVal.at(iIter);
  }
  for(Int_t iIter=0; iIter<mNonFemto3DParNum; iIter++) {
    fitPar[mFemto3DParNum+iIter] = mNonFemto3DVal.at(iIter);
  }

  h3dFitNumerator = new TH3F("h3dFitNumerator","h3dFitNumerator",
			     mQoutNbins, mQoutVal[0],mQoutVal[1],
			     mQsideNbins, mQsideVal[0],mQsideVal[1],
			     mQlongNbins, mQlongVal[0],mQlongVal[1]);

  Double_t mNonFemtoCorrFactor = 1.;
  Double_t wTheor, wFitNum;
  //Bin-by-bin loop
  /*
  for(Int_t iOutBin=1; iOutBin<=mQoutNbins; iOutBin++) {
    for(Int_t iSideBin=1; iSideBin<=mQsideNbins; iSideBin++) {
      for(Int_t iLongBin=1; iLongBin<=mQlongNbins; iLongBin++) {
  */
  for(Int_t iOutBin=mQoutFitRangeBin[0]; iOutBin<=mQoutFitRangeBin[1]; iOutBin++) {
    for(Int_t iSideBin=mQsideFitRangeBin[0]; iSideBin<=mQsideFitRangeBin[1]; iSideBin++) {
      for(Int_t iLongBin=mQlongFitRangeBin[0]; iLongBin<=mQlongFitRangeBin[1]; iLongBin++) {
	
	wTheor = TheorCorrFctn3D(fitPar, iOutBin, iSideBin, iLongBin);
	if(mCorrectOnMc) {
	  mNonFemtoCorrFactor = h3dMcCorrFctn->GetBinContent(iOutBin, iSideBin, iLongBin);
	  wTheor *= mNonFemtoCorrFactor;
	}
	wFitNum = wTheor * h3dDenom->GetBinContent(iOutBin, iSideBin, iLongBin);
	h3dFitNumerator->SetBinContent(iOutBin, iSideBin, iLongBin, wFitNum);
	h3dFitNumerator->SetBinError(iOutBin, iSideBin, iLongBin, 0.);

      } //for(Int_t iLongBin=1; iLongBin<=mQlongNbins; iLongBin++)
    } //for(Int_t iSideBin=1; iSideBin<=mQsideNbins, iSideBin++)
  } //for(Int_t iOutBin=1; iOutBin<=mQoutNbins, iOutBin++)


  if(mProjectionType == 0) { //By bins
    Project3DDistribution(h3dFitNumerator, 
			  h3dFitQoutProj, h3dFitQsideProj, h3dFitQlongProj,
			  mOutProjBins[0], mOutProjBins[1], 
			  mSideProjBins[0], mSideProjBins[1],
			  mLongProjBins[0], mLongProjBins[1]);
  }
  else { //By values
    Project3DDistribution(h3dFitNumerator, 
			  h3dFitQoutProj, h3dFitQsideProj, h3dFitQlongProj,
			  mOutProjVals[0], mOutProjVals[1], 
			  mSideProjVals[0], mSideProjVals[1],
			  mLongProjVals[0], mLongProjVals[1]);
  }
}

//_________________
void HbtFitter::MakeExpProjections() {
  
  if(!h3dNumer) {
    std::cout << "MakeExpProjections: 3D numerator is not set" << std::endl;
    return;
  }
  if(!h3dDenom) {
    std::cout << "MakeExpProjections: 3D denominator is not set" << std::endl;
    return;
  }

  //h3dNumer->Draw();
  if(mProjectionType == 0) { //By bins
    Project3DCorrFctn(h3dNumer, h3dDenom,
		      h3dQoutProjNum, h3dQsideProjNum, h3dQlongProjNum,
		      h3dQoutProjDen, h3dQsideProjDen, h3dQlongProjDen,
		      mOutProjBins[0], mOutProjBins[1], 
		      mSideProjBins[0], mSideProjBins[1],
		      mLongProjBins[0], mLongProjBins[1]);
  }
  else { //By values
    Project3DCorrFctn(h3dNumer, h3dDenom,
		      h3dQoutProjNum, h3dQsideProjNum, h3dQlongProjNum,
		      h3dQoutProjDen, h3dQsideProjDen, h3dQlongProjDen,
		      mOutProjVals[0], mOutProjVals[1], 
		      mSideProjVals[0], mSideProjVals[1],
		      mLongProjVals[0], mLongProjVals[1]);
  }

  if(!h3dQoutProjNum || !h3dQsideProjNum || !h3dQlongProjNum) {
    std::cout << "MakeExpProjections: Numerator projections are not defined" << std::endl;
    std::cout << "Out proj address: " << h3dQoutProjNum 
	      << " Side proj address: " << h3dQsideProjNum 
	      << " Long proj address: " << h3dQlongProjNum << std::endl;
    return;
  }

  if(!h3dQoutProjDen || !h3dQsideProjDen || !h3dQlongProjDen) {
    std::cout << "MakeExpProjections: Denominator projections are not defined" << std::endl;
    std::cout << "Out proj address: " << h3dQoutProjDen 
	      << " Side proj address: " << h3dQsideProjDen 
	      << " Long proj address: " << h3dQlongProjDen << std::endl;
    return;
  }

  //Out
  TH1F* hCorrOut = new TH1F(*h3dQoutProjNum);
  hCorrOut->Divide(h3dQoutProjDen);
  h3dCorrFctnOutProj = hCorrOut;
  //Side
  TH1F* hCorrSide = new TH1F(*h3dQsideProjNum);
  hCorrSide->Divide(h3dQsideProjDen);
  h3dCorrFctnSideProj = hCorrSide;
  //Long 
  TH1F* hCorrLong = new TH1F(*h3dQlongProjNum);
  hCorrLong->Divide(h3dQlongProjDen);
  h3dCorrFctnLongProj = hCorrLong;
}

//_________________
void HbtFitter::MakeMcProjections() {

  if(!h3dMcNumer) {
    std::cout << "MakeMcProjections: 3D numerator is not set" << std::endl;
    return;
  }
  if(!h3dMcDenom) {
    std::cout << "MakeMcProjections: 3D denominator is not set" << std::endl;
    return;
  }

  if(mProjectionType == 0) { //By bins
    Project3DCorrFctn(h3dMcNumer, h3dMcDenom,
		      h3dMcQoutProjNum, h3dMcQsideProjNum, h3dMcQlongProjNum,
		      h3dMcQoutProjDen, h3dMcQsideProjDen, h3dMcQlongProjDen,
		      mOutProjBins[0], mOutProjBins[1], 
		      mSideProjBins[0], mSideProjBins[1],
		      mLongProjBins[0], mLongProjBins[1]);
  }
  else { //By values
  Project3DCorrFctn(h3dMcNumer, h3dMcDenom,
		    h3dMcQoutProjNum, h3dMcQsideProjNum, h3dMcQlongProjNum,
		    h3dMcQoutProjDen, h3dMcQsideProjDen, h3dMcQlongProjDen,
		    mOutProjVals[0], mOutProjVals[1], 
		    mSideProjVals[0], mSideProjVals[1],
		    mLongProjVals[0], mLongProjVals[1]);
  }

  if(!h3dMcQoutProjNum || !h3dMcQsideProjNum || !h3dMcQlongProjNum) {
    std::cout << "MakeMcProjections: Numerator projections are not defined" << std::endl;
    std::cout << "Out proj address: " << h3dMcQoutProjNum 
	      << " Side proj address: " << h3dMcQsideProjNum 
	      << " Long proj address: " << h3dMcQlongProjNum << std::endl;
    return;
  }

  if(!h3dMcQoutProjDen || !h3dMcQsideProjDen || !h3dMcQlongProjDen) {
    std::cout << "MakeMcProjections: Denominator projections are not defined" << std::endl;
    std::cout << "Out proj address: " << h3dMcQoutProjDen 
	      << " Side proj address: " << h3dMcQsideProjDen 
	      << " Long proj address: " << h3dMcQlongProjDen << std::endl;
    return;
  }

  //Out
  TH1F* hCorrOut = new TH1F(*h3dMcQoutProjNum);
  hCorrOut->Divide(h3dMcQoutProjDen);
  h3dMcCorrFctnOutProj = hCorrOut;
  //Side
  TH1F* hCorrSide = new TH1F(*h3dMcQsideProjNum);
  hCorrSide->Divide(h3dMcQsideProjDen);
  h3dMcCorrFctnSideProj = hCorrSide;
  //Long 
  TH1F* hCorrLong = new TH1F(*h3dMcQlongProjNum);
  hCorrLong->Divide(h3dMcQlongProjDen);
  h3dMcCorrFctnLongProj = hCorrLong; 
}

//_________________
void HbtFitter::Project3DCorrFctn(TH3F* hNumer, TH3F* hDenom, 
				  TH1F *&hOutProjNum, TH1F *&hSideProjNum, 
				  TH1F *&hLongProjNum,
				  TH1F *&hOutProjDen, TH1F *&hSideProjDen, 
				  TH1F *&hLongProjDen,
				  Int_t mOutBinLow, Int_t mOutBinHi,
				  Int_t mSideBinLow, Int_t mSideBinHi,
				  Int_t mLongBinLow, Int_t mLongBinHi) {
  
  Project3DDistribution(hNumer, hOutProjNum, hSideProjNum, hLongProjNum,
			mOutBinLow, mOutBinHi, 
			mSideBinLow, mSideBinHi,
			mLongBinLow, mLongBinHi);

  if(!hOutProjNum || !hSideProjNum || !hLongProjNum) {
    std::cout << "Project3DCorrFctn: One of the numerator projections is not set" << std::endl;
  }
  
  Project3DDistribution(hDenom, hOutProjDen, hSideProjDen, hLongProjDen,
			mOutBinLow, mOutBinHi, 
			mSideBinLow, mSideBinHi,
			mLongBinLow, mLongBinHi);

  if(!hOutProjDen || !hSideProjDen || !hLongProjDen) {
    std::cout << "Project3DCorrFctn: One of the denominator projections is not set" << std::endl;
  }
}

//_________________
void HbtFitter::Project3DCorrFctn(TH3F* hNumer, TH3F* hDenom, 
				  TH1F *&hOutProjNum, TH1F *&hSideProjNum, 
				  TH1F *&hLongProjNum,
				  TH1F *&hOutProjDen, TH1F *&hSideProjDen, 
				  TH1F *&hLongProjDen,
				  Double_t mQoutLow, Double_t mQoutHi,
				  Double_t mQsideLow, Double_t mQsideHi,
				  Double_t mQlongLow, Double_t mQlongHi) {
  
  Project3DDistribution(hNumer, 
			hOutProjNum, hSideProjNum, hLongProjNum,
			mQoutLow, mQoutHi, 
			mQsideLow, mQsideHi,
			mQlongLow, mQlongHi);
  if(!hOutProjNum || !hSideProjNum || !hLongProjNum) {
    std::cout << "Project3DCorrFctn: One of the numerator projections is not set" << std::endl;
  }

  Project3DDistribution(hDenom, 
			hOutProjDen, hSideProjDen, hLongProjDen,
			mQoutLow, mQoutHi, 
			mQsideLow, mQsideHi,
			mQlongLow, mQlongHi);
  if(!hOutProjDen || !hSideProjDen || !hLongProjDen) {
    std::cout << "Project3DCorrFctn: One of the denominator projections is not set" << std::endl;
  }
}

//_________________
void HbtFitter::Project3DDistribution(TH3F* hQdist, 
				      TH1F *&hQoutProj, TH1F *&hQsideProj, TH1F *&hQlongProj,
				      Int_t mOutBinLow, Int_t mOutBinHi,
				      Int_t mSideBinLow, Int_t mSideBinHi,
				      Int_t mLongBinLow, Int_t mLongBinHi) {

  if(!hQdist) {
    std::cout << "Project3DDistribution: 3D histogram is not set" << std::endl;
    return;
  }

  //hQdist->Draw();

  //Obtaining number of bins in each axis
  Int_t mNbinsX = hQdist->GetNbinsX();
  Int_t mNbinsY = hQdist->GetNbinsY();
  Int_t mNbinsZ = hQdist->GetNbinsZ();  

  if(!mNbinsX || !mNbinsY || !mNbinsZ) {
    std::cout << "Project3DDistribution: 3D distribution is empty" << std::endl;
    return;
  }

  TH1F* hXproj;
  TH1F* hYproj;
  TH1F* hZproj;

  TString name, title;
  //Project 3D histogram using number of bins
  //Out
  name = hQdist->GetName();
  title = hQdist->GetTitle();
  name += "_out";
  title += " out projection";
  hXproj = (TH1F*)hQdist->ProjectionX("hQoutProj",
				      mSideBinLow, mSideBinHi,
				      mLongBinLow, mLongBinHi);
  hXproj->SetNameTitle(name.Data(), title.Data());

  //Side
  name = hQdist->GetName();
  title = hQdist->GetTitle();
  name += "_side";
  title += " side projection";
  hYproj = (TH1F*)hQdist->ProjectionY("hQsideProj",
				      mOutBinLow, mOutBinHi,
				      mLongBinLow, mLongBinHi);
  hYproj->SetNameTitle(name.Data(), title.Data());

  //Long
  name = hQdist->GetName();
  title = hQdist->GetTitle();
  name += "_long";
  title += " long projection";
  hZproj = (TH1F*)hQdist->ProjectionZ("hQlongProj",
				      mOutBinLow, mOutBinHi,
				      mSideBinLow, mSideBinHi);
  hZproj->SetNameTitle(name.Data(), title.Data());

  hQoutProj = hXproj;
  hQsideProj = hYproj;
  hQlongProj = hZproj;
}

//_________________
void HbtFitter::Project3DDistribution(TH3F* hQdist, 
				      TH1F *&hQoutProj, TH1F *&hQsideProj, TH1F *&hQlongProj,
				      Double_t mQoutLow, Double_t mQoutHi,
				      Double_t mQsideLow, Double_t mQsideHi,
				      Double_t mQlongLow, Double_t mQlongHi) {

  if(!hQdist) {
    std::cout << "Project3DDistribution: 3D histogram is not set" << std::endl;
    return;
  }

  //hQdist->Draw();

  //Obtaining number of bins in each axis
  Int_t mNbinsX = hQdist->GetNbinsX();
  Int_t mNbinsY = hQdist->GetNbinsY();
  Int_t mNbinsZ = hQdist->GetNbinsZ();  

  if(!mNbinsX || !mNbinsY || !mNbinsZ) {
    std::cout << "Project3DDistribution: 3D distribution has no bins" << std::endl;
    return;
  }

  //Obtaining the ranges
  Double_t xLow = hQdist->GetXaxis()->GetXmin();
  Double_t xHi = hQdist->GetXaxis()->GetXmax();
  Double_t yLow = hQdist->GetYaxis()->GetXmin();
  Double_t yHi = hQdist->GetYaxis()->GetXmax();
  Double_t zLow = hQdist->GetZaxis()->GetXmin();
  Double_t zHi = hQdist->GetZaxis()->GetXmax();

  TH1F* hXproj;
  TH1F* hYproj;
  TH1F* hZproj;

  TString name, title;

  //Project 3D histogram using the ranges
  //Out
  name = hQdist->GetName();
  title = hQdist->GetTitle();
  name += "_out";
  title += " out projection";
  hQdist->SetAxisRange(xLow, xHi, "x");
  hQdist->SetAxisRange(mQsideLow, mQsideHi, "y");
  hQdist->SetAxisRange(mQlongLow, mQlongHi, "z");
  hXproj = (TH1F*)hQdist->Project3D("xe");
  hXproj->SetNameTitle(name.Data(), title.Data());;

  //Side
  name = hQdist->GetName();
  title = hQdist->GetTitle();
  name += "_side";
  title += " side projection";
  hQdist->SetAxisRange(mQoutLow, mQoutHi, "x");
  hQdist->SetAxisRange(yLow, yHi, "y");
  hQdist->SetAxisRange(mQlongLow, mQlongHi, "z");
  hYproj = (TH1F*)hQdist->Project3D("ye");
  hYproj->SetNameTitle(name.Data(), title.Data());

  //Long
  name = hQdist->GetName();
  title = hQdist->GetTitle();
  name += "_long";
  title += " long projection";
  hQdist->SetAxisRange(mQoutLow, mQoutHi, "x");
  hQdist->SetAxisRange(mQsideLow, mQsideHi, "y");
  hQdist->SetAxisRange(zLow, zHi, "z");
  hZproj = (TH1F*)hQdist->Project3D("ze");
  hZproj->SetNameTitle(name.Data(), title.Data());
  
  hQoutProj = hXproj;
  hQsideProj = hYproj;
  hQlongProj = hZproj;
}