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);
  double bwMag = sqrt(bwReal*bwReal + bwImag*bwImag);
  return bwMag;
}

double bwSmear(double m,double m0,double Gamma,double m_smear,double m_centerShift, double *bwRealSmear, double *bwImagSmear){
  double gCenterArray[100] = {
    -2.97, -2.91, -2.85, -2.79, -2.73, -2.67, -2.61, -2.55, -2.49, -2.43,
    -2.37, -2.31, -2.25, -2.19, -2.13, -2.07, -2.01, -1.95, -1.89, -1.83,
    -1.77, -1.71, -1.65, -1.59, -1.53, -1.47, -1.41, -1.35, -1.29, -1.23,
    -1.17, -1.11, -1.05, -0.99, -0.93, -0.87, -0.81, -0.75, -0.69, -0.63,
    -0.57, -0.51, -0.45, -0.39, -0.33, -0.27, -0.21, -0.15, -0.09, -0.03,
    0.03, 0.09, 0.15, 0.21, 0.27, 0.33, 0.39, 0.45, 0.51, 0.57,
    0.63, 0.69, 0.75, 0.81, 0.87, 0.93, 0.99, 1.05, 1.11, 1.17,
    1.23, 1.29, 1.35, 1.41, 1.47, 1.53, 1.59, 1.65, 1.71, 1.77,
    1.83, 1.89, 1.95, 2.01, 2.07, 2.13, 2.19, 2.25, 2.31, 2.37,
    2.43, 2.49, 2.55, 2.61, 2.67, 2.73, 2.79, 2.85, 2.91, 2.97
  };
  double gValArray[100] = {
    0.000291608, 0.000347864, 0.000413481, 0.000489709, 0.000577906, 0.000679536, 0.000796168, 0.000929466, 0.00108118, 0.00125314,
    0.00144724, 0.00166538, 0.00190953, 0.00218159, 0.00248347, 0.00281695, 0.00318374, 0.00358535, 0.00402311, 0.0044981,
    0.0050111, 0.00556254, 0.00615248, 0.00678053, 0.00744584, 0.00814705, 0.00888226, 0.00964901, 0.0104443, 0.0112645,
    0.0121054, 0.0129624, 0.0138302, 0.0147031, 0.0155748, 0.016439, 0.0172887, 0.0181171, 0.0189169, 0.019681,
    0.0204025, 0.0210743, 0.0216901, 0.0222436, 0.0227293, 0.0231421, 0.0234778, 0.0237327, 0.0239042, 0.0239904,
    0.0239904, 0.0239042, 0.0237327, 0.0234778, 0.0231421, 0.0227293, 0.0222436, 0.0216901, 0.0210743, 0.0204025,
    0.019681, 0.0189169, 0.0181171, 0.0172887, 0.016439, 0.0155748, 0.0147031, 0.0138302, 0.0129624, 0.0121054,
    0.0112645, 0.0104443, 0.00964901, 0.00888226, 0.00814705, 0.00744584, 0.00678053, 0.00615248, 0.00556254, 0.0050111,
    0.0044981, 0.00402311, 0.00358535, 0.00318374, 0.00281695, 0.00248347, 0.00218159, 0.00190953, 0.00166538, 0.00144724,
    0.00125314, 0.00108118, 0.000929466, 0.000796168, 0.000679536, 0.000577906, 0.000489709, 0.000413481, 0.000347864, 0.000291608
  };


  double bwSmearVal = 0;
  double bwRealSmearVal = 0;
  double bwImagSmearVal = 0;
  if (m_smear > 0){ 
    for( unsigned int i = 0; i < 100; i++ ){
      double gcVal  = gCenterArray[i];
      double gVal   = sqrt(gValArray[i]);
      double deltaM = gcVal*m_smear + m_centerShift;
      
      bwSmearVal += gVal*bwMagFun(m+deltaM,m0,Gamma);
      bwRealSmearVal += gVal*bwRealFun(m+deltaM,m0,Gamma);
      bwImagSmearVal += gVal*bwImagFun(m+deltaM,m0,Gamma);
    }
  }else{
    bwSmearVal = bwMagFun(m,m0,Gamma);
    bwRealSmearVal = bwRealFun(m,m0,Gamma);
    bwImagSmearVal = bwImagFun(m,m0,Gamma);
  }
  *bwRealSmear = bwRealSmearVal;
  *bwImagSmear = bwImagSmearVal;
  return pow(bwSmearVal,2);
}


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

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

    double bwReal = bwRealFun(m,m0,Gamma);
    double bwMag = bwMagFun(m,m0,Gamma);

    double cosVal = bwReal/bwMag;

    //Fit value
    fitVal = TMath::ACos(cosVal);  

    return TMath::RadToDeg()*fitVal;
  }
};

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

    double m0    = parPtr[0];
    double Gamma = parPtr[1];
    double m_smear = parPtr[2];
    double m_centerShift = parPtr[3];

    double bwSmearVal = 0;
    double bwRealSmear = 0;
    double bwImagSmear = 0;
    return bwSmearVal = bwSmear(m,m0,Gamma,m_smear,m_centerShift,&bwRealSmear,&bwImagSmear); 
  }
};

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

    double m0    = parPtr[0];
    double Gamma = parPtr[1];
    double m_smear = parPtr[2];
    double m_centerShift = parPtr[3];

    double bwRealSmear = 0;
    double bwImagSmear = 0;

    double bwSmearVal = bwSmear(m,m0,Gamma,m_smear,m_centerShift,&bwRealSmear,&bwImagSmear); 

    double cosVal = bwRealSmear/sqrt(bwSmearVal);

    double phiVal = TMath::ACos(cosVal)*TMath::RadToDeg();

    return phiVal;
  }
};




void bwPhase(int iShow){

  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);


  double m0 = 892;
  double Gamma = 52;
  double mDelta = 4*Gamma;
  double m_smear = 9.0;
  double m_centerShift = 0.0;

  TH1D *hBWSmear = new TH1D("hBWSmear","",300,m0-mDelta,m0+mDelta);
  TH1D *hBW = new TH1D("hBW","",300,m0-mDelta,m0+mDelta);

  TH1D *hBWSmear1 = new TH1D("hBWSmear1","",300,m0-mDelta,m0+mDelta);
  TH1D *hBW1 = new TH1D("hBW1","",300,m0-mDelta,m0+mDelta);

  TH1D *hBWSmear2 = new TH1D("hBWSmear2","",300,m0-mDelta,m0+mDelta);
  TH1D *hBW2 = new TH1D("hBW2","",300,m0-mDelta,m0+mDelta);

  TH1D *hBWPhase = new TH1D("hBWPhase","",300,m0-mDelta,m0+mDelta);
  TH1D *hBWPhase1 = new TH1D("hBWPhase1","",300,m0-mDelta,m0+mDelta);
  TH1D *hBWPhase2 = new TH1D("hBWPhase2","",300,m0-mDelta,m0+mDelta);
  TH1D *hBWPhaseSmear = new TH1D("hBWPhaseSmear","",300,m0-mDelta,m0+mDelta);
  TH1D *hBWPhaseSmear1 = new TH1D("hBWPhaseSmear1","",300,m0-mDelta,m0+mDelta);
  TH1D *hBWPhaseSmear2 = new TH1D("hBWPhaseSmear2","",300,m0-mDelta,m0+mDelta);
  TH1D *hBWPhaseAve = new TH1D("hBWPhaseAve","",300,m0-mDelta,m0+mDelta);
  TH1D *hBWPhaseDiff = new TH1D("hBWPhaseDiff","",300,m0-mDelta,m0+mDelta);
  TH1D *hBWPhaseDiff12 = new TH1D("hBWPhaseDiff12","",300,m0-mDelta,m0+mDelta);
  TH1D *hBWPhaseDiffSmear12 = new TH1D("hBWPhaseDiffSmear12","",300,m0-mDelta,m0+mDelta);

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

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

  MyFitFun3 * fptr3 = new MyFitFun3();
  nPar = 4;
  TF1 *f3 = new TF1("f3",fptr3,&MyFitFun3::Evaluate,0.0,400.0,nPar,
		    "MyFitFun3","Evaluate");

  double m0Delta = 50;

  f2->SetParameter(0,m0);
  f2->SetParameter(1,Gamma);
  f2->SetParameter(2,m_smear);
  f2->SetParameter(3,m_centerShift);
  hBWSmear->Eval(f2);

  f2->SetParameter(0,m0+m0Delta);
  hBWSmear1->Eval(f2);

  f2->SetParameter(0,m0-m0Delta);
  hBWSmear2->Eval(f2);

  f2->SetParameter(0,m0);
  f2->SetParameter(2,-1);
  hBW->Eval(f2);

  f2->SetParameter(0,m0+m0Delta);
  hBW1->Eval(f2);

  f2->SetParameter(0,m0-m0Delta);
  hBW2->Eval(f2);

  f3->SetParameter(0,m0);
  f3->SetParameter(1,Gamma);
  f3->SetParameter(2,m_smear);
  f3->SetParameter(3,m_centerShift);
  hBWPhaseSmear->Eval(f3);

  f1->SetParameter(0,m0);
  f1->SetParameter(1,Gamma);
  hBWPhase->Eval(f1);



  //Phase 1
  f1->SetParameter(0,m0+m0Delta);
  hBWPhase1->Eval(f1);


  f3->SetParameter(0,m0+m0Delta);
  hBWPhaseSmear1->Eval(f3);

  //Phase 2
  f1->SetParameter(0,m0-m0Delta);
  hBWPhase2->Eval(f1);

  f3->SetParameter(0,m0-m0Delta);
  hBWPhaseSmear2->Eval(f3);

  hBWPhaseAve->Add(hBWPhase1,hBWPhase2);
  hBWPhaseAve->Scale(0.5);

  hBWPhaseDiff->Add(hBWPhase,hBWPhaseSmear,1.0,-1.0);
  hBWPhaseDiff12->Add(hBWPhase1,hBWPhase2,1.0,-1.0);
  hBWPhaseDiffSmear12->Add(hBWPhaseSmear1,hBWPhaseSmear2,1.0,-1.0);

  double intBWVal = hBW->Integral();
  double intBWSmearVal = hBWSmear->Integral();
  hBW->Scale(1.0/intBWVal);
  hBWSmear->Scale(1.0/intBWSmearVal);

  double intBW1Val = hBW1->Integral();
  double intBWSmear1Val = hBWSmear1->Integral();
  hBW1->Scale(1.0/intBW1Val);
  hBWSmear1->Scale(1.0/intBWSmear1Val);

  double intBW2Val = hBW2->Integral();
  double intBWSmear2Val = hBWSmear2->Integral();
  hBW2->Scale(1.0/intBW2Val);
  hBWSmear2->Scale(1.0/intBWSmear2Val);

  if(iShow == 91){ 
    hBWPhase1->Draw("");
    hBWPhase2->Draw("same");
    hBWPhaseDiff12->Draw("same");
  }

  hBWPhase1->GetXaxis()->SetTitle("Mass/MeV");
  hBWPhase1->GetYaxis()->SetTitle("Degree");
  hBWPhase1->SetLineColor(1);
  hBWPhase1->SetLineWidth(2);
  hBWPhase2->SetLineColor(1);
  hBWPhase2->SetLineWidth(2);

  hBWPhaseSmear1->SetLineColor(1);
  hBWPhaseSmear1->SetLineWidth(2);
  hBWPhaseSmear2->SetLineColor(1);
  hBWPhaseSmear2->SetLineWidth(2);

  hBWPhaseDiff12->SetLineColor(14);
  hBWPhaseDiff12->SetLineWidth(2);
  hBWPhaseDiffSmear12->SetLineColor(4);
  hBWPhaseDiffSmear12->SetLineWidth(2);


  if(iShow == 92){ 
    hBWPhase1->Draw("chist");
    hBWPhase2->Draw("chistsame");
    hBWPhaseDiff12->Draw("chistsame");
    hBWPhaseSmear1->SetLineColor(2);
    hBWPhaseSmear2->SetLineColor(2);
    hBWPhaseSmear1->Draw("chistsame");
    hBWPhaseSmear2->Draw("chistsame");
    hBWPhaseDiffSmear12->SetLineColor(4);
    hBWPhaseDiffSmear12->Draw("same");

    Double_t xl1=0.58, yl1 = 0.63, xl2=xl1+0.23, yl2=yl1+0.25;
    TLegend *legend = new TLegend(xl1,yl1,xl2,yl2);
    legend->SetFillColor(0);
    legend->SetLineColor(0);
    legend->SetTextSize(0.035);
    legend->AddEntry(hBWPhase1,"Phase: No smear","l");
    legend->AddEntry(hBWPhaseSmear1,"Phase: 9 MeV smear","l");
    legend->AddEntry(hBWPhaseDiff12,"#Delta Phase: No smear","l");
    legend->AddEntry(hBWPhaseDiffSmear12,"#Delta Phase: 9 MeV smear","l");
    legend->Draw();

    //xl1=0.12; yl1 = 0.33; xl2=xl1+0.23; yl2=yl1+0.25;
    //TLegend *legend2 = new TLegend(xl1,yl1,xl2,yl2);
    //legend2->SetFillColor(0);
    //legend2->SetLineColor(0);
    //legend2->SetTextSize(0.035);
    //legend2->AddEntry(hBWPhaseDiff12,"#Delta Phase: No smear","l");
    //legend2->AddEntry(hBWPhaseDiffSmear12,"#Delta Phase: 9 MeV smear","l");
    //legend2->Draw();

  }
  hBW1->GetXaxis()->SetTitle("Mass/MeV");
  hBW1->SetLineColor(1);
  hBW1->SetLineWidth(2);
  hBW2->SetLineColor(1);
  hBW2->SetLineWidth(2);
  hBWSmear1->SetLineColor(2);
  hBWSmear1->SetLineWidth(2);
  hBWSmear2->SetLineColor(2);
  hBWSmear2->SetLineWidth(2);
  if(iShow == 93){ 
    hBW1->Draw("chist");
    hBW2->Draw("chistsame");
    hBWSmear1->Draw("chistsame");
    hBWSmear2->Draw("chistsame");

    Double_t xl1=0.13, yl1 = 0.63, xl2=xl1+0.23, yl2=yl1+0.25;
    TLegend *legend = new TLegend(xl1,yl1,xl2,yl2);
    legend->SetFillColor(0);
    legend->SetLineColor(0);
    legend->SetTextSize(0.035);
    legend->AddEntry(hBW1,"No smear","l");
    legend->AddEntry(hBWSmear1,"9 MeV smear","l");
    legend->Draw();

  }

  hBW->SetLineColor(1);
  hBW->SetLineWidth(2);
  hBW->GetXaxis()->SetTitle("Mass/MeV");
  if(iShow == 1){ 
    hBW->Draw("chist");
  }

  hBWSmear->SetLineColor(2);
  hBWSmear->SetLineWidth(2);
  hBWSmear->GetXaxis()->SetTitle("Mass/MeV");
  if(iShow == 2){ 
    hBWSmear->Draw("chist");
  }

  if(iShow == 12){ 
    hBW->SetLineColor(1);
    hBWSmear->SetLineColor(2);
    hBW->Draw("chist");
    hBWSmear->Draw("chistsame");

    Double_t xl1=0.13, yl1 = 0.63, xl2=xl1+0.23, yl2=yl1+0.25;
    TLegend *legend = new TLegend(xl1,yl1,xl2,yl2);
    legend->SetFillColor(0);
    legend->SetLineColor(0);
    legend->SetTextSize(0.06);
    legend->AddEntry(hBW,"No smear","l");
    legend->AddEntry(hBWSmear,"9 MeV smear","l");
    legend->Draw();
  }

  hBWPhase->SetLineColor(1);   
  hBWPhase->SetLineWidth(2);   
  hBWPhase->GetXaxis()->SetTitle("Mass/MeV");
  hBWPhase->GetYaxis()->SetTitle("Phase/degree");

  hBWPhaseSmear->SetLineColor(2);   
  hBWPhaseSmear->SetLineWidth(2);
  if(iShow == 3){ 
    hBWPhase->Draw("chist");
    hBWPhaseSmear->Draw("chistsame");   
    Double_t xl1=0.53, yl1 = 0.63, xl2=xl1+0.23, yl2=yl1+0.25;
    TLegend *legend = new TLegend(xl1,yl1,xl2,yl2);
    legend->SetFillColor(0);
    legend->SetLineColor(0);
    legend->SetTextSize(0.06);
    legend->AddEntry(hBWPhase,"No smear","l");
    legend->AddEntry(hBWPhaseSmear,"9 MeV smear","l");
    legend->Draw();

  }
  if(iShow == 4){ 
    //hBWPhaseDiff->Scale(3.14159/180.0);
    hBWPhaseDiff->SetMaximum(90.0);
    hBWPhaseDiff->SetMinimum(-90.0);
    hBWPhaseDiff->Draw("chist");
  }
}