#include <iostream>
#include <fstream>
#include <iomanip> //setprecision

#include "TPaveLabel.h"
#include "TCanvas.h"
#include "TCut.h"
#include "TTree.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TColor.h"
#include "TRandom.h"
#include "TStyle.h"
#include "TPaveLabel.h"
#include "TPad.h"
#include "TMath.h"
#include "TGraphAsymmErrors.h"
#include "TGraph.h"
#include "TPaveText.h"
#include "TStopwatch.h"
#include "TDatime.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TRandom3.h"
#include "TEventList.h"
#include "TMath.h"
#include "Math/Minimizer.h"

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooPolynomial.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooExponential.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooSimultaneous.h"
#include "RooGenericPdf.h"
#include "RooExtendPdf.h"
#include "RooNumConvPdf.h"
#include "RooFFTConvPdf.h"
#include "RooProduct.h"
#include "RooFitResult.h"
#include "RooArgList.h"
#include "RooAbsCollection.h"
#include "RooPlot.h"
#include "RooHist.h"
#include "RooConstVar.h"
#include "RooCategory.h"
#include "RooBinning.h"
#include "RooJohnsonSU.h"

#include "RooChebychev.h"
#include "RooTrace.h"
#include "RooMinuit.h"
#include "MinuitParams.C"

using namespace RooFit ;
using namespace std;
using namespace ROOT::Math;


void fit() {
  //printing out values of TMinuit public members
  MinuitParams();

  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  RooTrace::active(kTRUE);
  
 RooRealVar B_VTX_mass("B_VTX_mass", "B_VTX_mass mass",  (Double_t)4930.,  (Double_t)5630.);
 RooRealVar BDT_fakeMuon_1("BDT_fakeMuon_1", "BDT_fakeMuon_1",  (Double_t)-1.,  (Double_t)1.);
 RooRealVar BDT_fakeMuon_2("BDT_fakeMuon_2", "BDT_fakeMuon_2",  (Double_t)-1.,  (Double_t)1.);


  SumW2Error(kTRUE);

  TString CUTSTRING = " BDT_fakeMuon_1>-0.435 && BDT_fakeMuon_2>-0.435 ";
  
  //Creating RooDataSet from ntuple of mass data point
  TFile *ntuple_signal_ctl = new TFile("MassMCsignal.root");
  TTree* tree_signal_ctl = (TTree*) ntuple_signal_ctl->Get("T_mv");
  RooDataSet *data_signal_ctl = new RooDataSet("data_signal_ctl", "data_signal_ctl", tree_signal_ctl, RooArgSet(B_VTX_mass,BDT_fakeMuon_1,BDT_fakeMuon_2),CUTSTRING);
  
  //Creating fit parameters
  Double_t N_sigNONRAD = data_signal_ctl->sumEntries();
  RooRealVar MC_fgaussNONRAD ("MC_fgaussNONRAD", "MC_fgaussNONRAD", (Double_t)0.06,  (Double_t)0.,  (Double_t)1.);
  
  RooRealVar   MC_nsignalNONrad_ctl("MC_nsignalNONrad_ctl", "gauss factor",                  (Double_t)N_sigNONRAD,    (Double_t)0.,  (Double_t)500000.);
  
  RooRealVar   mean   ("mean"  , "mean gauss",  (Double_t)5255.,  (Double_t)4930.,  (Double_t)5630.); //OK"point" of 1D fit convergence   (previous value 5267)
  RooRealVar   sigma  ("sigma" , "sigma gauss",  (Double_t)39.9,  (Double_t)0,  (Double_t)1000);  //OK"point" of 1D fit convergence (previous value 43.58)
  RooRealVar   meanJ        ("meanJ" , "",   (Double_t)5276.24,  (Double_t)4930.,  (Double_t)5630.);//OK
  RooRealVar   gammaJ        ("gammaJ" , "",   (Double_t)-0.19,  (Double_t)-100.,  (Double_t)0.);  //OK"point" of 1D fit convergence (previous value -0.13)
  RooRealVar   sigmaJ        ("sigmaJ" , "",   (Double_t)39.4,  (Double_t)0.,  (Double_t)1000.); //OK"point" of 1D fit convergence  (previous value 37.5)
  RooRealVar   deltaJ        ("deltaJ" , "",   (Double_t)1.8,  (Double_t)0.,  (Double_t)100.);  //OK"point" of 1D fit convergence (previous value 1.3)
  
  // summing two models into one composite PDF
  RooJohnsonSU JohnsonSU("JohnsonSU","Johnson SU distribution",B_VTX_mass,meanJ,sigmaJ,gammaJ,deltaJ);
  RooGaussian  gauss1 ("gauss1", "gaussian PDF",  B_VTX_mass, mean, sigma);

  RooAddPdf    JohnsonSUplusGauss  ("JohnsonSUplusGauss" , "JohnsonSUplusGauss for signal", RooArgList(gauss1,JohnsonSU), RooArgList(MC_fgaussNONRAD));
   
  // making the PDF extended
  RooExtendPdf eJohnsonSUplusGauss("eJohnsonSUplusGauss", "extended gauss p.d.f", JohnsonSUplusGauss, MC_nsignalNONrad_ctl);
  eJohnsonSUplusGauss.setTraceCounter(1000, kTRUE);
  eJohnsonSUplusGauss.verboseDirty(kTRUE);
  
  // running unbinned M.L. fit 
  RooAbsReal* nll = eJohnsonSUplusGauss.createNLL(*data_signal_ctl);
  RooMsgService::instance().setGlobalKillBelow(RooFit::DEBUG) ;

  RooMinuit minuit(*nll);
  minuit.setPrintLevel(5) ; // Print MINUIT messages
  minuit.setVerbose(kTRUE) ; // Print RooMinuit messages with parameter 
  minuit.setErrorLevel(5) ;
  minuit.setPrintEvalErrors(5) ;
  minuit.setStrategy(2) ;
  minuit.setWarnLevel(5) ;
  // saving the output to a log file 
  minuit.setLogFile("MinuitLog_MyMachine.log") ;
 // minuit.setLogFile("MinuitLog_LXplus.log") ;
 
  minuit.migrad();
  
  RooTrace::dump(cout,kTRUE);
 // RooFitResult * result = eJohnsonSUplusGauss.fitTo(*data_signal_ctl, Verbose(3), Save(),NumCPU(1));
 // RooTrace::dump(cout,kTRUE);
  
}