double bwRealFun(double m,double m0,double Gamma){
  double bwReal = (m*m - m0*m0)/(pow(m*m - m0*m0,2)+Gamma*Gamma*m0*m0);
  return bwReal;
}
double bwImagFun(double m,double m0,double Gamma){
  double bwImag = Gamma*m0/(pow(m*m - m0*m0,2)+Gamma*Gamma*m0*m0);
  return bwImag;
}
double bwMagFun(double m,double m0,double Gamma){
  double bwReal = bwRealFun(m,m0,Gamma);
  double bwImag = bwImagFun(m,m0,Gamma);
  //bwImag = 0;
  //bwReal = 0;
  double bwMag = sqrt(bwReal*bwReal + bwImag*bwImag);
  return bwMag;
}
class  MyFitFunBW {
public:
  double Evaluate(double *xPtr, double *parPtr) {
    double fitVal = 0.0;
    double m = *xPtr;

    double m0    = parPtr[0];
    double Gamma = parPtr[1];

    double bwMag = bwMagFun(m,m0,Gamma);
    fitVal = bwMag;
    return fitVal;
  }
};
double pCmVal(double m,double m1,double m2){
  double pcM = sqrt(fabs((m*m - pow(m1+m2,2))*(m*m - pow(m1-m2,2))))/(2*m);
  return pcM;
}
double flatteRealFun(double m,double mR,double g1,double g2){
  double mK = 0.493677;
  double mEta = 0.547862;
  double mPi0 = 0.1349768;
  double Gamma0 = g1*pCmVal(mR,mEta,mPi0);
  double GammaEtaPi0 = g1*pCmVal(m,mEta,mPi0);
  double GammaKK = g2*pCmVal(m,mK,mK);

  double flatteReal = 0;

  //flatteReal = mR*sqrt(Gamma0*GammaKK)*(mR*mR - m*m)/
  //(pow(mR*mR - m*m,2)+mR*mR*pow(GammaEtaPi0 + GammaKK,2));

  flatteReal = mR*(mR*mR - m*m)/
    (pow(mR*mR - m*m,2)+mR*mR*pow(GammaEtaPi0 + GammaKK,2));

  //flatteReal = mR*sqrt(Gamma0*GammaEtaPi0)*(mR*mR - m*m)/
  //(pow(mR*mR - m*m,2)+mR*mR*pow(GammaEtaPi0 + GammaKK,2));

  if (m < mK + mK) flatteReal = 0;  

  return flatteReal;
}

double flatteImagFun(double m,double mR,double g1,double g2){
  double mK = 0.493677;
  double mEta = 0.547862;
  double mPi0 = 0.1349768;
  double Gamma0 = g1*pCmVal(mR,mEta,mPi0);
  double GammaEtaPi0 = g1*pCmVal(m,mEta,mPi0);
  double GammaKK = g2*pCmVal(m,mK,mK);

  double flatteImag = 0;
  flatteImag = mR*sqrt(Gamma0*GammaKK)*mR*(GammaEtaPi0 + GammaKK)/
  (pow(mR*mR - m*m,2)+mR*mR*pow(GammaEtaPi0 + GammaKK,2));

  //flatteImag = mR*(GammaEtaPi0 + GammaKK)/
  //(pow(mR*mR - m*m,2)+mR*mR*pow(GammaEtaPi0 + GammaKK,2));

  //flatteImag = mR*sqrt(Gamma0*GammaEtaPi0)*mR*(GammaEtaPi0 + GammaKK)/
  //(pow(mR*mR - m*m,2)+mR*mR*pow(GammaEtaPi0 + GammaKK,2));
  
  if (m < mK + mK) flatteImag = 0;

  return flatteImag;
}
double flatteMagFun(double m,double m0,double g1,double g2){
  double flatteReal = flatteRealFun(m,m0,g1,g2);
  double flatteImag = flatteImagFun(m,m0,g1,g2);

  //flatteReal = 0;
  //flatteImag = 0;

  double flatteMag = sqrt(flatteReal*flatteReal + flatteImag*flatteImag);

  return flatteMag;
}



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

    double amp = parPtr[0];
    double m0 = parPtr[1];
    double g1 = parPtr[2];
    double g2 = parPtr[3];

    double flatteMag = amp*flatteMagFun(m,m0,g1,g2);
    
    fitVal = flatteMag;
    return fitVal;
  }
};

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

  double m0 = 0.980;
  double Gamma = 0.052;

  MyFitFunBW * fptr1 = new MyFitFunBW();
  int nPar = 2;
  TF1 *f1 = new TF1("f1",fptr1,&MyFitFunBW::Evaluate,0.0,4000.0,nPar,
                  "MyFitFunBW","Evaluate");

  MyFitFunFlatte * fptr2 = new MyFitFunFlatte();
  int nPar2 = 4;
  TF1 *f2 = new TF1("f2",fptr2,&MyFitFunFlatte::Evaluate,0.0,4000.0,nPar2,
                  "MyFitFunFlatte","Evaluate");

  TH1D *hBW = new TH1D("hBW","",300,m0-3*Gamma,m0+10*Gamma);



  f1->FixParameter(0,m0);
  f1->FixParameter(1,Gamma);
  //hBW->Eval(f1);

  //hBW->Scale(pow(10,8));

  hBW->FillRandom("f1",100000);

  


  f2->SetParameter(0,100);
  //f2->SetParameter(1,m0);
  f2->SetParameter(1,0.99);
  f2->FixParameter(2,0.212);
  f2->FixParameter(3,0.106);
  //f2->SetParameter(3,0.096721);

  hBW->Fit("f2","B","R",0.05,1.5);

  //hBW->Fit("gaus","B","R",0,100000);

  if (iShow == 1) {
    hBW->Draw();
  } 
}