#ifndef CINT
#include <RooGlobalFunc.h>
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TAxis.h"
#include <RooLandau.h>
#include <RooCBShape.h>
#include <RooFFTConvPdf.h>
#include "TROOT.h"
using namespace RooFit;


void RooFitConv()

{
  // std::string FileName(“124-175-3000ms+200ms-10Hf-ProtonGated-900+2000keV”);
  // std::string DirectoryName(“Correlator/DSSD_Decay+Clover”);
  // std::string HistName(“DSSD_Decay_!_Zn1”);
  
  TString inName1 = "/Users/mansisaxena/57ZnBetaP/e18018_proof/Root-Files-B/124-175-3000ms+200ms-10Hf-ProtonGated-900+2000keV.root";
    
  TFile* file = new TFile(inName1.Data(),"READ");
  if (file->IsOpen()) cout <<  "File " << inName1.Data() << " opened successfully" << endl;
  else {
    cout << "File " << inName1.Data() << " not opened successfully! Exiting..." <<endl;
    exit(EXIT_FAILURE);
  }
  double fitRange = 3E8;
  
  TString dirName = "Correlator/DSSD_Decay+Clover";
  file->cd(dirName.Data());
  TString histName = Form("DSSD_Decay_!_Zn1");
  TH1D * decayEnergy = (TH1D*) gDirectory->Get(histName.Data());
  decayEnergy->Rebin(10);
  decayEnergy->Sumw2(kTRUE);
  

  RooRealVar t("t", "t", 1800, 2500); //1800 2500
  // t.setBins(10000, "cache");
  RooDataHist* data = new RooDataHist("data", "data", t, Import(*decayEnergy));
  //landau function
  RooRealVar ml("ml", "mean landau", 2200, 1800, 2500); // 2100 1800 2400
  RooRealVar sl("sl", "sigma landau", 60, 1, 150);
  RooLandau landau("lx", "lx", t, ml, sl);
  //crystal ball function
  RooRealVar mc("mc","mc", 2200, 1800, 2500); // 2100 1800 2400
  RooRealVar sc("sc","sc", 50, 1, 150);
  RooRealVar alpha("alpha", "alpha", -1.20);//, 0, 100); // 1 0 100
  RooRealVar nn("nn", "nn", 5.1);//, 0, 100);
  RooCBShape crystal("crystal", "crystal", t, mc, sc, alpha, nn);
  //gaussian function
  RooRealVar mg("mg","mg",decayEnergy->GetMean(), 1800, 2500); //1800 2400
  RooRealVar sg("sg", "sg", decayEnergy->GetRMS(), 10, 150);
  RooGaussian gauss("gauss", "gauss", t, mg, sg);

  
  
  // std::stringstream fname;
  // fname << "/Users/mansisaxena/57ZnBetaP/e18018_proof/Root-Files-B/" << FileName << ".root";
  // TFile* f = TFile::Open(fname.str().c_str());
  // f->cd(DirectoryName.c_str());
  // TH1F* h =(TH1F*)f->FindObjectAny(HistName.c_str()); 


  //  t.setRange("signalRange",4200,5200) ;
  //convoluting and adding functions
  RooFFTConvPdf lxc("lxc", "Landau (X) CBall", t, landau, crystal);
  // Construct landau (x) gauss
  RooFFTConvPdf lxg("lxg","landau (X) gauss",t,landau,gauss);
  RooRealVar coeff("coeff", "coeff", .5, 0.01, 1);
  //landau (x) CB +  gauss
  RooAddPdf lxcg("lxcg", "Landau Convoluted Crystal Ball + Gaussian", lxc, gauss, coeff);

  //plotting and fitting functions to the data
  // RooPlot* frame = t.frame(Title("Landau Convoluted Crystal Ball + Gauss"));
  RooPlot* frame = t.frame(Title("Crystal Ball "));
  data->plotOn(frame);
  lxcg.fitTo(*data);
  crystal.fitTo(*data);
  // lxcg.plotOn(frame, LineColor(kBlue));
  // Overlay the background component of model with a dashed line
  //  lxcg.plotOn(frame, Components(lxc), LineStyle(kDashed),LineColor(kRed));
  //lxcg.plotOn(frame, Components(gauss), LineStyle(kDashed),LineColor(kRed));
  crystal.plotOn(frame, LineStyle(kDashed),LineColor(kGreen));
  // lxg.plotOn(frame, LineStyle(kDashed),LineColor(kRed));

  // Fit Landau x Gauss to data
  // landau.fitTo(*data,Extended()) ;
  gauss.fitTo(*data, Extended());
  // lxg.fitTo(*data, Extended());
    
  // Plot data, landau pdf, landau (X) gauss pdf
  RooPlot* frameL = t.frame(Title("landau (x) gauss convolution")) ;
  data->plotOn(frameL) ;
  gauss.plotOn(frameL,LineStyle(kDashed), LineColor(kRed)) ;
  // landau.plotOn(frameL,LineStyle(kDashed), LineColor(kBlue)) ;
  //landau.plotOn(frameL, LineColor(kRed)) ;
  landau.plotOn(frameL,LineStyle(kDashed), LineColor(kBlue)) ;
  // lxg.plotOn(frameL,Components(landau),LineStyle(kDashed),LineColor(kRed));
  // lxg.plotOn(frameL, Components(gauss), LineStyle(kDashed),LineColor(kRed));
  lxg.plotOn(frameL, LineStyle(kDashed), LineColor(kGreen));
  
    
   
 //integrating functions and calculating MPV
t.setRange("integral", 4400, 5200);
 RooAbsReal* Int = lxcg.createIntegral(t, Range("integral")); // landau x CB + Gauss
 RooAbsReal* lxcInt = lxc.createIntegral(t, Range("integral")); // landau x CB
 RooAbsReal* gaussInt = gauss.createIntegral(t, Range("integral")); //gauss
 double LxCInt = lxcInt->getVal();
 double GaussInt = gaussInt->getVal();
 double FullInt = Int->getVal();
 cout << "Integral of LxC = " << LxCInt << endl;
 cout << "Integral of gauss = " << GaussInt << endl;
 cout << "Full Integral = " << FullInt << endl;
 
 RooAbsReal* dldt = lxc.derivative(t) ;
Double_t tpeak ;
 Bool_t ok = RooBrentRootFinder(RooRealBinding(*dldt, t)).findRoot(tpeak, 1800, 2400, 0) ;
 cout << "MPV = " << tpeak << endl;
 
 //editing canvas and drawing it
 TCanvas* c1 = new TCanvas("gausscrystal", "gausscrystal", 1000, 1000);
 c1->Divide(2);
 c1->cd(1);
 gPad->SetLeftMargin(0.15);
 frame->GetYaxis()->SetTitleOffset(1.2);
 frame->GetYaxis()->SetTitle(" Counts");
 frame->GetXaxis()->SetTitle(" Decay Energy");
 frame->Draw();
  c1->cd(2);
 gPad->SetLeftMargin(0.15);
 frameL->GetYaxis()->SetTitleOffset(1.2);
 frameL->GetYaxis()->SetTitle(" Counts");
 frameL->GetXaxis()->SetTitle(" Decay Energy");
 frameL->Draw();
  
}