/// \file
/// \ingroup tutorial_roofit
/// \notebook -js
/// Organisation and simultaneous fits: using simultaneous p.d.f.s to describe simultaneous fits to multiple datasets
///
/// \macro_image
/// \macro_output
/// \macro_code
/// \author 07/2008 - Wouter Verkerke

#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
using namespace RooFit;

void rf501_simultaneouspdf()
{
   // C r e a t e   m o d e l   f o r   p h y s i c s   s a m p l e
   // -------------------------------------------------------------

   // Construct signal pdf
   RooRealVar x("x", "x", -8, 8);
   RooRealVar mean("mean", "mean", 0, -8, 8);
   RooRealVar sigma("sigma", "sigma", 2., 0.1, 10);
   RooGaussian gx("gx", "gx", x, mean, sigma);


   // Construct signal pdf, share sigma
   RooRealVar y("y", "y", -8, 8);
   RooRealVar meany("meany", "meany", 1, -8, 8);
   RooGaussian gy("gy", "gy", y, meany, sigma);

   // Define category to distinguish physics and control samples events
   RooCategory sample("sample", "sample");
   sample.defineType("X");
   sample.defineType("Y");

   // Construct a simultaneous pdf using category sample as index
   RooSimultaneous simPdf("simPdf", "simultaneous pdf", sample);

   // Associate model with the physics state and model_ctl with the control state
   simPdf.addPdf(gx, "X");
   simPdf.addPdf(gy, "Y");

   simPdf.Print("T");
   
   
   // G e n e r a t e   e v e n t s   f o r   b o t h   s a m p l e s
   // ---------------------------------------------------------------

   // Generate events in x and y from model
   RooDataSet *dataX = gx.generate(RooArgSet(x), 1000);
   RooDataSet *dataY = gy.generate(RooArgSet(y), 2000);

   TH1* histX = dataX->createHistogram("histX", x, Binning(100, -8, 8));
   TH1* histY = dataY->createHistogram("histY", y, Binning(100, -8, 8));

   // Construct combined dataset in (x,sample)
   RooRealVar weight("weight", "weight", 1.);
   RooArgSet datasetVars(x, y, weight, sample);
   RooDataSet combData("combData", "combined data", datasetVars, WeightVar(weight));
   
   sample.setLabel("X");
   for (int i=1; i <= histX->GetNbinsX()+1; ++i) {
      x.setVal(histX->GetBinCenter(i));
      combData.add(datasetVars, histX->GetBinContent(i));
   }
   sample.setLabel("Y");
   for (int i=1; i <= histY->GetNbinsX()+1; ++i) {
      y.setVal(histY->GetBinCenter(i));
      combData.add(datasetVars, histY->GetBinContent(i));
   }
   combData.Print("V");
   combData.get(0)->Print("V");
   //combData.get(100)->Print("V");

   // P e r f o r m   a   s i m u l t a n e o u s   f i t
   // ---------------------------------------------------

   // Perform simultaneous fit of model to data and model_ctl to data_ctl
   simPdf.fitTo(combData);

   // P l o t   m o d e l   s l i c e s   o n   d a t a    s l i c e s
   // ----------------------------------------------------------------

   // Make a frame for the physics sample
   RooPlot *frame1 = x.frame(Title("X"));

   // Plot all data tagged as physics sample
   combData.plotOn(frame1, Cut("sample==sample::X"));

   // Plot "physics" slice of simultaneous pdf.
   // NBL You _must_ project the sample index category with data using ProjWData
   // as a RooSimultaneous makes no prediction on the shape in the index category
   // and can thus not be integrated
   simPdf.plotOn(frame1, Slice(sample, "X"), ProjWData(sample, combData));

   // The same plot for the control sample slice
   RooPlot *frame2 = y.frame(Title("Y"));
   combData.plotOn(frame2, Cut("sample==sample::Y"));
   simPdf.plotOn(frame2, Slice(sample, "Y"), ProjWData(sample, combData));

   TCanvas *c = new TCanvas("rf501_simultaneouspdf", "rf403_simultaneouspdf", 800, 400);
   c->Divide(2);
   c->cd(1);
   gPad->SetLeftMargin(0.15);
   frame1->GetYaxis()->SetTitleOffset(1.4);
   frame1->Draw();
   c->cd(2);
   gPad->SetLeftMargin(0.15);
   frame2->GetYaxis()->SetTitleOffset(1.4);
   frame2->Draw();
}