class  MyFitFun2 {
public:
  double Evaluate(double *xPtr, double *parPtr) {
    double fitVal = 0.0;
    double xVal = *xPtr;
    
    
    //KStar
    double mag       = fabs(parPtr[0]);
    double center    = parPtr[1];
    double width     = parPtr[2];
    double gausSmear = parPtr[3];
    double voigtFun  = TMath::Voigt(xVal-center,gausSmear,width);

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

    double bkg = fabs(p0Val) + p1Val*xVal + p2Val*xVal*xVal;

    double sCenter = parPtr[7];
    double sWidth = parPtr[8];
    double sFun = 1.0/(exp((-xVal + sCenter)*sWidth) + 1);

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


    //fitVal =  sFun*gMag*gFun + sFun*bkg + sFun*bwFun;                         
    //fitVal =  sFun*gMag*gFun + sFun*bkg + sFun*voigtFun;

    //sFun=1;
    fitVal =  sFun2*mag*voigtFun*sFun + sFun2*bkg*sFun;


    return fitVal;
  }
};


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

  char  hId[120];
  char  hId2[120];

  char  fId[120];

  char  tString1[120];


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



  sprintf(hId2,"hKpStarVsKpKmPi0W");
  if (kmStarFit == 1) sprintf(hId2,"hKmStarVsKpKmPi0W");



  TH2D *hMassXY; hFile->GetObject(hId2,hMassXY);



  TH1D *hMX = (TH1D*) hMassXY->ProjectionX(hId);

  TF1 *f1 = new TF1("f1","gaus(0)+pol2(3)");
  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 *hCenter = new TH1D("hCenter","",nBinsX,xMin,xMax);
  TH1D *hSigma = new TH1D("hSigma","",nBinsX,xMin,xMax);

  hCenter->Sumw2();
  hSigma->Sumw2();

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

  int nPar = 11;
  MyFitFun2 * fptr2 = new MyFitFun2();
  TF1 *f2 = new TF1("f2",fptr2,&MyFitFun2::Evaluate,0.0,4.0,nPar,
		    "MyFitFun2","Evaluate");
  
  TH1D *hMassY[nBinsX];
  TH1D *hMassYSignal[nBinsX];
  TH1D *hMassYBack[nBinsX];

  int lBin = 14;
  int uBin = nBinsX;
  if (iShow > 0){
    lBin = iShow; 
    uBin = iShow; 
  }
  cout<<"uBin = "<<uBin<<endl;
  for (int iBinX = lBin; iBinX <= uBin; iBinX++) {
    sprintf(hId,"hMassY_%d",iBinX);
    hMassY[iBinX-1] = (TH1D*) hMassXY->ProjectionY(hId,iBinX,iBinX,"");
    //hMassY[iBinX-1] = (TH1D*) hMassXY->ProjectionY(hId);

    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.8,1.2);
    double maxVal = hMassY[iBinX-1]->GetMaximum();


    f2->SetParameter(0,1);
    f2->SetParLimits(0,1,1000);
    //f2->SetParameter(1,0.89547);
    f2->SetParameter(1,0.89547);
    f2->SetParLimits(1,0.89547-0.001,0.89547+0.001);
    f2->SetParameter(2,0.0514);
    f2->SetParLimits(2,0.0514-0.00026,0.0514+0.00026);
    f2->SetParameter(3,0.001);
    f2->SetParLimits(3,0.0,0.007);

    f2->SetParameter(4,50.0);
    f2->SetParameter(5,0.0);
    f2->SetParameter(6,0.0);

    f2->SetParameter(7,0.91);
    f2->SetParLimits(7,0.88,1.0);
    f2->SetParameter(8,-2.93068e+02);
    f2->SetParLimits(8,-2000,-50);

    f2->SetParameter(9,0.75);
    f2->SetParLimits(9,0.6,0.85);

    f2->FixParameter(10,293);
    f2->SetParLimits(10,200,400);
    ////////////
    //f2->FixParameter(0,100);
    //f2->FixParameter(4,0);
    f2->SetParLimits(5,-500,9999999);



    f2->FixParameter(6,0);


    f2->SetParameter(9,0.6);

    double fitLow = 0.8;
    double fitHigh = 0.99;

    if (iBinX <= 27){
      cout<<"<= BIN15 !!"<<endl;
      
      f2->SetParameter(1,0.84);
      f2->SetParLimits(1,0.84,0.895);
      //f2->SetParameter(2,0.0514);
      f2->SetParLimits(2,0.02,0.0514+0.00026);
      //f2->SetParameter(3,0.001);
      //f2->SetParLimits(3,0.0,0.007);

      f2->SetParLimits(5,-200,9999999);
      fitLow = 0.75;
    }





    cout<<"t:4"<<endl;
    hMassY[iBinX-1]->Fit("f2","BQ","R",fitLow,fitHigh);
    cout<<"t:5"<<endl;

    double par2 = f2->GetParameter(2);
    double par3 = f2->GetParameter(3);

    double par7 = f2->GetParameter(7);
    double par8 = f2->GetParameter(8);



    f2->FixParameter(2,par2);
    f2->FixParameter(3,par3);

    f2->FixParameter(7,par7);
    f2->FixParameter(8,par8);

    f2->ReleaseParameter(10);




    fitLow = 0.6;
    //if (iBinX <= 15) fitLow = 0.75;

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

    //hMassY[iBinX-1]->Draw();
    //return;

    //hMassY[iBinX-1]->Draw();
    //continue;

    double center = f2->GetParameter(1);
    double centerErr = f2->GetParError(1);
    double sigma1 = f2->GetParameter(2);
    double sigmaErr = f2->GetParError(2);
    double sigma2 = f2->GetParameter(3);
    double sigma = sqrt(pow(sigma1/2.35,2)+pow(sigma2,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);

    f2->SetParameter(0,0.0);
    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 = 0;

    cout<<"xLow, xHigh = "<<lowVal<<" < "<<highVal<<endl;
    
    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);
    
    hCenter->SetBinContent(iBinX,center);
    hCenter->SetBinError(iBinX,centerErr);
    hSigma->SetBinContent(iBinX,sigma);
    hSigma->SetBinError(iBinX,sigmaErr);

    hMassX->SetBinContent(iBinX,signalVal);
    hMassX->SetBinError(iBinX,sigmaSignal);
    cout<<"signal = "<<signalVal<<" +/- "<<sigmaSignal<<endl;


    hMassY[iBinX-1]->GetXaxis()->SetTitle("Mass[K^{+}#pi^{0}]");
    if (kmStarFit == 1) hMassY[iBinX-1]->GetXaxis()->SetTitle("Mass[K^{-}#pi^{0}]");
    hMassY[iBinX-1]->GetYaxis()->SetTitle("Q-weighted counts");
    
    hMassY[iBinX-1]->GetXaxis()->SetTitleOffset(1.2);
    hMassY[iBinX-1]->GetXaxis()->SetTitleFont(62);
    hMassY[iBinX-1]->GetYaxis()->SetTitleFont(62);
    hMassY[iBinX-1]->GetXaxis()->SetLabelFont(62);
    hMassY[iBinX-1]->GetYaxis()->SetLabelFont(62);
    hMassY[iBinX-1]->GetXaxis()->SetLabelSize(0.04);
    hMassY[iBinX-1]->GetYaxis()->SetLabelSize(0.04);


    hMassY[iBinX-1]->GetXaxis()->SetRangeUser(0.6,1.0);

    hMassY[iBinX-1]->SetLineColor(1);
    hMassY[iBinX-1]->SetMarkerStyle(21);

    hMassY[iBinX-1]->SetMinimum(0.0);
    hMassY[iBinX-1]->Draw("e1");
    hMassYBack[iBinX-1]->Draw("csame");
    hMassYSignal[iBinX-1]->Draw("samee");
    hMassYSignal[iBinX-1]->Draw("samehist");
    
    TLatex *t1 = new TLatex();
    t1->SetNDC();
    double xLocationFraction = 0.1;
    double yLocationFraction = 0.92;

    //ASDF
    int iEnergy = (iBinX-14)*10 + 1355;

    sprintf(tString1,"M[K^{+}K^{-}#pi^{0}] = %d MeV",iEnergy);
    t1->DrawLatex(xLocationFraction,yLocationFraction,tString1);




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

  }

  if (iShow < 0) {
 
    hMX->Draw();
    hMassX->GetXaxis()->SetRangeUser(1.3,1.5);
    hMassX->Draw("same");
    //hMassX->Draw("");
  } else {
    //hMassY[iShow]->Draw("e");
    //hMassYBack[iShow]->Draw("csame");
    //hMassYSignal[iShow]->Draw("same");
  }

  hMX->SetName("hMX");

  if (iShow < 0) {
    sprintf(fId,"./myOutKp.root");
    if (kmStarFit == 1)  sprintf(fId,"./myOutKm.root");
    TFile *outFile = new TFile(fId,"RECREATE");
    hMX->Write();
    hMassX->Write();
    outFile->Close();
  }
  

}