/// \file
 /// \ingroup tutorial_roofit
 /// \notebook
 ///  'SPECIAL PDFS' RooFit tutorial macro #707
 ///
 ///  Using non-parametric (multi-dimensional) kernel estimation p.d.f.s
 ///
 /// \macro_image
 /// \macro_output
 /// \macro_code
 /// \author 07/2008 - Wouter Verkerke
 
 #include "RooRealVar.h"
 #include "RooDataSet.h"
 #include "RooGaussian.h"
 #include "RooConstVar.h"
 #include "RooPolynomial.h"
 #include "RooKeysPdf.h"
 #include "RooNDKeysPdf.h"
 #include "RooProdPdf.h"
 #include "TCanvas.h"
 #include "TAxis.h"
 #include "TH1.h"
 #include "RooPlot.h"
 #include "RooJohnson.h"
 #include "TFile.h"
 using namespace RooFit ;
 
 void kde_ksks_extract()
 {
    // C r e a t e   l o w   s t a t s   1 - D   d a t a s e t
    // -------------------------------------------------------
 
    
    TFile *f = new TFile("selected.root");
    TTree *tree = (TTree*)f->Get("DstD0PiKSKS");
    //cout << "Entries = " << tree->GetEntries() << endl;
  
    RooRealVar max_sig("max_sig","", 0, 3000, "");
    RooDataSet data("data","data", tree, RooArgSet(max_sig));
    
    
    //Johnson
    RooRealVar mu_J("#mu_{J}","mean_johnson",0,0, 3000);
    RooRealVar lambda_J ("#lambda_{J}", "lambda_johnson", 0.015, 0, 0.5);
    RooRealVar gamma_J("#gamma_{J}", "gamma", 0.1,0,10);
    RooRealVar delta_J("#delta_{J}", "delta",37, -10., 10.);
    RooJohnson sig("johnson", "Johnson PDF", max_sig, mu_J, lambda_J, delta_J, gamma_J );
 
    // C r e a t e   1 - D   k e r n e l   e s t i m a t i o n   p d f
    // ---------------------------------------------------------------
 
    // Create adaptive kernel estimation pdf. In this configuration the input data
    // is mirrored over the boundaries to minimize edge effects in distribution
    // that do not fall to zero towards the edges
    RooKeysPdf kest1("kest1","kest1",max_sig,data,RooKeysPdf::MirrorBoth) ;
    // //RooNDKeysPdf ktest1("b02ggMbcDe_pdf","b02ggMbcDe_pdf", RooArgSet(Mbc,deltaE), *data,"am",0.7);
    
    // TH1F* h100 = static_cast<TH1F*>(ktest1.createHistogram("h100", max_sig, Binning(100,0,3000)));
    // TH1* hist_max_sig = (TH1*)h100;
    // RooDataHist max_sig_hist("max_sig_hist","",RooArgList(max_sig),hist_max_sig);
    // RooHistPdf ktest1_pdf("ktest1_pdf","",RooArgSet(max_sig),max_sig_hist);
 
    // An adaptive kernel estimation pdf on the same data without mirroring option
    // for comparison
    RooKeysPdf kest2("kest2","kest2",max_sig,data,RooKeysPdf::NoMirror) ;
 
   
    // Adaptive kernel estimation pdf with increased bandwidth scale factor
    // (promotes smoothness over detail preservation)
    RooKeysPdf kest3("kest1","kest1",max_sig,data,RooKeysPdf::MirrorBoth,2) ;
 
    // Plot kernel estimation pdfs with and without mirroring over data
    RooPlot* frame = max_sig.frame(Title("Adaptive kernel estimation pdf with and w/o mirroring"),Bins(100)) ;
    data.plotOn(frame) ;
    kest1.plotOn(frame) ;
    kest2.plotOn(frame,LineStyle(kDashed),LineColor(kRed)) ;
 
 
    // Plot kernel estimation pdfs with regular and increased bandwidth
    RooPlot* frame2 = max_sig.frame(Title("Adaptive kernel estimation pdf with regular, increased bandwidth")) ;
    kest1.plotOn(frame2) ;
    kest3.plotOn(frame2,LineColor(kMagenta)) ;
 
    // // C r e a t e   l o w   s t a t s   2 - D   d a t a s e t
    // // -------------------------------------------------------
    
    // RooRealVar max_sig_y("max_sig","", 0, 2000, "");
    
    // //Johnson
    // RooRealVar mu_Jy("#mu_{Jy}","mean_johnson",0,0, 2000);
    // RooRealVar lambda_Jy("#lambda_{Jy}", "lambda_johnson", 0.013, 0, 0.5);
    // RooRealVar gamma_Jy("#gamma_{Jy}", "gamma", 0.1,0,10);
    // RooRealVar delta_Jy("#delta_{Jy}", "delta",37, -10., 10.);
    // RooJohnson sigy("johnsony", "Johnson PDF", max_sig_y, mu_Jy, lambda_Jy, delta_Jy, gamma_Jy);
    
    // RooProdPdf johnsonxy("pxy","pxy",RooArgSet(sig,sigy));

    // RooDataSet data2("data2","data2", tree, RooArgSet(max_sig,max_sig_y));

    TCanvas* c = new TCanvas("rf707_kernelestimation","rf707_kernelestimation",800,800) ;
    c->Divide(1,2) ;
    c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4) ; frame->Draw() ;
    c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.8) ; frame2->Draw() ;

   // C r e a t e   w o r k s p a c e ,   i m p o r t   d a t a   a n d   m o d e l
   // -----------------------------------------------------------------------------
 
   // Create a new empty workspace
   RooWorkspace *w = new RooWorkspace("w", "workspace");
 
   // Import model and all its components into the workspace
   w->import(ktest1);
 
   // Import data into the workspace
   w->import(data);
 
   // Print workspace contents
   w->Print();
 
   // S a v e   w o r k s p a c e   i n   f i l e
   // -------------------------------------------
 
   // Save the workspace into a ROOT file
   w->writeToFile("sig_kde.root");
 
   // Workspace will remain in memory after macro finishes
   gDirectory->Add(w);
 }