#include "TLorentzVector.h"
#include <vector>
//#include "MyEvent.h"

#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TTree.h"
#include "TH1D.h"
#include "TRandom.h"
#include "RooGenericPdf.h"
#include "RooTFnBinding.h"
using namespace RooFit;



Double_t DoubleSidedCB2(double x, double mu, double width, double a1, double p1, double a2, double p2)
{
  double u   = (x-mu)/width;
  double A1  = TMath::Power(p1/TMath::Abs(a1),p1)*TMath::Exp(-a1*a1/2);
  double A2  = TMath::Power(p2/TMath::Abs(a2),p2)*TMath::Exp(-a2*a2/2);
  double B1  = p1/TMath::Abs(a1) - TMath::Abs(a1);
  double B2  = p2/TMath::Abs(a2) - TMath::Abs(a2);

  double result(1);
  if      (u<-a1) result *= A1*TMath::Power(B1-u,-p1);
  else if (u<a2)  result *= TMath::Exp(-u*u/2);
  else            result *= A2*TMath::Power(B2+u,-p2);
  return result;
}


double DoubleSidedCB(double* x, double *par)
{
  return(par[0] * DoubleSidedCB2(x[0], par[1],par[2],par[3],par[4],par[5],par[6]));
}

int DCB(){

TH1F *hmass = new TH1F("hmass","hmass",50,4.7,5.7);

TFile *f = new TFile("BDTresult_TagKMuMu_Mu9IP6_BDTtuple_v2.0_cmgTuple_v5.4/cuts_tagPt9_eta1p5_sdxy6_mu2Pt2_KPt1/MC_all_resonant.root","READ");
TTree *t1 = (TTree*)f->Get("mytreefit");

TH1F *real_mass = new TH1F("real_mass","real_mass",100,4.7,5.7);

Float_t Bmass,Mll;
Double_t xgb;
t1->SetBranchAddress("Bmass",&Bmass);
t1->SetBranchAddress("Mll",&Mll);
t1->SetBranchAddress("xgb",&xgb);

Int_t nentries = (Int_t)t1->GetEntries();

for (Int_t i=0; i<nentries; i++){
    t1->GetEntry(i);
    if(Mll>2.9 && Mll<3.2 && xgb>-1){real_mass->Fill(Bmass);}
}


TF1 *dcb = new TF1("dcb","DoubleSidedCB",4.7,5.7,7);
double par[7];
par[0]=35000;
par[1]=5.27;
par[2]=0.08;
par[3]=1;
par[4]=1;
par[5]=1;
par[6]=1;
dcb->SetParameters(&par[0]);
TCanvas *C = new TCanvas("C","C",700,700);
real_mass->Draw("PE1");
gStyle->SetOptStat(0);
real_mass->SetTitle(";m(KJ/#psi) [GeV];");
real_mass->Fit("dcb");
TLegend *LEG = new TLegend(0.6,0.8,0.89,0.89);
LEG->AddEntry(dcb,"Double Sided Crystal Ball");
LEG->SetBorderSize(0);
LEG->Draw("sames");


//Roofit
RooRealVar bMass("m(KJ/#psi)","m(KJ/#psi) [GeV]",4.7,5.7);
RooDataHist hist("hist","hist",bMass,Import(*real_mass));

//DCB parameters
RooRealVar mu("mu","mu",par[1],4.7,5.7);
RooRealVar width("width","width",par[2],0.,1.);
RooRealVar a1("a1","a1",par[3],0.,100.);
RooRealVar p1("p1","p1",par[4],0.,100.);
RooRealVar a2("a2","a2",par[5],0.,100.);
RooRealVar p2("p2","p2",par[6],0.,100.);

RooCrystalBall dcbPdf("dcbPdf","DoubleSidedCB",bMass,mu,width,a1,p1,a2,p2);


auto result = dcbPdf.fitTo(hist, RooFit::Save(true));
result->Print();

auto pl = bMass.frame();

hist.plotOn(pl);
dcbPdf.plotOn(pl);
pl->Draw(); 

return 0;
}