class  MyFitFun1 {
public:
  double Evaluate(double *xPtr, double *parPtr) {
    double fitVal = 0.0;
    double xVal = *xPtr;

    //Voigtian fit parameters
    double gMag   = parPtr[0];
    double center = parPtr[1];
    double sigma  = parPtr[2];


    //Voigtian function
    double gFun = TMath::Gaus(xVal,center,sigma);
    //double gFun = 1.0;

    //Background fit parameters
    double p0Val = parPtr[3];
    double p1Val = parPtr[4];
    double p2Val = parPtr[5];

    //Background function
    double bkg = p0Val + p1Val*xVal + p2Val*xVal*xVal;

    //Fit value 

    fitVal =  gMag*gFun + bkg;

    return fitVal;
  }
};

class  MyFitFun1B {
public:
  double Evaluate(double *xPtr, double *parPtr) {
    double fitVal = 0.0;
    double xVal = *xPtr;

    //Voigtian fit parameters
    double gMag   = parPtr[0];
    double center = parPtr[1];
    double sigma  = parPtr[2];


    //Voigtian function
    double gFun = TMath::Gaus(xVal,center,sigma);
    //gFun = gMag*gFun;

    double sFun = 1.0/(exp((-xVal + 1.05)*100) + 1);

    //double gFun = 1.0;

    //Background fit parameters
    double p0Val = parPtr[3];
    double p1Val = parPtr[4];
    double p2Val = parPtr[5];

    //Background function
    //double bkg = sFun*gFun;

    //Fit value 
    gFun = 1;
    gMag = 200;
    fitVal =  sFun*gMag*gFun;

    return fitVal;
  }
};

class  MyFitFun2 {
public:
  double Evaluate(double *xPtr, double *parPtr) {
    double fitVal = 0.0;
    double xVal = *xPtr;

    double gMag   = parPtr[0];
    double center = parPtr[1];
    double sigma  = parPtr[2];

    double gFun = TMath::Gaus(xVal,center,sigma);


    //Background fit parameters
    double p0Val = parPtr[3];
    double p1Val = parPtr[4];
    double p2Val = parPtr[5];

    //Background function
    double bkg = p0Val + p1Val*xVal + p2Val*xVal*xVal;

    //a0(980)
    double gMag2   = parPtr[6];
    double center2 = parPtr[7];
    double sigma2  = parPtr[8];

    double gFun2 = TMath::Gaus(xVal,center2,sigma2);

    double sFun = 1.0/(exp((-xVal + 1.0)*200) + 1);

    //Fit value 
    //sFun = 0; 
    //gFun2 = 1;
    //gMag2 = 100;
    fitVal =  sFun*gMag*gFun + sFun*bkg + sFun*gMag2*gFun2;
    //fitVal = sFun*gMag2*gFun2;
    
    return fitVal;
  }
};


void fitPhiPeaks(int iShow){
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);

  char  hId[120];

  TFile *hFile = new TFile("./all_kpkmpi0.root");//CHANGE ME

  TH2D *hMassXY; hFile->GetObject("hKpKmVsKpKmPi0CutD4",hMassXY);//CHANGE ME

  MyFitFun1 * fptr1 = new MyFitFun1();
  int nPar = 6;
  TF1 *f1 = new TF1("f1",fptr1,&MyFitFun1::Evaluate,0.0,4.0,nPar,
                        "MyFitFun1","Evaluate");

  nPar = 9;
  MyFitFun2 * fptr2 = new MyFitFun2();
  TF1 *f2 = new TF1("f2",fptr2,&MyFitFun2::Evaluate,0.0,4.0,nPar,
                        "MyFitFun2","Evaluate");


  TF1 *f0 = new TF1("f0","gaus");
  f1->SetLineWidth(2);
  f1->SetLineColor(4);

  int    nBinsX = hMassXY->GetNbinsX();
  double xMin = hMassXY->GetXaxis()->GetXmin();
  double xMax = hMassXY->GetXaxis()->GetXmax();

  TH1D *hMassX = new TH1D("hMassX","",nBinsX,xMin,xMax);
  TH1D *hMassX1 = new TH1D("hMassX1","",nBinsX,xMin,xMax);
  TH1D *hMassX2 = new TH1D("hMassX2","",nBinsX,xMin,xMax);
  TH1D *hMassXB = new TH1D("hMassXB","",nBinsX,xMin,xMax);

  int    nBinsY = hMassXY->GetNbinsY();
  double yMin = hMassXY->GetYaxis()->GetXmin();
  double yMax = hMassXY->GetYaxis()->GetXmax();
  
  TH1D *hMassY[nBinsX];
  TH1D *hMassYSignal[nBinsX];
  TH1D *hMassYSignal1[nBinsX];
  TH1D *hMassYSignal2[nBinsX];
  TH1D *hMassYBack[nBinsX];
  TH1D *hMassYBack1[nBinsX];
  TH1D *hMassYBack2[nBinsX];

  int lBin = 1;
  int uBin = nBinsX;
  if (iShow > 0){
    lBin = iShow; 
    uBin = iShow; 
  }
  for (int iBinX = lBin; iBinX <= uBin; iBinX++) {
    sprintf(hId,"hMassY_%d",iBinX);
    hMassY[iBinX-1] = (TH1D*) hMassXY->ProjectionY(hId,iBinX,iBinX,"");
    sprintf(hId,"hMassYSignal_%d",iBinX);
    hMassYSignal[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax);

    sprintf(hId,"hMassYSignal1_%d",iBinX);
    hMassYSignal1[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax);

    sprintf(hId,"hMassYSignal2_%d",iBinX);
    hMassYSignal2[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax);

    sprintf(hId,"hMassYBack_%d",iBinX);
    hMassYBack[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax);

    sprintf(hId,"hMassYBack1_%d",iBinX);
    hMassYBack1[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax);

    sprintf(hId,"hMassYBack2_%d",iBinX);
    hMassYBack2[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax);

    hMassY[iBinX-1]->GetXaxis()->SetRangeUser(0.96,1.2);
    double maxVal = hMassY[iBinX-1]->GetMaximum();
    f0->SetParameter(0,maxVal);
    f0->SetParameter(1,1.02);
    f0->SetParameter(2,0.005);
    hMassY[iBinX-1]->Fit("f0","B","R",1.0,1.03);
    double par0 = f0->GetParameter(0);
    double par1 = f0->GetParameter(1);
    double par2 = f0->GetParameter(2);
    f1->SetParameter(0,par0);
    f1->FixParameter(2,7.12905e-03);
    f1->FixParameter(1,par1);
    f1->SetParameter(3,0);
    f1->SetParameter(4,0);
    f1->SetParameter(5,0);
    f1->SetParLimits(2,0.005,0.01);
    double fitHigh = 1.2;
    if (iBinX <= 66) fitHigh = 1.02;
    if (iBinX == 67) fitHigh = 1.03;
    if (iBinX == 68) fitHigh = 1.04;

    if (iBinX == 68 || iBinX == 69){
      f1->FixParameter(4,0);
      f1->FixParameter(5,0);
    }
    hMassY[iBinX-1]->Fit("f1","B","R",0.98,fitHigh);

    double f1Par0 = f1->GetParameter(0);
    double f1Par1 = f1->GetParameter(1);
    double f1Par2 = f1->GetParameter(2);
    double f1Par3 = f1->GetParameter(3);
    double f1Par4 = f1->GetParameter(4);
    double f1Par5 = f1->GetParameter(5);

    cout<<"f1Par1 = "<<f1Par1<<endl;

    f2->SetParameter(0,f1Par0);
    f2->SetParLimits(0,0.0,10000000);
    f2->FixParameter(1,f1Par1);
    f2->FixParameter(2,7.12905e-03);
    f2->FixParameter(3,f1Par3);
    f2->FixParameter(4,f1Par4);
    f2->FixParameter(5,f1Par5);

    f2->SetParameter(6,0);
    f2->SetParLimits(6,0.0,10000000);
    f2->FixParameter(7,0.980);
    f2->SetParameter(8,0.04);

    hMassY[iBinX-1]->Fit("f2","","R",0.98,fitHigh);


    double center = f1->GetParameter(1);
    double sigma = f1->GetParameter(2);
    double lowVal = center - 3*sigma; 
    double highVal = center + 3*sigma; 
    int binLow = hMassY[iBinX-1]->GetXaxis()->FindBin(lowVal);
    int binHigh = hMassY[iBinX-1]->GetXaxis()->FindBin(highVal);

    double f2par0 = f2->GetParameter(0);
    double f2par1 = f2->GetParameter(1);
    double f2par2 = f2->GetParameter(2);
    double f2par3 = f2->GetParameter(3);
    double f2par4 = f2->GetParameter(4);
    double f2par5 = f2->GetParameter(5);
    double f2par6 = f2->GetParameter(6);
    double f2par7 = f2->GetParameter(7);
    double f2par8 = f2->GetParameter(8);

    f2->SetParameter(0,0.0);
    f2->SetParameter(3,0.0);
    f2->SetParameter(4,0.0);
    f2->SetParameter(5,0.0);
    hMassYSignal1[iBinX-1]->Eval(f2);//a0(980)

    f2->SetParameter(0,f2par0);
    f2->SetParameter(6,0);
    hMassYSignal2[iBinX-1]->Eval(f2);//phi

    f2->SetParameter(0,0.0);
    f2->SetParameter(6,0.0);
    f2->SetParameter(3,f2par3);
    f2->SetParameter(4,f2par4);
    f2->SetParameter(5,f2par5);
    hMassYBack[iBinX-1]->Eval(f2);


    hMassYSignal[iBinX-1]->Add(hMassY[iBinX-1],hMassYBack[iBinX-1],1.0,-1.0);

    hMassYBack[iBinX-1]->SetLineWidth(2);
    hMassYBack[iBinX-1]->SetLineColor(2);

    double signalVal = hMassYSignal[iBinX-1]->Integral(binLow,binHigh);
    double allVal = hMassY[iBinX-1]->Integral(binLow,binHigh);
    double backVal = hMassYBack[iBinX-1]->Integral(binLow,binHigh);
    double varSignal = allVal + backVal;
    double sigmaSignal = sqrt(varSignal);

    double signalVal1 = hMassYSignal1[iBinX-1]->Integral(binLow,binHigh);

    double signalVal2 = hMassYSignal2[iBinX-1]->Integral(binLow,binHigh);
    
    if (varSignal > 0)  sigmaSignal = sqrt(varSignal);
    hMassYSignal[iBinX-1]->GetXaxis()->SetRangeUser(lowVal,highVal);
    hMassYSignal[iBinX-1]->SetFillStyle(3001);
    hMassYSignal[iBinX-1]->SetFillColor(3);
    hMassYSignal[iBinX-1]->SetLineColor(3);

    hMassYSignal1[iBinX-1]->SetLineColor(4);
    hMassYSignal2[iBinX-1]->SetLineColor(7);
    hMassYSignal1[iBinX-1]->SetLineWidth(2);
    hMassYSignal2[iBinX-1]->SetLineWidth(2);
    

    hMassX->SetBinContent(iBinX,signalVal);
    hMassX->SetBinError(iBinX,sigmaSignal);

    hMassX1->SetBinContent(iBinX,signalVal1);

    hMassX2->SetBinContent(iBinX,signalVal2);

    hMassXB->SetBinContent(iBinX,backVal);
    hMassXB->SetBinError(iBinX,sqrt(backVal));

    cout<<"signal = "<<signalVal<<" +/- "<<sigmaSignal<<endl;
    hMassY[iBinX-1]->Draw("e");
    hMassYBack[iBinX-1]->Draw("csame");
    hMassYSignal[iBinX-1]->Draw("same");



    hMassYSignal1[iBinX-1]->Draw("csame");
    hMassYSignal2[iBinX-1]->Draw("csame");

    double binCenter = hMassX->GetBinCenter(iBinX);
    cout<<"binCenter = "<<binCenter<<endl;
  }

  hMassX1->SetLineColor(4);
  hMassX2->SetLineColor(7);
  hMassX->GetXaxis()->SetRangeUser(1.2,2.4);

  if (iShow < 0) {
    hMassX->Draw("");
    hMassX1->Draw("same");
    hMassX2->Draw("same");
    hMassXB->SetLineColor(2);
    hMassXB->Draw("same");

  } else {
    //hMassY[iShow]->Draw("e");
    //hMassYBack[iShow]->Draw("csame");
    //hMassYSignal[iShow]->Draw("same");
  }
}