#include <RooSimultaneous.h>
#include <TFile.h>
#include <TH1.h>
#include <TString.h>
#include <RooFitResult.h>
#include <RooRealSumPdf.h>
#include <RooWorkspace.h>
#include <RooDataSet.h>
#include <RooRealVar.h>


// The fixed version of getPropagatedError from ROOT 6.28 that also works for
// the RooRealSumPdf directly. Can be removed once ROOT 6.28 is used.
double getPropagatedError628(RooAbsReal& absReal, const RooFitResult &fr, const RooArgSet &nset={})
{
  // Calling getParameters() might be costly, but necessary to get the right
  // parameters in the RooAbsReal. The RooFitResult only stores snapshots.
  RooArgSet allParamsInAbsReal;
  absReal.getParameters(&nset, allParamsInAbsReal);

  RooArgList paramList;
  for(auto * rrvFitRes : static_range_cast<RooRealVar*>(fr.floatParsFinal())) {

     auto rrvInAbsReal = static_cast<RooRealVar const*>(allParamsInAbsReal.find(*rrvFitRes));

     // If this RooAbsReal is a RooRealVar in the fit result, we don't need to
     // propagate anything and can just return the error in the fit result
     if(rrvFitRes->namePtr() == absReal.namePtr()) return rrvFitRes->getError();

     // Strip out parameters with zero error
     if (rrvFitRes->getError() <= rrvFitRes->getVal() * std::numeric_limits<double>::epsilon()) continue;

     // Ignore parameters in the fit result that this RooAbsReal doesn't depend on
     if(!rrvInAbsReal) continue;

     // Checking for float equality is a bad. We check if the values are
     // negligibly far away from each other, relative to the uncertainty.
     if(std::abs(rrvInAbsReal->getVal() - rrvFitRes->getVal()) > 0.01 * rrvFitRes->getError()) {
        std::stringstream errMsg;
        errMsg << "RooAbsReal::getPropagatedError(): the parameters of the RooAbsReal don't have"
               << " the same values as in the fit result! The logic of getPropagatedError is broken in this case.";

        throw std::runtime_error(errMsg.str());
     }

     paramList.add(*rrvInAbsReal);
  }

  std::vector<double> plusVar;
  std::vector<double> minusVar;
  plusVar.reserve(paramList.size());
  minusVar.reserve(paramList.size());

  // Create std::vector of plus,minus variations for each parameter
  TMatrixDSym V(paramList.size() == fr.floatParsFinal().size() ?
      fr.covarianceMatrix() :
      fr.reducedCovarianceMatrix(paramList)) ;

  for (Int_t ivar=0 ; ivar<paramList.getSize() ; ivar++) {

    auto& rrv = static_cast<RooRealVar&>(paramList[ivar]);

    double cenVal = rrv.getVal() ;
    double errVal = sqrt(V(ivar,ivar)) ;

    // Make Plus variation
    rrv.setVal(cenVal+errVal) ;
    plusVar.push_back(absReal.getVal(nset)) ;

    // Make Minus variation
    rrv.setVal(cenVal-errVal) ;
    minusVar.push_back(absReal.getVal(nset)) ;

    rrv.setVal(cenVal) ;
  }

  // Re-evaluate this RooAbsReal with the central parameters just to be
  // extra-safe that a call to `getPropagatedError()` doesn't change any state.
  // It should not be necessarry because thanks to the dirty flag propagation
  // the RooAbsReal is re-evaluated anyway the next time getVal() is called.
  // Still there are imaginable corner cases where it would not be triggered,
  // for example if the user changes the RooFit operation more after the error
  // propagation.
  absReal.getVal(nset);

  TMatrixDSym C(paramList.getSize()) ;
  std::vector<double> errVec(paramList.getSize()) ;
  for (int i=0 ; i<paramList.getSize() ; i++) {
    errVec[i] = std::sqrt(V(i,i)) ;
    for (int j=i ; j<paramList.getSize() ; j++) {
      C(i,j) = V(i,j) / std::sqrt(V(i,i)*V(j,j));
      C(j,i) = C(i,j) ;
    }
  }

  // Make std::vector of variations
  TVectorD F(plusVar.size()) ;
  for (unsigned int j=0 ; j<plusVar.size() ; j++) {
    F[j] = (plusVar[j]-minusVar[j])/2 ;
  }

  // Calculate error in linear approximation from variations and correlation coefficient
  double sum = F*(C*F) ;

  return sqrt(sum) ;
}

// Util::ComponentToHistogram
void testJonas(){
    // Load everything
    TFile *tfile( TFile::Open("test626.root", "r") );

    RooWorkspace *w = tfile->Get<RooWorkspace>("w");


    RooSimultaneous *m_pdf = static_cast<RooSimultaneous*>(w->pdf("simPdf"));

    TString regionCategoryLabel("UserRegion_cuts");
    RooAbsPdf *m_regionPdf = static_cast<RooAbsPdf*>(m_pdf->getPdf(regionCategoryLabel.Data()));

    RooAbsData* obsData = static_cast<RooAbsData*>(w->data("obsData"));

    RooDataSet* regionData = static_cast<RooDataSet*>(obsData->reduce("channelCat==channelCat::UserRegion_cuts"));
    
    TString m_regionVariableName("obs_x_UserRegion_cuts");
    RooRealVar *variable = static_cast<RooRealVar*>( static_cast<RooArgSet*>(m_regionPdf->getObservables(*regionData))->find(m_regionVariableName));
    
    RooRealSumPdf* RRSPdf = (RooRealSumPdf*)w->pdf("UserRegion_cuts_model");
    
    RooFitResult* fitResult = static_cast<RooFitResult*>( w->genobj("RooExpandedFitResult_beforeFit") );
    w->loadSnapshot("snapshot_paramsVals_RooExpandedFitResult_beforeFit");

    // Build a TH1-based histogram from a pdf, an observable and a fit result

    // Build a histogram according to the binning of the variable
    auto hist = RRSPdf->createHistogram(Form("hist_%s", RRSPdf->GetName()), *variable);

    // Remember original variable value to reset later
    const double origVal = variable->getVal();

    // Now loop over the bins and fill them.
    // Skip the under- and overflow bins because they don't exist in RooFit.
    for(int i=1; i < hist->GetNbinsX()+1; ++i) {
        variable->setBin(i-1);

        // Now get the value and the error
        auto sum = RRSPdf->getVal();

        // In ROOT 6.28, we can direcly use RooAbsReal::getPropagatedError()
        // auto err = RRSPdf->getPropagatedError(*fitResult);
        auto err = getPropagatedError628(*RRSPdf, *fitResult);

        // and put them in the histogram
        hist->SetBinContent(i, sum);
        hist->SetBinError(i, err);
    }

    // Reset original variable value
    variable->setVal(origVal);

    hist->Print("all");
    
}