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

  char  hId[120];

  TFile *hFile = new TFile("./C_V1.root");

  TH2D *hMassXY; hFile->GetObject("hKpKmVsKpKmPi0",hMassXY);

  TF1 *f1 = new TF1("f1","gaus(0)+pol2(3)");
  TF1 *f1B = new TF1("f1B","gaus(0)+pol0(3)");
  TF1 *f0 = new TF1("f0","gaus");
  f1->SetLineWidth(2);
  f1->SetLineColor(4);

  //hMassXY->Draw("colz");
  //return;

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

  TH1D *hMassX = new TH1D("hMassX","",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();
  
  //hMassX->Draw();
  //return;

  TH1D *hMassY[nBinsX];
  TH1D *hMassYSignal[nBinsX];
  TH1D *hMassYBack[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,"hMassYBack_%d",iBinX);
    hMassYBack[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->SetParameter(1,1.02);
    f1->SetParameter(2,0.005);
    f1->FixParameter(1,par1);

    f1B->SetParameter(0,par0);
    //f1B->SetParameter(1,1.02);
    f1B->SetParameter(2,0.005);
    f1B->FixParameter(1,par1);

    //f1->SetParameter(2,par2);
    f1->SetParameter(3,0);
    f1->SetParameter(4,0);
    f1->SetParameter(5,0);
    //f1->SetParLimits(1,1.0,1.03);
    f1->SetParLimits(2,0.005,0.01);
    double fitHigh = 1.13;
    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 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);

    f1->SetParameter(0,0.0);
    hMassYBack[iBinX-1]->Eval(f1);
    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 backVal = hMassYBack[iBinX-1]->Integral();
    double varSignal = allVal + backVal;
    double sigmaSignal = 0;
    
    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);
    

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

    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");

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

  if (iShow < 0) {
    hMassX->GetXaxis()->SetTitle("mass(K^{+}K^{-}#pi^{0})");
    hMassX->GetXaxis()->SetRangeUser(1.22,1.35);
    hMassX->Draw("");
    hMassXB->SetLineColor(2);
    hMassXB->Draw("same");
  } else {
    //hMassY[iShow]->Draw("e");
    //hMassYBack[iShow]->Draw("csame");
    //hMassYSignal[iShow]->Draw("same");
  }
}
