#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooTruthModel.h"
#include "RooFormulaVar.h"
#include "RooRealSumPdf.h"
#include "RooPolyVar.h"
#include "RooProduct.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
  
using namespace RooFit;
void toy(){


   //if (!gROOT->GetClass("TGenPhaseSpace")) gSystem.Load("libPhysics");

   TLorentzVector W(0.0, 0.0, 0.0, 4.42);

   //(Momentum, Energy units are Gev/C, GeV)
   Double_t masses[3] = { 3.686, 0.139, 0.139} ;

   TGenPhaseSpace event;
   event.SetDecay(W, 3, masses);

  double x_min=3.8;
  double x_max=4.3;
  TFile* output = new TFile("toymc.root", "recreate");

   TH2F *h2 = new TH2F("h2","h2", 100,x_min,x_max, 100,x_min,x_max);

   for (Int_t n=0;n<100000;n++) {
      Double_t weight = event.Generate();

      TLorentzVector *Ppsip= event.GetDecay(0);

      TLorentzVector *pPip    = event.GetDecay(1);
      TLorentzVector *pPim    = event.GetDecay(2);

      TLorentzVector pPpsipPip = *Ppsip + *pPip;
      TLorentzVector pPpsipPim = *Ppsip + *pPim;

      h2->Fill(pPpsipPip.M() ,pPpsipPim.M() ,weight);
   }

   output->Write();
   output->Close();
}

using namespace RooFit;

// using intOrder = 1 generation is wrong
void draw(){

    Double_t masses[3] = { 3.686, 0.139, 0.139} ;
    double x_min=3.8;
    double x_max=4.3;

    RooRealVar  x("X1", "Y1", x_min, x_max);
    RooRealVar  y("X2", "Y2", x_min, x_max);

    RooRealVar p0("p0", "poly 0", 0.0) ;
    RooRealVar p1("p1", "poly 1", 0.0) ;

    RooPolynomial fyconst("fyconst", "fyconstent",y,p0);
    RooPolynomial fxconst("fxconst", "fxconstent",x,p1);
  
    RooProdPdf sigxy("sigxy","sigxy",fxconst,fyconst);




   TFile *f_phsppi0pi0 = new TFile("toymc.root");

   TH2F *h2 = (TH2F*)f_phsppi0pi0->Get("h2");

   RooDataHist sig2d_hist("sig2d_hist",      "sig2d_hist",  RooArgSet(x,y),        h2);
   RooHistPdf  sig2d_shape("sig2d_shape",     "sig2d_shape",      RooArgSet(x,y),        sig2d_hist , 0);


   RooProdPdf sigshpe("sigshpe","sigshpe",sigxy,sig2d_shape) ;


   RooMsgService::instance().addStream(DEBUG,Topic(Generation)) ;

   RooDataSet* data_sig = sigshpe.generate(RooArgSet(x,y),100000);


   h2->SetTitle("Original Histogram");
   h2->Draw("COLZ");

   // test drawing of histPdf
    new TCanvas();
    TH2F * hist=data_sig->createHistogram(x,y,100,100);
    hist->Draw();
}

void toymc() {
   toy();
   draw();
}