// 15-10-08
//
//  use  < RooMCStudy >  to run a large number of toy MC expriments
//
// Input root file: m3.root
//
//-------------------------------------
{
    gSystem->Load("libRooFit") ;
    using namespace RooFit ;


    // input root file
    TFile *f = new TFile("m3.root");

    // Get the histograms
    TH1D *hs = f->Get("hs_1");
    TH1D *hwj = f->Get("hwj_1");
    TH1D *hzj = f->Get("hzj_1");

    double nexp_tt = hs->Integral(0,1000); // 1265.15 (100/pb)
    double nexp_wj = hwj->Integral(0,1000); // 332.411 (100/pb)
    double nexp_zj = hzj->Integral(0,1000); // 143.391 (100/pb)

    RooRealVar    m3("m3",  "m3 (GeV)",    0, 1000);
    RooDataHist  h_s("hs",  "signal data", m3, hs);
    RooDataHist h_wj("hwj", "wj data",     m3, hwj);

    
    // Build PDFs from histograms
    RooHistPdf pdf_tt("tt", "Signal pdf", m3, h_s, 2);
    RooHistPdf pdf_wj("wj", "W+jets pdf", m3, h_wj, 2);


    // declare fit parameters (unconstrained)
    RooRealVar ntt("ntt","fitted number of tt signal events",   nexp_tt, -1000 , 10000);
    RooRealVar nwj("nwj","fitted number of W+jets bgnd events", nexp_wj, -1000 , 10000);

    // create the model for data generation and fitting
    RooAddPdf model("model", "sig+wj", RooArgList(pdf_tt,pdf_wj),
		    RooArgList(ntt, nwj) ) ;
    

    // Use RooMCStudy to perform extended, unbinned maximum likelihood fit 
    //   (fit number of events)

    RooMCStudy mgr(model, m3, Silence(), 
		   FitOptions(Extended(kTRUE), PrintLevel(-1)), 
		   Extended(kTRUE) );

    
    // specify how many toy MC expriments we want to do
    const int ntoy = 10; //10000; 
    mgr.generateAndFit( ntoy  );


    // display the fit results
    RooPlot* sframe = ntt.frame(600, 2000, 50);
    mgr.plotParamOn(sframe);
    RooPlot* bframe = nwj.frame(-200, 1000, 50);
    mgr.plotParamOn(bframe);

    // --- Plot frame on canvas ---
    TCanvas c1("c1","c1",800,400);
    c1.Divide(2,1);
    c1.cd(1);    sframe->Draw() ;
    c1.cd(2);    bframe->Draw() ;
    
}