Double_t fTest(Double_t *x,Double_t *parPtr) {
  Double_t fitVal,xVal;

  xVal = *x;
  fitVal = 0.0;
  
  double gMag   = parPtr[0];
  double center = parPtr[1];
  double sigma  = parPtr[2];

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

  double bkg = parPtr[3] + parPtr[4]*xVal + parPtr[5]*pow(xVal,2);

  double sCenter1 = parPtr[6];
  double sWidth1 = parPtr[7];
  double sFun1 = 1.0/(exp((-xVal + sCenter1)*sWidth1) + 1);


  double sCenter2 = parPtr[8];
  double sWidth2 = parPtr[9];
  double sFun2 = 1.0/(exp((xVal - sCenter2)*sWidth2) + 1);

  double gMag2   = parPtr[10];
  double center2 = parPtr[11];
  double width2  = parPtr[12];


  //double gFun2 = TMath::Gaus(xVal,center2,sigma2)*gMag2;
  double gFun2 = TMath::Voigt(xVal-center2,sigma,width2)*gMag2;

  //sFun1 = 1;
  //sFun2 = 1;

  double gMag3   = parPtr[13];
  double center3 = parPtr[14];
  double width3  = parPtr[15];
  double gFun3 = TMath::Gaus(xVal,center3,width3)*gMag3;

  //double gMag4   = parPtr[16];
  //double center4 = parPtr[17];
  //double width4  = parPtr[18];
  //double gFun3 = TMath::Gaus(xVal,center3,width3)*gMag3;


  fitVal = fabs(gFun) + sFun1*bkg*sFun2 + fabs(gFun2) + fabs(gFun3);

  return fitVal;

}

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

  // .L fitCascadePeaks.C
  // fitCascadePeaks(-1)
  // .L fitCasStarAll.C
  // fitCasStarAll(1,6)

  char  hId[120];

  //TFile *hFile = new TFile("./all_pi0kpkpxim__B4_F1_M23.root");
  //TFile *hFile = new TFile("./all_pi0kpkpxim__B4_M23.root");
  TFile *hFileMC = new TFile("./mcTest.root");
  //TFile *hFile = new TFile("./allSp18.root");
  TFile *hFile = new TFile("./all18.root");
 

  //TH2D *hMassXY;     hFile->GetObject("hCasStarVsCas_MC_Kin_3",hMassXY);
  //TH2D *hMassXY;     hFile->GetObject("hCasStarVsCas_MC_Kin_4KCut",hMassXY);
  //TH2D *hMassXY;     hFile->GetObject("hCasStarVsCas_MC_Kin_4T",hMassXY);
  TH2D *hMassXY;     hFile->GetObject("hCasStarVsCas_MC_Kin_4B",hMassXY);
  //TH2D *hMassXY;     hFile->GetObject("hCasStarVsCas_MC_Kin_4",hMassXY);
  TH2D *hMassXYMC; hFileMC->GetObject("hCasStarVsCas_MC_Kin_4",hMassXYMC);

  //hMassXY->RebinX(2);//1
  //hMassXY->RebinY(2);//2
  TF1 *f2 = new TF1("f2",fTest,0,100,16);
  TF1 *f1 = new TF1("f1","gaus(0)+pol2(3)");
  TF1 *f0 = new TF1("f0","gaus");
  f2->SetLineWidth(2);
  f2->SetLineColor(1);

  //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 *hMassX2 = new TH1D("hMassX2","",nBinsX,xMin,xMax);

  int    nBinsY = hMassXY->GetNbinsY();
  double yMin = hMassXY->GetYaxis()->GetXmin();
  double yMax = hMassXY->GetYaxis()->GetXmax();
  
  //hMassX->Draw();
  //return;

  TH1D *hMassY[nBinsX];
  TH1D *hMassYMC[nBinsX];
  TH1D *hMassYSignal[nBinsX];
  TH1D *hMassYSignal2[nBinsX];
  TH1D *hMassYBack[nBinsX];

  int lBin = 1;
  int uBin = nBinsX;
  if (iShow > 0){
    lBin = iShow; 
    uBin = iShow; 
  }
  //hMassXY->Draw();
  //return;
  for (int iBinX = lBin; iBinX <= uBin; iBinX++) {
    sprintf(hId,"hMassY_%d",iBinX);
    hMassY[iBinX-1] = (TH1D*) hMassXY->ProjectionY(hId,iBinX,iBinX,"");
    sprintf(hId,"hMassYMC_%d",iBinX);
    hMassYMC[iBinX-1] = (TH1D*) hMassXYMC->ProjectionY(hId,iBinX,iBinX,"");
    sprintf(hId,"hMassYSignal_%d",iBinX);
    hMassYSignal[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);

    hMassY[iBinX-1]->GetXaxis()->SetRangeUser(1.2,1.8);
    double maxVal = hMassY[iBinX-1]->GetMaximum();
    double maxValMC = hMassYMC[iBinX-1]->GetMaximum();
    //f0->SetParameter(0,maxVal);
    //f0->FixParameter(1,1.323);
    //f0->FixParameter(2,4.03432e-03);
    //f0->SetParLimits(2,0.003,0.01);
    //hMassY[iBinX-1]->Fit("f0","LB","R",1.31,1.33);

    //f2->SetParameter(0,maxVal);
    //f2->FixParameter(1,1.323);
    //f2->FixParameter(2,0.006);
    //hMassY[iBinX-1]->Fit("f2","LB","R",1.31,1.33);
    


    //double par0 = f0->GetParameter(0);
    //double par1 = f0->GetParameter(1);
    //double par2 = f0->GetParameter(2);
    bool yesMC = true;
    if (maxValMC < 1) yesMC = false;
    f2->SetParameter(0,maxVal/2);
    f2->FixParameter(1,1.323);
    f2->FixParameter(2,4.03432e-03);
    f2->SetParameter(3,20);
    f2->SetParameter(4,0);
    f2->FixParameter(5,0);
    f2->SetParameter(6,1.3);
    //f2->SetParameter(6,1.26);
    f2->SetParameter(7,200);
    //f2->SetParameter(8,1.5);
    //f2->SetParameter(9,200);
    f2->SetParLimits(6,1.26,1.33);
    f2->SetParLimits(7,50,1000);
    f2->SetParLimits(8,1.45,1.6);
    f2->SetParLimits(9,100,1000);

    f2->SetParameter(10,10);
    f2->SetParameter(11,1.3837);
    f2->SetParLimits(11,1.3837-0.001,1.3837+0.001);
    f2->SetParameter(12,0.036);
    f2->SetParLimits(12,0.036-0.005,0.035+0.005);

    f2->SetParameter(13,10);
    f2->SetParameter(14,1.285);
    f2->SetParLimits(14,1.28,1.288);
    f2->SetParameter(15,0.01);
    f2->SetParLimits(15,0.01-0.005,0.01+0.005);

    //f1->SetParLimits(2,0.00,par2);
    //double fitHigh = 1.11;
    //if (iBinX <= 66) fitHigh = 1.02;
    //if (iBinX == 67) fitHigh = 1.03;
    //hMassY[iBinX-1]->Fit("f1","BL","R",1.28,1.45);
    if (yesMC == true){
      hMassYMC[iBinX-1]->Fit("f2","BL","R",1.26,1.36);
      double fracVal = f2->GetParameter(3);
      //f2->FixParameter(3,fracVal);
    } else {
      //f2->FixParameter(3,0.2);
      cout<<"mc false"<<endl;
    }
    
    //hMassYMC[iBinX-1]->Draw();
    //return;
    //hMassY[iBinX-1]->Fit("f2","BL","R",1.28,1.36);
    double fitLow = 1.3;
    double fitHigh = 1.44;

    hMassY[iBinX-1]->Fit("f2","BL","R",1.26,fitHigh);
    hMassY[iBinX-1]->Draw("e");

    //hMassY[iBinX-1]->Draw();
    //return;
    //return;
    double center = f2->GetParameter(1);
    double sigma = f2->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 center2 = f2->GetParameter(11);
    double sigma2 = f2->GetParameter(12);
    double lowVal2 = center2 - sigma2; 
    double highVal2 = center2 + sigma2; 
    int binLow2 = hMassY[iBinX-1]->GetXaxis()->FindBin(lowVal2);
    int binHigh2 = hMassY[iBinX-1]->GetXaxis()->FindBin(highVal2);

    f2->SetParameter(0,0.0);
    f2->SetParameter(10,0.0);
    hMassYBack[iBinX-1]->Sumw2();
    hMassYBack[iBinX-1]->Eval(f2);
    hMassYSignal[iBinX-1]->Add(hMassY[iBinX-1],hMassYBack[iBinX-1],1.0,-1.0);
    hMassYSignal2[iBinX-1]->Add(hMassY[iBinX-1],hMassYBack[iBinX-1],1.0,-1.0);
    hMassYBack[iBinX-1]->SetLineWidth(2);
    hMassYBack[iBinX-1]->SetLineColor(2);

    double mcVal = hMassYMC[iBinX-1]->Integral(binLow,binHigh);
    if (mcVal > 0){
      
      hMassYSignal[iBinX-1]->Scale(100.0/mcVal);
      hMassYSignal2[iBinX-1]->Scale(100.0/mcVal);
      //hMassYSignal[iBinX-1]->Scale(1.0);
    } else {
      //hMassYSignal[iBinX-1]->Scale(0.0);
    }

    double sigmaSignal = 0;
    double signalVal = hMassYSignal[iBinX-1]->IntegralAndError(binLow,binHigh,sigmaSignal,"");

    double sigmaSignal2 = 0;
    double signalVal2 = hMassYSignal[iBinX-1]->IntegralAndError(binLow2,binHigh2,sigmaSignal2,"");

    double allVal = hMassY[iBinX-1]->Integral(binLow,binHigh);
    double backVal = hMassYBack[iBinX-1]->Integral(binLow,binHigh);

    double allVal2 = hMassY[iBinX-1]->Integral(binLow2,binHigh2);
    double backVal2 = hMassYBack[iBinX-1]->Integral(binLow2,binHigh2);
    //double varSignal = allVal + backVal;
    //double sigmaSignal = 0;
    
    

    //if (varSignal > 0)  sigmaSignal = sqrt(varSignal);
    //if (varSignal < 0) signalVal = 0.0;
    hMassYSignal[iBinX-1]->GetXaxis()->SetRangeUser(lowVal,highVal);
    hMassYSignal[iBinX-1]->SetFillStyle(3001);
    hMassYSignal[iBinX-1]->SetFillColor(4);
    hMassYSignal[iBinX-1]->SetLineColor(4);

    hMassYSignal2[iBinX-1]->GetXaxis()->SetRangeUser(lowVal2,highVal2);
    hMassYSignal2[iBinX-1]->SetFillStyle(3001);
    hMassYSignal2[iBinX-1]->SetFillColor(3);
    hMassYSignal2[iBinX-1]->SetLineColor(3);
    
    double eTest;
    int yBinLow = hMassYSignal[iBinX-1]->GetXaxis()->FindBin(lowVal);
    int yBinHigh = hMassYSignal[iBinX-1]->GetXaxis()->FindBin(highVal);
    double yTest = hMassYSignal[iBinX-1]->IntegralAndError(yBinLow,yBinHigh,eTest,"");

    cout<<"yTest = "<<yTest<<" +/- "<<eTest<<endl;
    sigmaSignal = eTest;

    double yMin1 = hMassYSignal[iBinX-1]->GetMinimum();
    double yMin2 = hMassYSignal2[iBinX-1]->GetMinimum();
    double yMin = yMin1;
    if (yMin2 < yMin) yMin = yMin2; 
    //ASDF
    hMassX->SetBinContent(iBinX,signalVal);
    hMassX->SetBinError(iBinX,sigmaSignal);
    hMassX2->SetBinContent(iBinX,signalVal2);
    hMassX2->SetBinError(iBinX,sigmaSignal2);
    cout<<"signal = "<<signalVal<<" +/- "<<sigmaSignal<<endl;
    hMassY[iBinX-1]->SetLineColor(1);
    hMassY[iBinX-1]->SetMarkerColor(1);
    hMassY[iBinX-1]->SetMarkerStyle(21);
    hMassY[iBinX-1]->GetXaxis()->SetRangeUser(1.23,fitHigh);
    //hMassY[iBinX-1]->GetXaxis()->SetRangeUser(1.23,1.6);
    hMassY[iBinX-1]->SetMinimum(yMin);
    hMassY[iBinX-1]->Draw("e1");


    //hMassY2->GetXaxis()->SetRangeUser(lowVal2,highVal2);
    //hMassY2->Draw("e1same");
    hMassYBack[iBinX-1]->Draw("csamehist");
    hMassYSignal[iBinX-1]->Draw("samehist");
    hMassYSignal2[iBinX-1]->Draw("samehist");
    //hMassYMC[iBinX-1]->Draw("same");
    double binCenter = hMassX->GetBinCenter(iBinX);
    cout<<"binCenter = "<<binCenter<<endl;
    cout<<"mcVal = "<<mcVal<<endl;
    //hMassYMC[iBinX-1]->Draw("histsame");
  }

  

  if (iShow == -1) {
    hMassX->SetMinimum(0);
    hMassX->SetMaximum(400);
    hMassX->SetLineColor(4);
    hMassX->Draw();
  }
  if (iShow == -2) {
    hMassX2->SetMinimum(0);
    hMassX2->SetMaximum(200);
    hMassX2->SetLineColor(3);
    hMassX2->Draw("e1");
  }
  TFile *outFile = new TFile("./outFile.root","RECREATE");
  //outFile->cd();
  hMassX->Write();
  outFile->Close();
}