/////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #110
//
// MODIFIED - 03/2013
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooGenericPdf.h"
#include "RooConstVar.h"
#include "RooAbsReal.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
using namespace RooFit ;


void rf110_normintegration_2Dmod()
{
  // S e t u p   m o d e l 
  // ---------------------
RooRealIntegral::setCacheAllNumeric(3);
  // Create observables x,y
  RooRealVar x("x","x",-10,10) ;
  RooRealVar y("y","y",-10,10) ;

  // Create p.d.f. gaussx(x,-2,3) 
  RooGaussian gx("gx","gx",x,RooConst(-2),RooConst(3)) ;
  RooGaussian gy("gy","gy",y,RooConst(-2),RooConst(3)) ;
  RooProdPdf gxy("gxy","gxy", gx, gy);
//  gxy.forceNumInt(kTRUE);

  RooGenericPdf g("g","g","exp(-((@0+2)**2+(@1+2)**2)/(2*9))", RooArgList(x,y));

  // I n t e g r a t e   n o r m a l i z e d   p d f   o v e r   s u b r a n g e
  // ----------------------------------------------------------------------------
  // Create an integral of gx_Norm[x] over x in range "signal"
  // This is the fraction of of p.d.f. gx_Norm[x] which is in the
  // range named "signal"
  RooAbsReal* igxy_sig = gxy.createIntegral(RooArgSet(x,y),NormSet(RooArgSet(x,y)),Range("thisbin")) ;
  RooAbsReal* ig_sig = g.createIntegral(RooArgSet(x,y),NormSet(RooArgSet(x,y)),Range("thisbin")) ;

  // Create histograms for displaying results
  TH2F *hist_gxy = new TH2F("hist_gxy","RooGaussian X RooGaussian",20,-10,10,20,-10,10);
  TH2F *hist_g = new TH2F("hist_g","2D RooGenericPdf",20,-10,10,20,-10,10);

  for(int i=1; i<=20; i++) {
  for(int j=1; j<=20; j++) {
    // Define a range named "signal" in x
    x.setRange("thisbin",-11+i,-10+i) ;
    y.setRange("thisbin",-11+j,-10+j) ;
  
    // Create an integral of gx_Norm[x] over x in range "signal"
    // This is the fraction of of p.d.f. gx_Norm[x] which is in the
    // range named "signal"
    // Output for user friendlyness
    double ig_sit_val = ig_sig->getVal();
    cout << "Range values: (" << -11+i << "," << -10+i << " ; " << -11+j << "," << -10+j << ")" << endl;
    cout << "   gx_Int[x,y|signal]_Norm[x,y] = " << igxy_sig->getVal() << endl ;
    cout << "   g_Int[x,y|signal]_Norm[x,y]  = " << ig_sit_val << endl ;

    // Save values to histograms
    hist_gxy->SetBinContent(i,j,igxy_sig->getVal());
    hist_g->SetBinContent(i,j,ig_sig->getVal());
  }
  }

  // Plot Histograms
  TCanvas *TmpCanv = new TCanvas("TmpCanv","TmpCanv",1200,600) ;
  TmpCanv->Divide(2,1);
  TmpCanv->cd(1);	hist_gxy->Draw("COLZ");
  TmpCanv->cd(2);	hist_g->Draw("COLZ");

}