/////////////////////////////////////////////////////////////////////////
//
// '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_1Dmod()
{
  // S e t u p   m o d e l 
  // ---------------------

  // Create observables x,y
  RooRealVar x("x","x",-10,10) ;

  // Create p.d.f. gaussx(x,-2,3) 
  RooGaussian gx("gx","gx",x,RooConst(-2),RooConst(3)) ;
  RooGenericPdf g("g","g","exp(-(x+2)**2/(2*9))", RooArgList(x));

  // 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* igx_sig = gx.createIntegral(RooArgSet(x),NormSet(x),Range("thisbin")) ;
  RooAbsReal* ig_sig  = g.createIntegral(RooArgSet(x),NormSet(x),Range("thisbin")) ;

  // Create histograms for displaying results
  TH1F *hist_gx = new TH1F("hist_gx","RooGaussian",20,-10,10);
  TH1F *hist_g = new TH1F("hist_g","1D RooGenericPdf",20,-10,10);

  for(int i=1; i<=20; i++) {
     // Define a range named "signal" in x
     x.setRange("thisbin",-11+i,-10+i) ;
  
     // Output for user friendlyness
     double ig_sig_val = ig_sig->getVal();
     cout << "Range: (" << i << ", " << i+1 << ")" << endl;
     cout << "  gx_Int[x|signal]_Norm[x] = " << igx_sig->getVal() << endl ;
     cout << "  g_Int[x|signal]_Norm[x] = " << ig_sig_val << endl ;

     // Save values to histograms
     hist_gx->SetBinContent(i,igx_sig->getVal());
     hist_g->SetBinContent(i,ig_sig_val);
  }

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

}