#include "TROOT.h"
#include "RooRealVar.h"
#include "RooDataSet.h"


void TestRoot(){

    using namespace RooFit;

    Int_t NBins=88;
    Double_t Dijet_Data[] = {9903888.000000, 6232449.000000, 4023546.000000, 2657443.000000, 1787047.000000, 1227732.000000, 854711.000000, 601900.000000, 431655.000000, 312226.000000, 228546.000000, 168581.000000, 125222.000000, 93555.000000, 70896.000000, 53891.000000, 41409.000000, 31530.000000, 24481.000000, 19037.000000, 14836.000000, 11729.000000, 9144.000000, 7161.000000, 5653.000000, 4615.000000, 3701.000000, 2871.000000, 2260.000000, 1825.000000, 1489.000000, 1242.000000, 937.000000, 791.000000, 614.000000, 506.000000, 450.000000, 349.000000, 287.000000, 253.000000, 185.000000, 138.000000, 116.000000, 103.000000, 80.000000, 65.000000, 50.000000, 32.000000, 34.000000, 25.000000, 17.000000, 14.000000, 10.000000, 9.000000, 3.000000, 9.000000, 4.000000, 7.000000, 2.000000, 2.000000, 2.000000, 3.000000, 4.000000, 1.000000, 2.000000, 0.000000, 1.000000, 0.000000, 1.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000};
    Double_t xmin = 1200.;
    Double_t xmax=10000.;

    // Check all OK 
    /*
    for (int i=0; i<NBins;i++){
        cout << Dijet_Data[i] << ", ";
    }
    cout << endl;
    */




    RooRealVar wt("wt","wt",1);
    RooRealVar mass("mass","M_{JJ} [GeV]",xmin,xmax);
    int nbins = TMath::Nint((xmax - xmin)/100);
    mass.setBins(nbins);

    RooDataSet d("d","d",RooArgSet(mass,wt), WeightVar(wt)) ;
    d.Print("V");


    Double_t binWidth = (xmax-xmin)/NBins;
    Double_t SumData = 0;

    for(int i =1; i < NBins+1; i++){
        Double_t center = xmin + (i-0.5)*binWidth;
        //cout << i << "\t" << center <<  "\t" << Dijet_Data[i-1] <<endl;
        wt =Dijet_Data[i-1];
        SumData += Dijet_Data[i-1];
        d.add(RooArgSet(mass,wt),Dijet_Data[i-1]);

    }
    // Printing as a cross check
    
    d.Print("V");

    RooRealVar nbg("nbg","nbg",SumData,1e6,3e7);  // background yield
    
    RooRealVar *nsig; 

    RooRealVar p2("p2","p2", 9.82979,2.,15.);
    RooRealVar p3("p3","p3", -2.97813,-8,4.);
    RooRealVar p4("p4","p4", 9.39320e-01,-5.,15.);
    RooRealVar p5("p5","p5",  1.50531e-01,-5,5.);

    RooArgList setValues(nbg,p2,p3,p4,p5);

    RooGenericPdf bgThreeParamPDF("bgThreeParam", "TMath::Power((1 - @0/13000),@1)*TMath::Power((@0/13000.),(@2+(@3*TMath::Log(@0/13000.))+(@4*TMath::Log(@0/13000.)*TMath::Log(@0/13000.))))",RooArgList(mass,p2,p3,p4,p5));
    RooAddPdf *ThreeParamBackgroundOnlyPDF  = new RooAddPdf("ThreeParamBackgroundOnlyPDF","mass fit",RooArgList(bgThreeParamPDF),RooArgList(nbg));
    RooBinSamplingPdf *integratedPDF =  new RooBinSamplingPdf("integratedPDF", "IntegratedThreeParamBackgroundOnlyPDF", mass, *ThreeParamBackgroundOnlyPDF, 1e-7);


    //RooAbsReal* nll = integratedPDF->createNLL(d,RooFit::Offset(kTRUE)); //,RooFit::Extended(kTRUE));
    RooAbsReal* nll = ThreeParamBackgroundOnlyPDF->createNLL(d,RooFit::Offset(kTRUE));
    
    nll->Print("a");
    RooMinimizer *minim = new RooMinimizer(*nll);
    minim->setStrategy(2);
    minim->Print();
    minim->setPrintLevel(1);
    minim->setEps(1e-5);
    minim->optimizeConst(2);
    int status1 = minim->minimize("Minuit2");

    RooAbsReal* nllBad = integratedPDF->createNLL(d,RooFit::Offset(kTRUE)); 
    
    nllBad->Print("a");
    RooMinimizer *minimBad = new RooMinimizer(*nll);
    minimBad->setStrategy(2);
    minimBad->Print();
    minimBad->setPrintLevel(1);
    minimBad->setEps(1e-5);
    minimBad->optimizeConst(2);
     status1 = minimBad->minimize("Minuit2");


}