//Dear emacs this is c++

#include <vector>
#include "TRandom3.h"
#include "RooStats/HistFactory/Measurement.h"
#include "RooStats/HistFactory/Channel.h"
#include "RooStats/HistFactory/Sample.h"
#include "RooStats/HistFactory/HistoToWorkspaceFactoryFast.h"
#include "RooStats/HistFactory/MakeModelAndMeasurementsFast.h"
#include "RooStats/ModelConfig.h"

//---------------------------------------------------------------------------
// Convert integer to string
std::string int_to_str(    int number){
  std::stringstream ss; //create a stringstream
  ss << number;    //add number to the stream
  return ss.str(); //return a string with the contents of the stream
}//end int_to_string

//---------------------------------------------------------------------------
// Do the actual fit here
void dofit(int nChan = 1){

  std::string meas_name = "testMeas";

  //Start creating your measurement 
  RooStats::HistFactory::Measurement meas( ("meas_"+meas_name).c_str() , ("meas_"+meas_name).c_str() );
    
  //Global meas parameters
  meas.SetPOI("mu");
  meas.SetLumi(1.0);
  meas.SetLumiRelErr(0.03);

  //Loop over number of channels
  for( int ichan=0; ichan<nChan; ichan++ ) { 
    
    RooStats::HistFactory::Channel chan( ("channel_"+meas_name+"_"+int_to_str(ichan)).c_str() );
    
    //-----------------------------------------------------------------------------------
    //Define signal
    RooStats::HistFactory::Sample sig( ("signal_"+int_to_str(ichan)).c_str(), ("sigHist_"+int_to_str(ichan)).c_str(), "histFile.root");
    sig.AddNormFactor( "mu"                                   , 3.0      ,  0.0      , 10.0 );
    sig.SetNormalizeByTheory( true );  
    
    chan.AddSample( sig );
    
    //-----------------------------------------------------------------------------------
    //Add the background
    RooStats::HistFactory::Sample bkg( ("background_"+int_to_str(ichan)).c_str(), ("bkgHist_"+int_to_str(ichan)).c_str(), "histFile.root");
    bkg.AddNormFactor( "bkg_mu"                                  , 1.5       , 0.0, 10.0 );
    bkg.SetNormalizeByTheory( true );
    
    chan.AddSample( bkg );
    
    meas.AddChannel( chan );
  }
    
  //-----------------------------------------------------------------------------------
  //Collect everything
  meas.CollectHistograms();
    
  //-----------------------------------------------------------------------------------    
  //Create the workspace
  RooWorkspace *ws = RooStats::HistFactory::HistoToWorkspaceFactoryFast::MakeCombinedModel( meas ); 
  //-----------------------------------------------------------------------------------    
  //Get output ... 

  //Create a vector to contain the asimov data
  std::vector<TH1*> hist;

  //Get asimov and signal / background pdfs
  RooAbsData *asimov  = ws->data("asimovData")                                 ;

  for( int ichan=0; ichan<nChan; ichan++ ) { 

    std::string funcName = int_to_str(ichan)+"_channel_testMeas_"+int_to_str(ichan)+"_nominal";

    RooAbsReal *sigFunc = ws->function( ("signal_"    +funcName).c_str() );
    RooAbsReal *bkgFunc = ws->function( ("background_"+funcName).c_str() );

    hist.push_back( asimov->createHistogram(  ("obs_x_channel_testMeas_"+int_to_str(ichan)).c_str() ) );
    hist.push_back( sigFunc->createHistogram( ("obs_x_channel_testMeas_"+int_to_str(ichan)).c_str() ) );
    hist.push_back( bkgFunc->createHistogram( ("obs_x_channel_testMeas_"+int_to_str(ichan)).c_str() ) );
  }

  TFile *tempFile = new TFile("tempFile.root", "RECREATE");
  for( std::vector<TH1*>::iterator it = hist.begin(); it != hist.end(); ++it){
    (*it)->Write();
  }
  tempFile->Write();
  tempFile->Close();
  
}//end dofit

//-----------------------------------------------------------------------------------    
//Main call here
int fittingScript( int nChan = 1 ){

  //Start by creating some pseudo data
  TH1D *sigHist_0 = new TH1D("sigHist_0", "sigHist", 100, -3.0, 3.0 );
  TH1D *bkgHist_0 = new TH1D("bkgHist_0", "bkgHist", 100, -3.0, 3.0 );

  TH1D *sigHist_1 = new TH1D("sigHist_1", "sigHist", 100, -3.0, 3.0 );
  TH1D *bkgHist_1 = new TH1D("bkgHist_1", "bkgHist", 100, -3.0, 3.0 );

  TRandom3 r;
  for( int i=0; i<10000; i++){
    //    sigHist_0->Fill( r.Gaus(0.0, 0.3    ) );
    sigHist_0->Fill( r.Uniform(-3.0, 3.0) );
    bkgHist_0->Fill( r.Uniform(-3.0, 3.0) );

    //    sigHist_1->Fill( r.Gaus(0.0, 0.15   ) );
    sigHist_1->Fill( r.Uniform(-3.0, 3.0) );
    bkgHist_1->Fill( r.Uniform(-3.0, 3.0) );
  }

  sigHist_0->Scale( 1.0 / sigHist_0->Integral());
  bkgHist_0->Scale( 1.0 / bkgHist_0->Integral());

  sigHist_1->Scale( 1.0 / sigHist_1->Integral());
  bkgHist_1->Scale( 1.0 / bkgHist_1->Integral());

  TFile *outFile = new TFile("histFile.root", "RECREATE");
  sigHist_0->Write();
  bkgHist_0->Write();

  sigHist_1->Write();
  bkgHist_1->Write();

  outFile->Write();
  outFile->Close();

  //Do fit here
  dofit( nChan );

  return 0;
}