//please look after line 135

#ifndef __CINT__
#include <RooGlobalFunc.h>
#endif
#include <Riostream.h>
#include <RooAbsReal.h>
#include <RooAbsCategory.h>
#include <RooRealProxy.h>
#include <RooCategoryProxy.h>
#include <RooRealVar.h>
#include <RooDataSet.h>
#include <RooDataHist.h>
#include <RooGaussian.h>
#include <RooAbsPdf.h>
#include <RooPlot.h>

#include <TStyle.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TF1.h>
#include <TMath.h> 

double mygauss(double * x, double * par);
double mygauss_func(double _x, double const1, double mean, double width, double const2, double e_index);

using namespace RooFit ;
using namespace std;

class my_gaus_exp : public RooAbsPdf {
public:
  my_gaus_exp() {} ; 
  my_gaus_exp(const char *name, const char *title,
        RooAbsReal& _x,
        RooAbsReal& _mean,
        RooAbsReal& _width,
        RooAbsReal& _const1,
        RooAbsReal& _const2,
        RooAbsReal& _e_index );
  my_gaus_exp(const my_gaus_exp& other, const char* name=0) ;
  virtual TObject* clone(const char* newname) const { return new my_gaus_exp(*this,newname); }
  inline virtual ~my_gaus_exp() { }

protected:

  RooRealProxy x ;
  RooRealProxy mean ;
  RooRealProxy width ;
  RooRealProxy const1 ;
  RooRealProxy const2 ;
  RooRealProxy e_index ;
  
  Double_t evaluate() const ;

};

 my_gaus_exp::my_gaus_exp(const char *name, const char *title, 
                        RooAbsReal& _x,
                        RooAbsReal& _mean,
                        RooAbsReal& _width,
                        RooAbsReal& _const1,
                        RooAbsReal& _const2,
                        RooAbsReal& _e_index ) :
   RooAbsPdf(name,title), 
   x("x","x",this,_x),
   mean("mean","mean",this,_mean),
   width("width","width",this,_width),
   const1("const1","const1",this,_const1),
   const2("const2","const2",this,_const2),
   e_index("e_index","e_index",this,_e_index)
 { 
 } 

  my_gaus_exp::my_gaus_exp(const my_gaus_exp& other, const char* name) :  
   RooAbsPdf(other,name), 
   x("x",this,other.x),
   mean("mean",this,other.mean),
   width("width",this,other.width),
   const1("const1",this,other.const1),
   const2("const2",this,other.const2),
   e_index("e_index",this,other.e_index)
 { 
 } 
 
 Double_t my_gaus_exp::evaluate() const 
 { 
   return  mygauss_func(x, const1, mean, width, const2, e_index);
 } 
 
void ex3();
int main()
{
   ex3();
}
void ex3()
{
  //please look after line 135
  TCanvas* can = new TCanvas("can","can",600,600) ; can->Divide(2,3);
  
  //create base variables
  RooRealVar t("t","t",0,50) ;
  RooPlot* frame = t.frame(Title("gauss+exp")) ;  
  
  RooRealVar const1  ("const1", "const1", 10., 1., 100.);
  RooRealVar mean    ("mean", "mean", 25., 1., 20.);
  RooRealVar mean2   ("mean2", "mean2", 35., 1., 50.);
  RooRealVar width   ("width", "width", 3., 1., 10.);
  RooRealVar width2  ("width2", "width2", 5., 1., 10.);
  RooRealVar const2  ("const2", "const2", 5., 1., 10.);
  RooRealVar e_index ("e_index", "e_index", -0.1, -1., -0.0001);
  
  //histogram and function from ROOT
  TH1F* hist = new TH1F("h", "h", 100, 0, 50);
  TF1*  func = new TF1 ("func",mygauss, 0, 50, 5);
  
  //set function parameters
  func->SetParName(0,"const1");  func->SetParameter(0, const1 .getValV());
  func->SetParName(1,"mean");    func->SetParameter(1, mean   .getValV());
  func->SetParName(2,"width");   func->SetParameter(2, width  .getValV());
  func->SetParName(3,"const2");  func->SetParameter(3, const2 .getValV());
  func->SetParName(4,"e_index"); func->SetParameter(4, e_index.getValV());

  cout << "Before fit:"<< endl
       << "const1 = "  << const1 .getValV() << endl
       << "mean = "    << mean   .getValV() << endl
       << "width = "   << width  .getValV() << endl
       << "const2 = "  << const2 .getValV() << endl
       << "e_index = " << e_index.getValV() << endl;
  
  //fill and fit test histogram
  hist->FillRandom("func", 2000);  hist->Fit("func", "QR"); can->cd(1); can->cd(2); hist->Draw();
  can->Print("pic1.png");   
  
  //import TH1F into roofit and fit it with my function
  RooDataHist* rf_hist = new RooDataHist("rf_hist","rf_hist",t,hist);
  my_gaus_exp mge("mge","mge", t, mean, width, const1, const2, e_index) ;
  //fit here
  mge.fitTo(*rf_hist);
  rf_hist->plotOn(frame) ;
  mge.plotOn(frame);
  can->cd(2); frame->Draw() ;
  can->Print("pic2.png");
  
  cout << "After fit:" << endl
       << "const1 = "  << const1 .getValV() << endl
       << "mean = "    << mean   .getValV() << endl
       << "width = "   << width  .getValV() << endl
       << "const2 = "  << const2 .getValV() << endl
       << "e_index = " << e_index.getValV() << endl;
       
  //set func values to the ones from fit
  func->SetParName(0,"const1");  func->SetParameter(0, const1 .getValV());
  func->SetParName(1,"mean");    func->SetParameter(1, mean   .getValV());
  func->SetParName(2,"width");   func->SetParameter(2, width  .getValV());
  func->SetParName(3,"const2");  func->SetParameter(3, const2 .getValV());
  func->SetParName(4,"e_index"); func->SetParameter(4, e_index.getValV());
  TH1F* hhist = new TH1F(*hist);  can->cd(3); hhist->Draw();  func->Draw("same");
  can->Print("pic3.png");

}

double mygauss(double * x, double * par)
{
  double arg = 0;
  if (par[2] <  0) par[2]=-par[2];
  if (par[2] != 0) arg = (x[0] - par[1])/par[2];

  return par[0] * TMath::Exp(-0.5*arg*arg)  
  + par[3] *  TMath::Exp(x[0] * par[4]); 
}

double mygauss_func(double _x, double const1, double mean, double width, double const2, double e_index)
{
   double x[1];
   double par[5];
   
   x[0] = _x;
   par[0] = const1;
   par[1] = mean;
   par[2] = width;
   par[3] = const2;
   par[4] = e_index;
   
   return mygauss(x, par);
}