#include "TFile.h"
#include "TCanvas.h"
#include "TTree.h"
#include "TH1F.h"
#include "TF1.h"
#include "TMinuit.h"
#include "TMath.h"

//DoubleSidedCB function def.
Double_t func_DSCB(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] * func_DSCB(x[0], par[1],par[2],par[3],par[4],par[5],par[6])));
}


int double_SCB(){
//input file
TFile *file = new TFile( "/eos/user/k/ktariq/invar_mass_SS_data_v34.root" );
TH1F * h_mee = (TH1F*)file->Get("invar_mass");

//parameters values
TF1 *dcb = new TF1("dcb","DoubleSidedCB",50,130,7);
double par[7];
par[0]=1200;
par[1]=90;
par[2]=3;
par[3]=1.5;
par[4]=0.5;
par[5]=1.5;
par[6]=0.5;
dcb->SetParameters(&par[0]);
dcb->SetNpx(250);

TCanvas *C = new TCanvas("C","C",800,600);
TFitResultPtr r = h_mee->Fit(dcb,"SILRV+", "ep");
r->Print("V");     // print full information of fit including covariance matrix
gStyle->SetOptFit(111);
//gStyle->SetOptStat(0);

//Z_signal events (in mean+/-4Sigma) range 
double mean_fit = dcb->GetParameter("p1");
double sigma_fit = dcb->GetParameter("p2");
double z_sig_max = mean_fit + 4*sigma_fit;
double z_sig_min = mean_fit - 4*sigma_fit;

//Integral and IntegralError
double Z_sig = dcb->Integral(z_sig_min,z_sig_max);
double Z_sig_err =  dcb->IntegralError(z_sig_min,z_sig_max);

cout<<"Range Z_Sig: "<<z_sig_min<<" to "<<z_sig_max<<endl;
cout<<"Z_sig: "<<Z_sig<<" +/- "<<Z_sig_err<<endl;
return 0;
}