#include "DataProcessing.h"
#include "FemtoUtils.h"
#include "HbtFitter.h"
#include <TLegend.h>
#include <TCanvas.h>
#include <TMath.h>

const Double_t MHC     = 0.197326960277;  // GeV * fm
const Double_t MHC_SQR = 0.038937929230;  // GeV^2 * fm^2
const Double_t INVSQRTPI = 0.564189584;   // 1/sqrt(pi)
const Double_t PION_MASS = 0.13957;
const Double_t PION_MASS_SQR = 0.0194797;
const Double_t KAON_MASS = 0.493677;
const Double_t KAON_MASS_SQR = 0.24371698;

vector<double> DoFit(TH1* histo, Int_t PartType=4,
		     Int_t SourceFunc=0, Int_t NonFemtoFunc=1,
		     Bool_t DrawFitLegend=false, Bool_t ExcludePoints=false);
//_________________
DataProcessing::DataProcessing() {

  mExpFileName = "~/work/star/data/expKaonHBT_auau200_2011.root";
  mOutFileName = "~/work/star/soft/star_software/processing/auau/oKaonHBT_pp200.root";

  mQinvBins = 40;
  mQinvLow = 0.;
  mQinvHi = 0.4;

  mExpFile = NULL;
  mOutFile = NULL;

  mDetSelection = 2;

  mPartType = 4;

  mCorrFctnColor = 1;
  mCorrFctnLineWidth = 2;
  mCorrFctnMarkerStyle = 20;
  mCorrFctnMarkerSize = 1.5;

  mYaxisLo = 0.6;
  mYaxisHi = 1.6;

  mMultBins = 9;
  mKtBins = 24;
  mKtRanges = 6;

  mSourceFunc = 0;
  mNonFemtoFunc = 0;
  mDrawLegend = true;
}

//_________________
DataProcessing::~DataProcessing() {
  /* noting to do */
}

//_________________
void DataProcessing::SetCorrFctnStyle(Int_t color, Int_t width, 
				      Int_t style, Double_t size) {
  mCorrFctnColor = color;
  mCorrFctnLineWidth = width;
  mCorrFctnMarkerStyle = style;
  mCorrFctnMarkerSize = size;
}

//_________________
void DataProcessing::SetHistoStyle(TH1 *h) {
  h->SetLineColor(mCorrFctnColor);
  h->SetLineWidth(mCorrFctnLineWidth);
  h->SetMarkerStyle(mCorrFctnMarkerStyle);
  h->SetMarkerSize(mCorrFctnMarkerSize);
  h->SetMarkerColor(mCorrFctnColor);
  h->GetYaxis()->SetRangeUser(mYaxisLo,mYaxisHi);
}

//_________________
void DataProcessing::OpenFiles() {

  std::cout << "Opening file with experimental data: " << mExpFileName;
  mExpFile = TFile::Open(mExpFileName);
  std::cout << "\t [DONE]" << std::endl;
  std::cout << "Creating output file: " << mOutFileName;
  mOutFile = new TFile(mOutFileName,"recreate");
  std::cout << "\t [DONE]" << std::endl;
}

//_________________
void DataProcessing::DoProcessing() {

  std::cout << "Start processing..." << std::endl;
  OpenFiles();
  ReadHistos();
  ProcessHistos();
  FitHistos();
  MakeTGraphs();
  Finish();
  std::cout << "Processing is finished" << std::endl;
}

//_________________
void DataProcessing::ReadHistos() {

  std::cout << "Retrieving histograms from files...";

  for(Int_t iCharge=0; iCharge<2; iCharge++) {
    //for(Int_t iDet=0; iDet<5; iDet++) {
    for(Int_t iMult=0; iMult<mMultBins; iMult++) {
      for(Int_t iKt=0; iKt<mKtBins; iKt++) {

	if(mPartType == 4) {
	  //Kaons
	  TH1F *hExpKaonKtMixNum =
	    (TH1F*)mExpFile->Get(Form("hKaonQinvMixKt_%d_%d_%d_Num_bin_%d",
				      iCharge,mDetSelection,iMult,iKt));
	  TH1F *hExpKaonKtMixDen =
	    (TH1F*)mExpFile->Get(Form("hKaonQinvMixKt_%d_%d_%d_Den_bin_%d",
				      iCharge,mDetSelection,iMult,iKt));
	  hQinvMixNum.push_back(hExpKaonKtMixNum);
	  hQinvMixDen.push_back(hExpKaonKtMixDen);
	}
	else if(mPartType == 3) {
	  //Pions
	  TH1F *hExpPionKtMixNum =
	    (TH1F*)mExpFile->Get(Form("hPionQinvMixKt_%d_%d_%d_Num_bin_%d",
				      iCharge,mDetSelection,iMult,iKt));
	  TH1F *hExpPionKtMixDen =
	    (TH1F*)mExpFile->Get(Form("hPionQinvMixKt_%d_%d_%d_Den_bin_%d",
				      iCharge,mDetSelection,iMult,iKt));
	  hQinvMixNum.push_back(hExpPionKtMixNum);
	  hQinvMixDen.push_back(hExpPionKtMixDen);
	}

      } //for(Int_t iKt=0; iKt<20; iKt++)
    } //for(Int_t iMult=0; iMult<10; iMult++)
      //} //for(Int_t iDet=0; iDet<5; iDet++)
  } //for(Int_t iCharge=0; iCharge<2; iCharge++)

  std::cout << "\t [DONE]" << std::endl;

  std::cout << Form("Number of histos in numerator:   %d", hQinvMixNum.size()) 
	    << std::endl
	    << Form("Number of histos in denominator: %d", hQinvMixDen.size()) 
	    << std::endl;
}

//_________________
void DataProcessing::ProcessHistos() {

  std::cout << "Processing multiplicity histograms...";
  ProcessMultHistos();
  std::cout << "\t [DONE]" << std::endl;
  std::cout << "Processing Mult vs kT histograms...";
  ProcessMultKtHistos();
  std::cout << "\t [DONE]" << std::endl;
}

//_________________
TH1F* DataProcessing::BuildMultCF(Int_t charge, Int_t BinLo, Int_t BinHi) {

  if(charge<0 || charge>2) {
    std::cout << "DataProcessing::BuildMultCF [WARNING] - Wrong charge: " << charge
	      << std::endl;
    return 0;
  }
  
  if(BinLo<0 || BinHi>=mMultBins) {
    std::cout << "DataProcessing::BuildMultCF [WARNING] - Wrong multiplicity bins: " 
	      << BinLo << "_" << BinHi << std::endl;
    return 0;
  }

  TH1F *hNumSum = new TH1F("hNumSum","",mQinvBins, mQinvLow, mQinvHi);
  TH1F *hDenSum = new TH1F("hDenSum","",mQinvBins, mQinvLow, mQinvHi);
  TH1F *hRatio;
  TString sName = "h";

  sName += "CentCF_";
  sName += charge;
  sName += "_";
  sName += BinLo;
  sName += "_";
  sName += BinHi;

  hRatio = new TH1F(sName.Data(),"",mQinvBins,mQinvLow,mQinvHi);

  if(charge<2) {
    for(Int_t iMult=BinLo; iMult<=BinHi; iMult++) {
      for(Int_t iKt=0; iKt<mKtBins; iKt++) {
	hNumSum->Add(hQinvMixNum[charge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
	hDenSum->Add(hQinvMixDen[charge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
      } //for(Int_t iKt=0; iKt<mKtBins; iKt++)
      std::cout << std::endl;
    } //for(Int_t iMult=BinLo; iMult<=BinHi; iMult++)
  } //if(charge<2)
  else {
    for(Int_t iCharge=0; iCharge<2; iCharge++) {
      for(Int_t iMult=BinLo; iMult<=BinHi; iMult++) {
	for(Int_t iKt=0; iKt<mKtBins; iKt++) {

	  hNumSum->Add(hQinvMixNum[iCharge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
	  hDenSum->Add(hQinvMixDen[iCharge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
	} //for(Int_t iKt=0; iKt<mKtBins; iKt++)
      } //for(Int_t iMult=BinLo; iMult<=BinHi; iMult++)
    } //for(Int_t iCharge=0; iCharge<2; iCharge++)
  } //if(charge==2)

  Make1DCorrFctn(hRatio,hNumSum,hDenSum);
  delete hNumSum;
  delete hDenSum;

  return hRatio;
}

//_________________
TH1F* DataProcessing::BuildMultKtCF(Int_t charge, Int_t MultBinLo, Int_t MultBinHi,
				    Int_t KtBinLo, Int_t KtBinHi) {

  if(charge<0 || charge>2) {
    std::cout << "DataProcessing::BuildMultKtCF [WARNING] - Wrong charge: " << charge
	      << std::endl;
    return 0;
  }
  
  if(MultBinLo<0 || MultBinHi>=mMultBins) {
    std::cout << "DataProcessing::BuildMultKtCF [WARNING] - Wrong multiplicity bins: " 
	      << MultBinLo << "_" << MultBinHi << std::endl;
    return 0;
  }

  if(KtBinLo<0 || KtBinHi>=mKtBins) {
    std::cout << "DataProcessing::BuildMultKtCF [WARNING] - Wrong kT bins: " 
	      << KtBinLo << "_" << KtBinHi << std::endl;
    return 0;
  }


  TH1F *hNumSum = new TH1F("hNumSum","",mQinvBins, mQinvLow, mQinvHi);
  TH1F *hDenSum = new TH1F("hDenSum","",mQinvBins, mQinvLow, mQinvHi);
  TH1F *hRatio;
  TString sName = "h";

  sName += "MultKtCF_";
  sName += charge;
  sName += "_";
  sName += MultBinLo;
  sName += "_";
  sName += MultBinHi;
  sName += "_";
  sName += KtBinLo;
  sName += "_";
  sName += KtBinHi;

  hRatio = new TH1F(sName.Data(),"",mQinvBins,mQinvLow,mQinvHi);

  if(charge<2) {
    for(Int_t iMult=MultBinLo; iMult<=MultBinHi; iMult++) {
      for(Int_t iKt=KtBinLo; iKt<=KtBinHi; iKt++) {

	hNumSum->Add(hQinvMixNum[charge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
	hDenSum->Add(hQinvMixDen[charge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
      } //for(Int_t iKt=KtBinLo; iKt<=KtBinHi; iKt++)
    } //for(Int_t iMult=MultBinLo; iMult<=MultBinHi; iMult++)
  } //if(charge<2)
  else {
    for(Int_t iCharge=0; iCharge<2; iCharge++) {
      for(Int_t iMult=MultBinLo; iMult<=MultBinHi; iMult++) {
	for(Int_t iKt=KtBinLo; iKt<=KtBinHi; iKt++) {

	  hNumSum->Add(hQinvMixNum[iCharge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
	  hDenSum->Add(hQinvMixDen[iCharge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
	} //for(Int_t iKt=KtBinLo; iKt<=KtBinHi; iKt++)
      } //for(Int_t iMult=MultBinLo; iMult<=MultBinHi; iMult++)
    } //for(Int_t iCharge=0; iCharge<2; iCharge++)
  } //if(charge==2)

  Make1DCorrFctn(hRatio,hNumSum,hDenSum);
  delete hNumSum;
  delete hDenSum;

  return hRatio;
}

//_________________
void DataProcessing::ProcessMultHistos() {

  Int_t multBinLo, multBinHi;
  for(Int_t iCharge=0; iCharge<3; iCharge++) {
    for(Int_t iMultBin=0; iMultBin<mMultBins; iMultBin++) {

      std::cout << Form("ProcessMultHistos - iCharge: %d iMultBin: %d ", 
			iCharge, iMultBin) << std::endl;  
      multBinLo = iMultBin; multBinHi = iMultBin;
      hMultMixCF.push_back( BuildMultCF(iCharge, multBinLo, multBinHi) );
      SetHistoStyle( hMultMixCF.back() );

    } //for(Int_t iMultBin=0; iMultBin<mMultBins; iMultBin++)
  } //for(Int_t iCharge=0; iCharge<3; iCharge++)
}

//_________________
void DataProcessing::ProcessMultKtHistos() {

  Int_t multBinLo, multBinHi;
  Int_t ktBinLo, ktBinHi;
  for(Int_t iCharge=0; iCharge<3; iCharge++) {
    for(Int_t iMultBin=0; iMultBin<mMultBins; iMultBin++) {
      for(Int_t iKtBin=0; iKtBin<mKtRanges; iKtBin++) {

      multBinLo = iMultBin; multBinHi = iMultBin;
      switch(iKtBin) {
      case 0: ktBinLo = 0; ktBinHi = 5; break;
      case 1: ktBinLo = 6; ktBinHi = 7; break;
      case 2: ktBinLo = 8; ktBinHi = 9; break;
	//case 3: ktBinLo = 10; ktBinHi = 18; break;
      case 3: ktBinLo = 10; ktBinHi = 13; break;
      case 4: ktBinLo = 14; ktBinHi = 17; break;
      case 5: ktBinLo = 18; ktBinHi = 23; break;
      default: ktBinLo = 0; ktBinHi = 23; break;
      };
      hMultKtMixCF.push_back( BuildMultKtCF(iCharge, multBinLo, multBinHi,
					    ktBinLo, ktBinHi) );
      SetHistoStyle( hMultKtMixCF.back() );

      } //for(Int_t iKtBin=0; iKtBin<mKtRanges; iKtBin++)
    } //for(Int_t iMultBin=0; iMultBin<mMultBins; iMultBin++)
  } //for(Int_t iCharge=0; iCharge<3; iCharge++)
}

//_________________
void DataProcessing::Finish() {

  for(Int_t iCharge=0; iCharge<3; iCharge++) {
    gMultDepRadiiGrErr[iCharge]->Write();
    for(Int_t iCent=0; iCent<9; iCent++) {
      gMultKtDepRadiiGrErr[iCharge][iCent]->Write();
      gMultKtDepRadiiMtGrErr[iCharge][iCent]->Write();
    }
  }
  mOutFile->Write();
  mOutFile->Close();
  mExpFile->Close();
}

//_________________
void DataProcessing::FitHistos() {

  vector<double> mFitParams;

  for(unsigned int i=0; i<hMultMixCF.size(); i++) {
    if(!mFitParams.empty()) {
      mFitParams.clear();
    }
    mFitParams = DoFit( hMultMixCF[i], mPartType, mSourceFunc, 
			mNonFemtoFunc, mDrawLegend );
    mMultMixRadii.push_back(mFitParams.at(0));
    mMultMixRadiiErr.push_back(mFitParams.at(1));
    mMultMixLambda.push_back(mFitParams.at(3));
    mMultMixLambdaErr.push_back(mFitParams.at(4));
  } //for(unsigned int i=0; i<hMultMixCF.size(); i++)

  for(unsigned int i=0; i<hMultKtMixCF.size(); i++) {
    if(!mFitParams.empty()) {
      mFitParams.clear();
    }
    mFitParams = DoFit( hMultKtMixCF[i], mPartType, mSourceFunc, 
			mNonFemtoFunc, mDrawLegend );
    mMultKtMixRadii.push_back(mFitParams.at(0));
    mMultKtMixRadiiErr.push_back(mFitParams.at(1));
    mMultKtMixLambda.push_back(mFitParams.at(3));
    mMultKtMixLambdaErr.push_back(mFitParams.at(4));
  } //for(unsigned int i=0; i<hMultKtMixCF.size(); i++)
}

//_________________
void DataProcessing::MakeTGraphs() {

  Double_t mMultRadiiArr[3][9];
  Double_t mMultRadiiErrArr[3][9];
  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 mRadiiTmp[9];
  Double_t mRadiiErrTmp[9];

  //Double_t mKtVal[4] = {0.21, 0.35, 0.45, 0.66};
  //Double_t mKtErrVal[4] = {0.01, 0.01, 0.01, 0.01};
  //Double_t mKtRadiiTmp[4];
  //Double_t mKtRadiiErrTmp[4];
  //Double_t mMtVal[4];
  Double_t mKtRadiiTmp[6];
  Double_t mKtRadiiErrTmp[6];
  Double_t mKtVal[6] = {0.18, 0.31, 0.41, 0.55, 0.74, 0.97};
  Double_t mKtErrVal[6] = {0.01, 0.01, 0.01, 0.01, 0.01, 0.01};
  Double_t mMtVal[6];
  if(mPartType == 3) {
    for(Int_t i=0; i<mKtRanges; i++) {
      mMtVal[i] = TMath::Sqrt(PION_MASS_SQR + mKtVal[i]);
    }
  }
  else if(mPartType == 4) {
    for(Int_t i=0; i<mKtRanges; i++) {
      mMtVal[i] = TMath::Sqrt(KAON_MASS_SQR + mKtVal[i]);
    }   
  }

  Int_t mMultDepColor[3] = {2, 4, 1};
  Int_t mMultDepMarkerStyle;
  Double_t mMultDepMarkerSize = 1.6;
  if(mPartType == 4) {
    mMultDepMarkerStyle = 20;
  }
  else {
    mMultDepMarkerStyle = 21;
  }

  Int_t mMultKtDepColor[7] = {2, 4, 6, 8, 9, 28, 30};
  Double_t mMultKtDepMarkerSize = 1.5;
  Int_t mMultKtDepMarkerStyle[3];
  if(mPartType == 4) {
    mMultKtDepMarkerStyle[0] = 20;
    mMultKtDepMarkerStyle[1] = 24;
    mMultKtDepMarkerStyle[2] = 34;
  }
  else if(mPartType == 3) {
    mMultKtDepMarkerStyle[0] = 21;
    mMultKtDepMarkerStyle[1] = 25;
    mMultKtDepMarkerStyle[2] = 28;    
  }
   
  //Multiplicity dependence
  for(Int_t iCharge=0; iCharge<3; iCharge++) {
    if(iCharge==2) {
      std::cout << std::endl;
      std::cout << "Radii values (centrality) for sum of charges:" << std::endl;
    }
    for(Int_t iCent=0; iCent<mMultBins; iCent++) {
      if(iCent==9) continue; //Skip integrated values
      
      mMultRadiiArr[iCharge][iCent] = mMultMixRadii.at(iCharge*mMultBins+iCent);
      mRadiiTmp[iCent] = mMultMixRadii.at(iCharge*mMultBins+iCent);
      mMultRadiiErrArr[iCharge][iCent] = mMultMixRadiiErr.at(iCharge*mMultBins+iCent);
      mRadiiErrTmp[iCent] = mMultMixRadiiErr.at(iCharge*mMultBins+iCent);
     
      if(iCharge==2) {
	std::cout << " " << mMultRadiiArr[iCharge][iCent];
      }
    } //for(Int_t iCent=0; iCent<mMultBins; iCent++)

    gMultDepRadiiGrErr[iCharge] = new TGraphErrors(mMultBins,mCentrality,mRadiiTmp,
						   mCentralityErr,mRadiiErrTmp);
    gMultDepRadiiGrErr[iCharge]->SetName(Form("gMultDepRadiiGrErr_%d",iCharge));
    gMultDepRadiiGrErr[iCharge]->SetMarkerColor(mMultDepColor[iCharge]);
    gMultDepRadiiGrErr[iCharge]->SetLineColor(mMultDepColor[iCharge]);
    gMultDepRadiiGrErr[iCharge]->SetMarkerStyle(mMultDepMarkerStyle);
    gMultDepRadiiGrErr[iCharge]->SetMarkerSize(mMultDepMarkerSize);
  } //for(Int_t iCharge=0; iCharge<3; iCharge++)

  //Multiplicity and kT depence
  for(Int_t iCharge=0; iCharge<3; iCharge++) {
    if(iCharge==2) {
      std::cout << std::endl;
      std::cout << "Radii values (kT) for sum of charges:" << std::endl;
    }
    for(Int_t iCent=0; iCent<mMultBins; iCent++) {
      if(iCharge==2) {
	std::cout << std::endl;
	std::cout << "Centrality bin: " << iCent << " "; 
      }
      if(iCent==9) continue; //Skip integrated values
      for(Int_t iKt=0; iKt<mKtRanges; iKt++) {

	mKtRadiiTmp[iKt] = mMultKtMixRadii.at(iCharge*mMultBins*mKtRanges +
					      iCent*mKtRanges + iKt);
	mKtRadiiErrTmp[iKt] = mMultKtMixRadiiErr.at(iCharge*mMultBins*mKtRanges +
						    iCent*mKtRanges + iKt);

	if(iCharge==2) {
	  std::cout << " " << mKtRadiiTmp[iKt];
	}
      } //for(Int_t iKt=0; iKt<mKtRanges; iKt++)

      gMultKtDepRadiiGrErr[iCharge][iCent] = 
	new TGraphErrors(mKtRanges, mKtVal, mKtRadiiTmp, mKtErrVal, mKtRadiiErrTmp);
      gMultKtDepRadiiGrErr[iCharge][iCent]->SetName(Form("gMultKtDepRadiiGrErr_%d_%d",
							 iCharge,iCent));
      gMultKtDepRadiiMtGrErr[iCharge][iCent] =
	new TGraphErrors(mKtRanges, mMtVal, mKtRadiiTmp, mKtErrVal, mKtRadiiErrTmp);
      gMultKtDepRadiiMtGrErr[iCharge][iCent]->SetName(Form("gMultKtDepRadiiMtGrErr_%d_%d",
							   iCharge,iCent));
      if(iCent<2) continue; //Skip first 2 point due to the mistake.
      gMultKtDepRadiiGrErr[iCharge][iCent]->SetMarkerColor(mMultKtDepColor[iCent]);
      gMultKtDepRadiiGrErr[iCharge][iCent]->SetLineColor(mMultKtDepColor[iCent]);
      gMultKtDepRadiiGrErr[iCharge][iCent]->SetMarkerSize(mMultKtDepMarkerSize);
      gMultKtDepRadiiGrErr[iCharge][iCent]->SetMarkerStyle(mMultKtDepMarkerStyle[iCharge]);

      gMultKtDepRadiiMtGrErr[iCharge][iCent]->SetMarkerColor(mMultKtDepColor[iCent]);
      gMultKtDepRadiiMtGrErr[iCharge][iCent]->SetLineColor(mMultKtDepColor[iCent]);
      gMultKtDepRadiiMtGrErr[iCharge][iCent]->SetMarkerSize(mMultKtDepMarkerSize);
      gMultKtDepRadiiMtGrErr[iCharge][iCent]->SetMarkerStyle(mMultKtDepMarkerStyle[iCharge]);
    } //for(Int_t iCent=0; iCent<mMultBins; iCent++)
  } //for(Int_t iCharge=0; iCharge<3; iCharge++)
}

//_________________
vector<double> DoFit(TH1* histo, Int_t PartType, Int_t SourceFunc, 
		     Int_t NonFemtoFunc, Bool_t DrawFitLegend, Bool_t mExcludePoints) {

  vector<double> mRetVect;

  Int_t mNbins = histo->GetXaxis()->GetNbins();
  Double_t mFitMinVal = 0.;
  Double_t mDeltaMaxFitVal = -0.009; // -0.005
  Double_t mFitMaxVal = histo->GetXaxis()->GetBinUpEdge(mNbins) + mDeltaMaxFitVal;
  Int_t mFemtoParNum = 2;       //lambda, R
  Int_t mNonFemtoParNum = 0;

  if(NonFemtoFunc == 0)
    mNonFemtoParNum = 1;
  else if(NonFemtoFunc == 1)
    mNonFemtoParNum = 2;
  else if(NonFemtoFunc >= 2 && NonFemtoFunc<=3)
    mNonFemtoParNum = 3;
  else if(NonFemtoFunc >= 4)
    mNonFemtoParNum = 2;

  Int_t mFullParNum = mFemtoParNum + mNonFemtoParNum;
  HbtFitter* mHbtFitter = new HbtFitter(PartType,SourceFunc,NonFemtoFunc);
  mHbtFitter->Set1DMaxVal(histo->GetXaxis()->GetBinUpEdge(mNbins));
  mHbtFitter->Set1DBinsNum(mNbins);
  TF1* mCorrFit = new TF1("mCorrFit",mHbtFitter,&HbtFitter::FitFunction,
                          mFitMinVal,mFitMaxVal,mFullParNum,"HbtFitter","FitFunction");
  mCorrFit->SetParameter(0, 0.5);
  mCorrFit->SetParLimits(0, 0., 1.05);
  mCorrFit->SetParameter(1, 15.);
  mCorrFit->SetParLimits(1, 0., 50.);
  mCorrFit->SetParameter(2, 1.);
  //mCorrFit->SetParLimits(2, 0.9, 1.1);
  mCorrFit->SetLineWidth(3);
  mCorrFit->SetLineColor(kRed);
  mCorrFit->SetNpx(40);
  histo->Fit("mCorrFit","MRQE","ep");
  mCorrFit->Draw("same");

  if(SourceFunc==0) {
    mRetVect.push_back(mCorrFit->GetParameter(1)*MHC);                 //0 R
    mRetVect.push_back(mCorrFit->GetParError(1)*MHC);                  //1 dR
  }
  else if(SourceFunc==1) {
    mRetVect.push_back(mCorrFit->GetParameter(1)*MHC*INVSQRTPI);       //0 R
    mRetVect.push_back(mCorrFit->GetParError(1)*MHC*INVSQRTPI);        //1 dR
  }
  mRetVect.push_back(mCorrFit->GetParameter(0));                       //2 Lam
  mRetVect.push_back(mCorrFit->GetParError(0));                        //3 dLam
  mRetVect.push_back(mCorrFit->GetChisquare()/mCorrFit->GetNDF());     //4 chi2/ndf
  mRetVect.push_back(mCorrFit->GetParameter(2));                       //5 Norm
  mRetVect.push_back(mCorrFit->GetParError(2));                        //6 dNorm

  if(DrawFitLegend) {
    TLegend* mLegend = new TLegend(0.4,0.5,0.9,0.9);
    mLegend->SetFillColor(kWhite);
    mLegend->SetHeader("Fit parameters:");
    mLegend->AddEntry("mCorrFit","Bowler-Sinyukov fit","l");
    mLegend->AddEntry((TObject*)0, Form( "#chi^{2}/n.d.f. = %4.3f / %2d",
					 mCorrFit->GetChisquare(),mCorrFit->GetNDF() ),"");
    mLegend->AddEntry((TObject*)0, Form("N = %4.3f \\pm %4.3f",
					mRetVect.at(5),mRetVect.at(6) ), "");
    mLegend->AddEntry((TObject*)0, Form("#lambda = %4.2f \\pm %4.2f",
					mRetVect.at(2),mRetVect.at(3) ),"");
    mLegend->AddEntry((TObject*)0, Form("R = %4.2f \\pm %4.2f",
					mRetVect.at(0),mRetVect.at(1)), "");
    mLegend->Draw();
  }
  return mRetVect;
}