// @(#)root/roostats:$Id$
// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
/*************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

//_________________________________________________
/*
BEGIN_HTML
<p>
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation 
of a credible interval using a Bayesian method. 
The class works only for one single parameter of interest and it integrates the likelihood function with the given prior
probability density function to compute the posterior probability. The result of the class is a one dimensional interval 
(class SimpleInterval ), which is obtained from inverting the cumulative posterior distribution. 
</p>

<p>
The interface allows one to construct the class by passing the data set, probability density function for the model, the prior
functions and then the parameter of interest to scan. The nuisance parameters can also be passed to be marginalized when 
computing the posterior. Alternatively, the class can be constructed by passing the data and the ModelConfig containing
all the needed information (model pdf, prior pdf, parameter of interest, nuisance parameters, etc..)
</p>
<p>
After configuring the calculator, one only needs to ask GetInterval(), which
will return an SimpleInterval object. By default the extrem of the integral are obtained by inverting directly the
cumulative posterior distribution. By using the method <tt>SetScanOfPosterior(nbins)</tt> the interval is then obtained by 
scanning  the posterior function in the given number of points. The firts method is in general faster but it requires an
integration one extra dimension  ( in the poi in addition to the nuisance parameters), therefore in some case it can be
less robust.   
</p>
The class can also return the posterior function (method <tt>GetPosteriorFunction</tt>) or if needed the normalized
posterior function (the posterior pdf) (method <tt>GetPosteriorPdf</tt>). A posterior plot is also obtained using 
the <tt>GetPosteriorPlot</tt> method.
<p>
The class allows to use different integration methods for integrating in the nuisances and in the poi. All the numerical 
integration methods of ROOT can be used via the method <tt>SetIntegrationType</tt> (see more in the documentation of
this method).

END_HTML
*/
//_________________________________________________


// include other header files

#include "RooAbsFunc.h"
#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooArgSet.h"
#include "RooBrentRootFinder.h"
#include "RooFormulaVar.h"
#include "RooGenericPdf.h"
#include "RooPlot.h"
#include "RooProdPdf.h"
#include "RooDataSet.h"

// include header file of this class 
#include "RooStats/BayesianCalculator.h"
#include "RooStats/ModelConfig.h"
#include "RooStats/RooStatsUtils.h"

#include "Math/IFunction.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/Integrator.h"
#include "Math/RootFinder.h"
#include "Math/BrentMinimizer1D.h"
#include "RooFunctor.h"
#include "RooFunctor1DBinding.h"
#include "RooTFnBinding.h"
#include "RooMsgService.h"

#include "TAxis.h"
#include "TF1.h"
#include "TH1.h"
#include "TMath.h"
#include "TCanvas.h"

#include <map>
#include <cmath>

//#include "TRandom.h"
#include "RConfigure.h"

ClassImp(RooStats::BayesianCalculator)

using namespace std;

namespace RooStats { 


// first some utility classes and functions 

#ifdef R__HAS_MATHMORE
   const ROOT::Math::RootFinder::EType kRootFinderType = ROOT::Math::RootFinder::kGSL_BRENT; 
#else 
   const ROOT::Math::RootFinder::EType kRootFinderType = ROOT::Math::RootFinder::kBRENT; 
#endif




struct  LikelihoodFunction { 
   LikelihoodFunction(RooFunctor & f, RooFunctor * prior = 0, double offset = 0) : 
      fFunc(f), fPrior(prior),
      fOffset(offset), fMaxL(0) {
      fFunc.binding().resetNumCall(); 
   }

   void SetPrior(RooFunctor * prior) { fPrior = prior; }

   double operator() (const double *x ) const { 
      double nll = fFunc(x) - fOffset;
      double likelihood =  std::exp(-nll);

      if (fPrior) likelihood *= (*fPrior)(x);

      int nCalls = fFunc.binding().numCall();
      if (nCalls > 0 && nCalls % 1000 == 0) { 
         ooccoutD((TObject*)0,Eval) << "Likelihood evaluation ncalls = " << nCalls
                                    << " x0 " << x[0] << "  nll = " << nll+fOffset;
         if (fPrior) ooccoutD((TObject*)0,Eval) << " prior(x) = " << (*fPrior)(x); 
         ooccoutD((TObject*)0,Eval) << " likelihood " << likelihood                                     
                                    << " max Likelihood " << fMaxL << std::endl;                                            
      }
 
      if  (likelihood > fMaxL ) { 
         fMaxL = likelihood; 
         if ( likelihood > 1.E10) { 
            ooccoutW((TObject*)0,Eval) << "LikelihoodFunction::()  WARNING - Huge likelihood value found for  parameters ";
            for (int i = 0; i < fFunc.nObs(); ++i) 
               ooccoutW((TObject*)0,Eval) << " x[" << i << " ] = " << x[i]; 
            ooccoutW((TObject*)0,Eval) << "  nll = " << nll << " L = " << likelihood << std::endl;
         }
      }

      return likelihood; 
   }

   // for the 1D case
   double operator() (double x) const { 
      // just call the previous method
      assert(fFunc.nObs() == 1); // check nobs = 1
      double tmp = x; 
      return (*this)(&tmp); 
   }

   RooFunctor & fFunc;     // functor representing the nll function 
   RooFunctor * fPrior;     // functor representing the prior function 
   double fOffset;         //  offset used to bring the nll in a reasanble range for computing the exponent
   mutable double fMaxL;
};


// Posterior CDF Function class 
// for integral of posterior function in nuisance and POI
// 1-Dim function as function of the poi 

class PosteriorCdfFunction : public ROOT::Math::IGenFunction { 

public:

   PosteriorCdfFunction(RooAbsReal & nll,  RooArgList & bindParams, RooAbsReal * prior = 0, const char * integType = 0, double nllMinimum = 0) : 
      fFunctor(nll, bindParams, RooArgList() ),              // functor 
      fPriorFunc(std::shared_ptr<RooFunctor>((RooFunctor*)0)),
      fLikelihood(fFunctor, 0, nllMinimum),         // integral of exp(-nll) function
      fIntegrator(ROOT::Math::IntegratorMultiDim::GetType(integType) ),  // integrator 
      fXmin(bindParams.getSize() ),               // vector of parameters (min values) 
      fXmax(bindParams.getSize() ),               // vector of parameter (max values) 
      fNorm(1.0), fNormErr(0.0), fOffset(0), fMaxPOI(0), 
      fHasNorm(false),  fUseOldValues(true), fError(false) 
   {     

      if (prior) { 
         fPriorFunc = std::shared_ptr<RooFunctor>(new RooFunctor(*prior, bindParams, RooArgList() ));
         fLikelihood.SetPrior(fPriorFunc.get() );
      }

      fIntegrator.SetFunction(fLikelihood, bindParams.getSize() );
         
      ooccoutD((TObject*)0,NumIntegration) << "PosteriorCdfFunction::Compute integral of posterior in nuisance and poi. " 
                                           << " nllMinimum is " << nllMinimum << std::endl; 
   
      std::vector<double> par(bindParams.getSize());
      for (unsigned int i = 0; i < fXmin.size(); ++i) { 
         RooRealVar & var = (RooRealVar &) bindParams[i]; 
         fXmin[i] = var.getMin(); 
         fXmax[i] = var.getMax();
         par[i] = var.getVal();
         ooccoutD((TObject*)0,NumIntegration) << "PosteriorFunction::Integrate" << var.GetName() 
                                              << " in interval [ " <<  fXmin[i] << " , " << fXmax[i] << " ] " << std::endl;
      } 

      fIntegrator.Options().Print(ooccoutD((TObject*)0,NumIntegration));

      // store max POI value because it will be changed when evaluating the function 
      fMaxPOI = fXmax[0];

      // compute first the normalization with  the poi 
      fNorm = (*this)( fMaxPOI );  
      if (fError) ooccoutE((TObject*)0,NumIntegration) << "PosteriorFunction::Error computing normalization - norm = " << fNorm << std::endl;
      fHasNorm = true; 
      fNormCdfValues.insert(std::make_pair(fXmin[0], 0) );
      fNormCdfValues.insert(std::make_pair(fXmax[0], 1.0) );

   }
   
   // copy constructor (needed for Cloning the object)
   // need special treatment because integrator 
   // has no copy constuctor
   PosteriorCdfFunction(const PosteriorCdfFunction & rhs) : 
      ROOT::Math::IGenFunction(),
      fFunctor(rhs.fFunctor),
      //fPriorFunc(std::shared_ptr<RooFunctor>((RooFunctor*)0)),
      fPriorFunc(rhs.fPriorFunc),
      fLikelihood(fFunctor, fPriorFunc.get(), rhs.fLikelihood.fOffset),  
      fIntegrator(ROOT::Math::IntegratorMultiDim::GetType( rhs.fIntegrator.Name().c_str() ) ),  // integrator 
      fXmin( rhs.fXmin), 
      fXmax( rhs.fXmax), 
      fNorm( rhs.fNorm), 
      fNormErr( rhs.fNormErr), 
      fOffset(rhs.fOffset), 
      fMaxPOI(rhs.fMaxPOI), 
      fHasNorm(rhs.fHasNorm), 
      fUseOldValues(rhs.fUseOldValues), 
      fError(rhs.fError), 
      fNormCdfValues(rhs.fNormCdfValues)
   { 
      fIntegrator.SetFunction(fLikelihood, fXmin.size() );
      // need special treatment for the smart pointer
      // if (rhs.fPriorFunc.get() ) { 
      //    fPriorFunc = std::shared_ptr<RooFunctor>(new RooFunctor(*(rhs.fPriorFunc) ) );
      //    fLikelihood.SetPrior( fPriorFunc.get() );
      // }
   }
                                   
    
   bool HasError() const { return fError; }


   ROOT::Math::IGenFunction * Clone() const { 
      ooccoutD((TObject*)0,NumIntegration) << " cloning function .........." << std::endl;
      return new PosteriorCdfFunction(*this); 
   }

   // offset value for computing the root
   void SetOffset(double offset) { fOffset = offset; }

private:

   // make assignment operator private
   PosteriorCdfFunction& operator=(const PosteriorCdfFunction &) { 
      return *this; 
   }

   double DoEval (double x) const {
 
      // evaluate cdf at poi value x by integrating poi from [xmin,x] and all the nuisances  
      fXmax[0] = x;
      if (x <= fXmin[0] ) return -fOffset; 
      // could also avoid a function evaluation at maximum
      if (x >= fMaxPOI && fHasNorm) return 1. - fOffset;  // cdf is bound to these values

      // computes the integral using a previous cdf estimate
      double  normcdf0 = 0; 
      if (fHasNorm && fUseOldValues) { 
         // look in the map of the stored cdf values the closes one
         std::map<double,double>::iterator itr = fNormCdfValues.upper_bound(x); 
         itr--;   // upper bound returns a poistion 1 up of the value we want
         if (itr != fNormCdfValues.end() ) { 
            fXmin[0] = itr->first; 
            normcdf0 = itr->second;
            // ooccoutD((TObject*)0,NumIntegration) << "PosteriorCdfFunction:   computing integral between in poi interval : " 
            //                                      << fXmin[0] << " -  " << fXmax[0] << std::endl; 
         }
      }

      fFunctor.binding().resetNumCall();  // reset number of calls for debug

      double cdf = fIntegrator.Integral(&fXmin[0],&fXmax[0]);  
      double error = fIntegrator.Error(); 
      double normcdf =  cdf/fNorm;  // normalize the cdf 

      ooccoutD((TObject*)0,NumIntegration) << "PosteriorCdfFunction: poi = [" << fXmin[0] << " , " 
                                           << fXmax[0] << "] integral =  " << cdf << " +/- " << error 
                                           << "  norm-integ = " << normcdf << " cdf(x) = " << normcdf+normcdf0 
                                           << " ncalls = " << fFunctor.binding().numCall() << std::endl; 

      if (TMath::IsNaN(cdf) || cdf > std::numeric_limits<double>::max()) { 
         ooccoutE((TObject*)0,NumIntegration) << "PosteriorFunction::Error computing integral - cdf = " 
                                              << cdf << std::endl;
         fError = true; 
      }

      if (cdf != 0 && error/cdf > 0.2 ) 
         oocoutW((TObject*)0,NumIntegration) << "PosteriorCdfFunction: integration error  is larger than 20 %   x0 = " << fXmin[0]  
                                              << " x = " << x << " cdf(x) = " << cdf << " +/- " << error << std::endl;

      if (!fHasNorm) { 
         oocoutI((TObject*)0,NumIntegration) << "PosteriorCdfFunction - integral of posterior = " 
                                             << cdf << " +/- " << error << std::endl; 
         fNormErr = error;
         return cdf; 
      }

      normcdf += normcdf0;

      // store values in the map
      if (fUseOldValues) { 
         fNormCdfValues.insert(std::make_pair(x, normcdf) );
      }

      double errnorm = sqrt( error*error + normcdf*normcdf * fNormErr * fNormErr )/fNorm;  
      if (normcdf > 1. + 3 * errnorm) {
         oocoutW((TObject*)0,NumIntegration) << "PosteriorCdfFunction: normalized cdf values is larger than 1" 
                                              << " x = " << x << " normcdf(x) = " << normcdf << " +/- " << error/fNorm << std::endl;
      }

      return normcdf - fOffset;  // apply an offset (for finding the roots) 
   }

   mutable RooFunctor fFunctor;                   // functor binding nll 
   mutable std::shared_ptr<RooFunctor> fPriorFunc;  // functor binding the prior 
   LikelihoodFunction fLikelihood;               // likelihood function
   mutable ROOT::Math::IntegratorMultiDim  fIntegrator; // integrator  (mutable because Integral() is not const
   mutable std::vector<double> fXmin;    // min value of parameters (poi+nuis) - 
   mutable std::vector<double> fXmax;   // max value of parameters (poi+nuis) - max poi changes so it is mutable
   double fNorm;      // normalization value (computed in ctor) 
   mutable double fNormErr;    // normalization error value (computed in ctor) 
   double fOffset;   // offset for computing the root 
   double fMaxPOI;  // maximum value of POI 
   bool fHasNorm; // flag to control first call to the function 
   bool fUseOldValues;  // use old cdf values
   mutable bool fError;     // flag to indicate if a numerical evaluation error occurred
   mutable std::map<double,double> fNormCdfValues; 
};

//__________________________________________________________________
// Posterior Function class 
// 1-Dim function as function of the poi 
// and it integrated all the nuisance parameters 

class PosteriorFunction : public ROOT::Math::IGenFunction { 

public: 


   PosteriorFunction(RooAbsReal & nll, RooRealVar & poi, RooArgList & nuisParams, RooAbsReal * prior = 0, const char * integType = 0, double
                     norm = 1.0,  double nllOffset = 0, int niter = 0) :
      fFunctor(nll, nuisParams, RooArgList() ),
      fPriorFunc(std::shared_ptr<RooFunctor>((RooFunctor*)0)),
      fLikelihood(fFunctor, 0, nllOffset), 
      fPoi(&poi),
      fXmin(nuisParams.getSize() ),
      fXmax(nuisParams.getSize() ), 
      fNorm(norm),
      fError(0)
   { 

      if (prior) { 
         fPriorFunc = std::shared_ptr<RooFunctor>(new RooFunctor(*prior, nuisParams, RooArgList() ));
         fLikelihood.SetPrior(fPriorFunc.get() );
      }

      ooccoutD((TObject*)0,NumIntegration) << "PosteriorFunction::Evaluate the posterior function by integrating the nuisances: " << std::endl;
      for (unsigned int i = 0; i < fXmin.size(); ++i) { 
         RooRealVar & var = (RooRealVar &) nuisParams[i]; 
         fXmin[i] = var.getMin(); 
         fXmax[i] = var.getMax();
         ooccoutD((TObject*)0,NumIntegration) << "PosteriorFunction::Integrate " << var.GetName() 
                                              << " in interval [" <<  fXmin[i] << " , " << fXmax[i] << " ] " << std::endl;
      }
      if (fXmin.size() == 1) { // 1D case  
         fIntegratorOneDim.reset( new ROOT::Math::Integrator(ROOT::Math::IntegratorOneDim::GetType(integType) ) );
                                          
         fIntegratorOneDim->SetFunction(fLikelihood);
         // interested only in relative tolerance
         //fIntegratorOneDim->SetAbsTolerance(1.E-300);
         fIntegratorOneDim->Options().Print(ooccoutD((TObject*)0,NumIntegration) );
      }
      else if (fXmin.size() > 1) { // multiDim case          
         fIntegratorMultiDim.reset(new ROOT::Math::IntegratorMultiDim(ROOT::Math::IntegratorMultiDim::GetType(integType) ) );
         fIntegratorMultiDim->SetFunction(fLikelihood, fXmin.size());
         ROOT::Math::IntegratorMultiDimOptions opt = fIntegratorMultiDim->Options();
         if (niter > 0) { 
            opt.SetNCalls(niter); 
            fIntegratorMultiDim->SetOptions(opt);
         }
         //fIntegratorMultiDim->SetAbsTolerance(1.E-300);
         // print the options
         opt.Print(ooccoutD((TObject*)0,NumIntegration) );
      }
   }
      
      
   ROOT::Math::IGenFunction * Clone() const { 
      assert(1); 
      return 0; // cannot clone this function for integrator 
   } 

   double Error() const { return fError;}
   
   
private: 
   double DoEval (double x) const { 

      // evaluate posterior function at a poi value x by integrating all nuisance parameters

      fPoi->setVal(x);
      fFunctor.binding().resetNumCall();  // reset number of calls for debug

      double f = 0; 
      double error = 0; 
      if (fXmin.size() == 1) { // 1D case  
         f = fIntegratorOneDim->Integral(fXmin[0],fXmax[0]); 
         error = fIntegratorOneDim->Error();
      }
      else if (fXmin.size() > 1) { // multi-dim case
         f = fIntegratorMultiDim->Integral(&fXmin[0],&fXmax[0]); 
         error = fIntegratorMultiDim->Error();
      } else { 
         // no integration to be done
         f = fLikelihood(x);
      }

      // debug 
      ooccoutD((TObject*)0,NumIntegration) << "PosteriorFunction:  POI value  =  " 
                                           << x << "\tf(x) =  " << f << " +/- " << error 
                                           << "  norm-f(x) = " << f/fNorm 
                                           << " ncalls = " << fFunctor.binding().numCall() << std::endl; 




      if (f != 0 && error/f > 0.2 ) 
         ooccoutW((TObject*)0,NumIntegration) << "PosteriorFunction::DoEval - Error from integration in " 
                                              << fXmin.size() <<  " Dim is larger than 20 % " 
                                              << "x = " << x << " p(x) = " << f << " +/- " << error << std::endl;

      fError = error / fNorm; 
      return f / fNorm;
   }

   mutable RooFunctor fFunctor; 
   mutable std::shared_ptr<RooFunctor> fPriorFunc;  // functor binding the prior 
   LikelihoodFunction fLikelihood; 
   RooRealVar * fPoi;
   std::unique_ptr<ROOT::Math::Integrator>  fIntegratorOneDim; 
   std::unique_ptr<ROOT::Math::IntegratorMultiDim>  fIntegratorMultiDim; 
   std::vector<double> fXmin; 
   std::vector<double> fXmax; 
   double fNorm;
   mutable double fError;
};

//////////////////////////////////////////////////////////////////////////////////////////////////////////
// Posterior function obtaining sampling  toy MC for th enuisance according to their pdf  
class PosteriorFunctionFromToyMC : public ROOT::Math::IGenFunction { 

public: 


   PosteriorFunctionFromToyMC(RooAbsReal & nll, RooAbsPdf & pdf, RooRealVar & poi, RooArgList & nuisParams, RooAbsReal * prior = 0, double
                              nllOffset = 0, int niter = 0, bool redoToys = true ) :
      fFunctor(nll, nuisParams, RooArgList() ),
      fPriorFunc(std::shared_ptr<RooFunctor>((RooFunctor*)0)),
      fLikelihood(fFunctor, 0, nllOffset), 
      fPdf(&pdf),
      fPoi(&poi),
      fNuisParams(nuisParams),
      fGenParams(0),
      fNumIterations(niter),
      fError(-1),
      fRedoToys(redoToys)
   { 
      if (niter == 0) fNumIterations = 100; // default value 

      if (prior) { 
         fPriorFunc = std::shared_ptr<RooFunctor>(new RooFunctor(*prior, nuisParams, RooArgList() ));
         fLikelihood.SetPrior(fPriorFunc.get() );
      }

      ooccoutI((TObject*)0,InputArguments) << "PosteriorFunctionFromToyMC::Evaluate the posterior function by randomizing the nuisances:  niter " << fNumIterations << std::endl;

      ooccoutI((TObject*)0,InputArguments) << "PosteriorFunctionFromToyMC::Pdf used for randomizing the nuisance is " << fPdf->GetName() << std::endl; 
      // check that pdf contains  the nuisance 
      RooArgSet * vars = fPdf->getVariables(); 
      for (int i = 0; i < fNuisParams.getSize(); ++i) { 
         if (!vars->find( fNuisParams[i].GetName() ) ) { 
            ooccoutW((TObject*)0,InputArguments) << "Nuisance parameter " << fNuisParams[i].GetName() 
                                                 << " is not part of sampling pdf. " 
                                                 << "they will be trated as constant " << std::endl;
         }
      }
      delete vars;

      if (!fRedoToys) { 
         ooccoutI((TObject*)0,InputArguments) << "PosteriorFunctionFromToyMC::Generate nuisance toys only one time (for all POI points)" << std::endl; 
         GenerateToys();
      }
   }

   virtual ~PosteriorFunctionFromToyMC() { if (fGenParams) delete fGenParams; }

   // generate first n-samples of the nuisance parameters 
   void GenerateToys() const {    
      if (fGenParams) delete fGenParams;
      fGenParams = fPdf->generate(fNuisParams, fNumIterations);
      if(fGenParams==0) {
         ooccoutE((TObject*)0,InputArguments) << "PosteriorFunctionFromToyMC - failed to generate nuisance parameters" << std::endl;
      }
   }

   double Error() const { return fError;}

   ROOT::Math::IGenFunction * Clone() const { 
      // use defsult copy constructor 
      //return new PosteriorFunctionFromToyMC(*this);
      //  clone not implemented  
      assert(1);
      return 0;
   }

private:
   // evaluate the posterior at the poi value x 
   double DoEval( double x) const { 

      int npar = fNuisParams.getSize();
      assert (npar > 0);  

      
      // generate the toys 
      if (fRedoToys) GenerateToys();
      if (!fGenParams) return 0;

      // evaluate posterior function at a poi value x by integrating all nuisance parameters

      fPoi->setVal(x);

      // loop over all of the generate data 
      double sum = 0; 
      double sum2 = 0; 

      for(int iter=0; iter<fNumIterations; ++iter) {

         // get the set of generated parameters and set the nuisance parameters to the generated values
         std::vector<double> p(npar); 
         for (int i = 0; i < npar; ++i) { 
            const RooArgSet* genset=fGenParams->get(iter);
            RooAbsArg * arg = genset->find( fNuisParams[i].GetName() );
            RooRealVar * var = dynamic_cast<RooRealVar*>(arg);
            assert(var != 0);
            p[i] = var->getVal();
            ((RooRealVar &) fNuisParams[i]).setVal(p[i]);
         }

         // evaluate now the likelihood function 
         double fval =  fLikelihood( &p.front() );

         // liklihood already must contained the pdf we have sampled 
         // so we must divided by it. The value must be normalized on all 
         // other parameters 
         RooArgSet arg(fNuisParams);
         double nuisPdfVal = fPdf->getVal(&arg); 
         fval /= nuisPdfVal;


         if( fval > std::numeric_limits<double>::max()  ) {
            ooccoutE((TObject*)0,Eval) <<  "BayesianCalculator::EvalPosteriorFunctionFromToy : " 
                        << "Likelihood evaluates to infinity " << std::endl;
            ooccoutE((TObject*)0,Eval) <<  "poi value =  " << x << std::endl; 
            ooccoutE((TObject*)0,Eval) <<  "Nuisance  parameter values :  ";
            for (int i = 0; i < npar; ++i)  
               ooccoutE((TObject*)0,Eval) << fNuisParams[i].GetName() << " = " << p[i] << " ";
            ooccoutE((TObject*)0,Eval) <<  " - return 0   " << std::endl;

            fError = 1.E30; 
            return 0;
         }
         if(  TMath::IsNaN(fval) ) {
            ooccoutE((TObject*)0,Eval) <<  "BayesianCalculator::EvalPosteriorFunctionFromToy : " 
                        << "Likelihood is a NaN " << std::endl;
            ooccoutE((TObject*)0,Eval) <<  "poi value =  " << x << std::endl; 
            ooccoutE((TObject*)0,Eval) <<  "Nuisance  parameter values :  ";
            for (int i = 0; i < npar; ++i)  
               ooccoutE((TObject*)0,Eval) << fNuisParams[i].GetName() << " = " << p[i] << " ";
            ooccoutE((TObject*)0,Eval) <<  " - return 0   " << std::endl;
            fError = 1.E30; 
            return 0;
         }


         
         sum += fval; 
         sum2 += fval*fval; 
      }
   
      // compute the average and variance 
      double val = sum/double(fNumIterations);
      double dval2 = std::max( sum2/double(fNumIterations) - val*val, 0.0);
      fError = std::sqrt( dval2 / fNumIterations);

      // debug 
      ooccoutD((TObject*)0,NumIntegration) << "PosteriorFunctionFromToyMC:  POI value  =  " 
                                           << x << "\tp(x) =  " << val << " +/- " << fError << std::endl;
      

      if (val != 0 && fError/val > 0.2 ) {
         ooccoutW((TObject*)0,NumIntegration) << "PosteriorFunctionFromToyMC::DoEval" 
                                              << " - Error in estimating posterior is larger than 20% ! " 
                                              << "x = " << x << " p(x) = " << val << " +/- " << fError << std::endl;
      }


      return val; 
   }

   mutable RooFunctor fFunctor; 
   mutable std::shared_ptr<RooFunctor> fPriorFunc;  // functor binding the prior 
   LikelihoodFunction fLikelihood; 
   mutable RooAbsPdf * fPdf;
   RooRealVar * fPoi;
   RooArgList fNuisParams;
   mutable RooDataSet * fGenParams;
   int fNumIterations;
   mutable double fError; 
   bool fRedoToys;                    // do toys every iteration

};

////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Implementation of BayesianCalculator 
/////////////////////////////////////////////////////////////////////////////////////////////////////

BayesianCalculator::BayesianCalculator() :
   fData(0),
   fPdf(0),
   fPriorPdf(0),
   fNuisancePdf(0),
   fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0), 
   fPosteriorFunction(0), fApproxPosterior(0),
   fLower(0), fUpper(0),
   fNLLMin(0),
   fSize(0.05), fLeftSideFraction(0.5), 
   fBrfPrecision(0.00005), 
   fNScanBins(-1),
   fNumIterations(0),
   fValidInterval(false)
{
   // default constructor
}

BayesianCalculator::BayesianCalculator( /* const char* name,  const char* title, */						   
						    RooAbsData& data,
                                                    RooAbsPdf& pdf,
						    const RooArgSet& POI,
						    RooAbsPdf& priorPdf,
						    const RooArgSet* nuisanceParameters ) :
   fData(&data),
   fPdf(&pdf),
   fPOI(POI),
   fPriorPdf(&priorPdf),
   fNuisancePdf(0),
   fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0),
   fPosteriorFunction(0), fApproxPosterior(0),
   fLower(0), fUpper(0), 
   fNLLMin(0),
   fSize(0.05), fLeftSideFraction(0.5), 
   fBrfPrecision(0.00005), 
   fNScanBins(-1),
   fNumIterations(0),
   fValidInterval(false)
{
   // Constructor from data set, model pdf, parameter of interests and prior pdf
   // If nuisance parameters are given they will be integrated according either to the prior or 
   // their constraint term included in the model

   if (nuisanceParameters) fNuisanceParameters.add(*nuisanceParameters); 
   // remove constant nuisance parameters 
   RemoveConstantParameters(&fNuisanceParameters);
}

BayesianCalculator::BayesianCalculator( RooAbsData& data,
                       ModelConfig & model) : 
   fData(&data), 
   fPdf(model.GetPdf()),
   fPriorPdf( model.GetPriorPdf()),
   fNuisancePdf(0),
   fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0),
   fPosteriorFunction(0), fApproxPosterior(0),
   fLower(0), fUpper(0), 
   fNLLMin(0),
   fSize(0.05), fLeftSideFraction(0.5), 
   fBrfPrecision(0.00005), 
   fNScanBins(-1),
   fNumIterations(0),
   fValidInterval(false)
{
   // Constructor from a data set and a ModelConfig 
   // model pdf, poi and nuisances will be taken from the ModelConfig
   SetModel(model);
}


BayesianCalculator::~BayesianCalculator()
{
   // destructor
   ClearAll(); 
}

void BayesianCalculator::ClearAll() const { 
   // clear all cached pdf objects
   if (fProductPdf) delete fProductPdf; 
   if (fLogLike) delete fLogLike; 
   if (fLikelihood) delete fLikelihood; 
   if (fIntegratedLikelihood) delete fIntegratedLikelihood; 
   if (fPosteriorPdf) delete fPosteriorPdf;    
   if (fPosteriorFunction) delete fPosteriorFunction; 
   if (fApproxPosterior) delete fApproxPosterior; 
   fPosteriorPdf = 0; 
   fPosteriorFunction = 0; 
   fProductPdf = 0;
   fLogLike = 0; 
   fLikelihood = 0; 
   fIntegratedLikelihood = 0; 
   fLower = 0;
   fUpper = 0;
   fNLLMin = 0; 
   fValidInterval = false;
}

void BayesianCalculator::SetModel(const ModelConfig & model) {
   // set the model to use 
   // The model pdf, prior pdf, parameter of interest and nuisances 
   // will be taken according to the model 

   fPdf = model.GetPdf();
   fPriorPdf =  model.GetPriorPdf(); 
   // assignment operator = does not do a real copy the sets (must use add method) 
   fPOI.removeAll();
   fNuisanceParameters.removeAll();
   fConditionalObs.removeAll();
   if (model.GetParametersOfInterest()) fPOI.add( *(model.GetParametersOfInterest()) );
   if (model.GetNuisanceParameters())  fNuisanceParameters.add( *(model.GetNuisanceParameters() ) );
   if (model.GetConditionalObservables())  fConditionalObs.add( *(model.GetConditionalObservables() ) );
   // remove constant nuisance parameters 
   RemoveConstantParameters(&fNuisanceParameters);

   // invalidate the cached pointers
   ClearAll(); 
}



RooAbsReal* BayesianCalculator::GetPosteriorFunction() const
{
   // Build and return the posterior function (not normalized) as a RooAbsReal
   // the posterior is obtained from the product of the likelihood function and the
   // prior pdf which is then intergated in the nuisance parameters (if existing).
   // A prior function for the nuisance can be specified either in the prior pdf object
   // or in the model itself. If no prior nuisance is specified, but prior parameters are then
   // the integration is performed assuming a flat prior for the nuisance parameters.        
   // NOTE: the return object is managed by the class, users do not need to delete it

   if (fIntegratedLikelihood) return fIntegratedLikelihood; 
   if (fLikelihood) return fLikelihood; 

   // run some sanity checks
   if (!fPdf ) {
      coutE(InputArguments) << "BayesianCalculator::GetPosteriorPdf - missing pdf model" << std::endl;
      return 0;
   }
   if (fPOI.getSize() == 0) {
      coutE(InputArguments) << "BayesianCalculator::GetPosteriorPdf - missing parameter of interest" << std::endl;
      return 0;
   }
   if (fPOI.getSize() > 1) { 
      coutE(InputArguments) << "BayesianCalculator::GetPosteriorPdf - current implementation works only on 1D intervals" << std::endl;
      return 0; 
   }


   RooArgSet* constrainedParams = fPdf->getParameters(*fData);
   // remove the constant parameters
   RemoveConstantParameters(constrainedParams);
   
   //constrainedParams->Print("V");

   // use RooFit::Constrain() to be sure constraints terms are taken into account
   fLogLike = fPdf->createNLL(*fData, RooFit::Constrain(*constrainedParams), RooFit::ConditionalObservables(fConditionalObs) );



   ccoutD(Eval) <<  "BayesianCalculator::GetPosteriorFunction : " 
                << " pdf value " <<  fPdf->getVal() 
                << " neglogLikelihood = " << fLogLike->getVal() << std::endl;
                
   if (fPriorPdf) 
      ccoutD(Eval)  << "\t\t\t priorPOI value " << fPriorPdf->getVal() << std::endl;

   // check that likelihood evaluation is not inifinity 
   double nllVal =  fLogLike->getVal();   
   if ( nllVal > std::numeric_limits<double>::max() ) {
      coutE(Eval) <<  "BayesianCalculator::GetPosteriorFunction : " 
                  << " Negative log likelihood evaluates to infinity " << std::endl 
                  << " Non-const Parameter values : ";
      RooArgList p(*constrainedParams);
      for (int i = 0; i < p.getSize(); ++i) {
         RooRealVar * v = dynamic_cast<RooRealVar *>(&p[i] );
         if (v!=0) ccoutE(Eval) << v->GetName() << " = " << v->getVal() << "   ";
      }
      ccoutE(Eval) << std::endl;
      ccoutE(Eval) << "--  Perform a full likelihood fit of the model before or set more reasanable parameter values"  
                   << std::endl; 
      coutE(Eval) << "BayesianCalculator::GetPosteriorFunction : " << " cannot compute posterior function "  << std::endl; 
      return 0;
   }



   // need do find minimum of log-likelihood in the range to shift function 
   // to avoid numerical errors when we compute the likelihood (overflows in the exponent)
   // N.B.: this works for only 1 parameter of interest otherwise Minuit should be used for finding the minimum
   RooFunctor * nllFunc = fLogLike->functor(fPOI);
   assert(nllFunc);
   ROOT::Math::Functor1D wnllFunc(*nllFunc);
   RooRealVar* poi = dynamic_cast<RooRealVar*>( fPOI.first() ); 
   assert(poi);

   // try to reduce some error messages
   //Bool_t silentMode = (RooMsgService::instance().globalKillBelow() >= RooFit::ERROR || RooMsgService::instance().silentMode()) ;
   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors);


   coutI(Eval) <<  "BayesianCalculator::GetPosteriorFunction : " 
               << " nll value " <<  nllVal << " poi value = " << poi->getVal() << std::endl;


   ROOT::Math::BrentMinimizer1D minim; 
   minim.SetFunction(wnllFunc,poi->getMin(),poi->getMax() );
   bool ret  = minim.Minimize(100,1.E-3,1.E-3);
   fNLLMin = 0; 
   if (ret) fNLLMin = minim.FValMinimum();

   coutI(Eval) << "BayesianCalculator::GetPosteriorFunction : minimum of NLL vs POI for POI =  " 
          << poi->getVal() << " min NLL = " << fNLLMin << std::endl;

   delete nllFunc;

   delete constrainedParams;


   if ( fNuisanceParameters.getSize() == 0 ||  fIntegrationType.Contains("ROOFIT") ) { 

      ccoutD(Eval) << "BayesianCalculator::GetPosteriorFunction : use ROOFIT integration  " 
                   << std::endl;

#ifdef DOLATER // (not clear why this does not work)
      // need to make in this case a likelihood from the nll and make the product with the prior
      TString likeName = TString("likelihood_times_prior_") + TString(fPriorPdf->GetName());   
      TString formula; 
      formula.Form("exp(-@0+%f+log(@1))",fNLLMin);
      fLikelihood = new RooFormulaVar(likeName,formula,RooArgList(*fLogLike,*fPriorPdf));
#else
      // here use RooProdPdf (not very nice) but working

      if (fLogLike) delete fLogLike; 
      if (fProductPdf) { 
         delete fProductPdf;
         fProductPdf = 0; 
      }

      // // create a unique name for the product pdf 
      RooAbsPdf *  pdfAndPrior = fPdf; 
      if (fPriorPdf) { 
         TString prodName = TString("product_") + TString(fPdf->GetName()) + TString("_") + TString(fPriorPdf->GetName() );   
         // save this as data member since it needs to be deleted afterwards
         fProductPdf = new RooProdPdf(prodName,"",RooArgList(*fPdf,*fPriorPdf));
         pdfAndPrior = fProductPdf;
      }

      RooArgSet* constrParams = fPdf->getParameters(*fData);
      // remove the constant parameters
      RemoveConstantParameters(constrParams);
      fLogLike = pdfAndPrior->createNLL(*fData, RooFit::Constrain(*constrParams),RooFit::ConditionalObservables(fConditionalObs) );
      delete constrParams;

      TString likeName = TString("likelihood_times_prior_") + TString(pdfAndPrior->GetName());   
      TString formula; 
      formula.Form("exp(-@0+%f)",fNLLMin);
      fLikelihood = new RooFormulaVar(likeName,formula,RooArgList(*fLogLike));
#endif

      
      // if no nuisance parameter we can just return the likelihood function
      if (fNuisanceParameters.getSize() == 0) { 
         fIntegratedLikelihood = fLikelihood; 
         fLikelihood = 0; 
      }
      else 
         // case of using RooFit for the integration
         fIntegratedLikelihood = fLikelihood->createIntegral(fNuisanceParameters);


   }

   else if ( fIntegrationType.Contains("TOYMC") ) { 
      // compute the posterior as expectation values of the likelihood function 
      // sampling on the nuisance parameters 

      RooArgList nuisParams(fNuisanceParameters); 

      bool doToysEveryIteration = true;
      // if type is 1-TOYMC or TOYMC-1
      if ( fIntegrationType.Contains("1") || fIntegrationType.Contains("ONE")  ) doToysEveryIteration = false;

      RooAbsPdf * samplingPdf = (fNuisancePdf) ? fNuisancePdf : fPdf;
      if (!fNuisancePdf) { 
         ccoutI(Eval) << "BayesianCalculator::GetPosteriorFunction : no nuisance pdf is provided, try using global pdf (this will be slower)"
                      << std::endl;         
      }
      fPosteriorFunction = new PosteriorFunctionFromToyMC(*fLogLike, *samplingPdf, *poi, nuisParams, fPriorPdf, fNLLMin,
                                                          fNumIterations, doToysEveryIteration ); 

      TString name = "toyposteriorfunction_from_"; 
      name += fLogLike->GetName();  
      fIntegratedLikelihood = new RooFunctor1DBinding(name,name,*fPosteriorFunction,*poi);
      
      // need to scan likelihood in this case
      if (fNScanBins <= 0) fNScanBins = 100;
      
   }

   else  { 

      // use ROOT integration method if there are nuisance parameters 

      RooArgList nuisParams(fNuisanceParameters); 
      fPosteriorFunction = new PosteriorFunction(*fLogLike, *poi, nuisParams, fPriorPdf, fIntegrationType, 1., fNLLMin, fNumIterations ); 
      
      TString name = "posteriorfunction_from_"; 
      name += fLogLike->GetName();  
      fIntegratedLikelihood = new RooFunctor1DBinding(name,name,*fPosteriorFunction,*poi);

   }


   if (RooAbsReal::numEvalErrors() > 0) 
      coutW(Eval) << "BayesianCalculator::GetPosteriorFunction : " << RooAbsReal::numEvalErrors() << " errors reported in evaluating log-likelihood function "
                   << std::endl;         
   RooAbsReal::clearEvalErrorLog();
   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors);
   
   return fIntegratedLikelihood;  

}

RooAbsPdf* BayesianCalculator::GetPosteriorPdf() const
{
   // Build and return the posterior pdf (i.e posterior function normalized to all range of poi)
   // Note that an extra integration in the POI is required for the normalization
   // NOTE: user must delete the returned object 
   
   RooAbsReal * plike = GetPosteriorFunction();
   if (!plike) return 0;

   
   // create a unique name on the posterior from the names of the components
   TString posteriorName = this->GetName() + TString("_posteriorPdf_") + plike->GetName(); 

   RooAbsPdf * posteriorPdf = new RooGenericPdf(posteriorName,"@0",*plike);

   return posteriorPdf;
}


RooPlot* BayesianCalculator::GetPosteriorPlot(bool norm, double precision ) const
{
   // return a RooPlot with the posterior  and the credibility region
   // NOTE: User takes ownership of the returned object 

   GetPosteriorFunction(); 

   // if a scan is requested approximate the posterior
   if (fNScanBins > 0) 
      ApproximatePosterior();

   RooAbsReal * posterior = fIntegratedLikelihood; 
   if (norm) { 
      // delete and re-do always posterior pdf (could be invalid after approximating it)
      if (fPosteriorPdf) delete fPosteriorPdf; 
      fPosteriorPdf = GetPosteriorPdf();
      posterior = fPosteriorPdf;
   }
   if (!posterior) return 0;

   if (!fValidInterval) GetInterval();

   RooAbsRealLValue* poi = dynamic_cast<RooAbsRealLValue*>( fPOI.first() );
   assert(poi);


   RooPlot* plot = poi->frame();
   if (!plot) return 0;

   // try to reduce some error messages
   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors);

   plot->SetTitle(TString("Posterior probability of parameter \"")+TString(poi->GetName())+TString("\""));  
   posterior->plotOn(plot,RooFit::Range(fLower,fUpper,kFALSE),RooFit::VLines(),RooFit::DrawOption("F"),RooFit::MoveToBack(),RooFit::FillColor(kGray),RooFit::Precision(precision));
   posterior->plotOn(plot);
   plot->GetYaxis()->SetTitle("posterior function");

   // reset the counts and default mode
   RooAbsReal::clearEvalErrorLog();
   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors);
   
   return plot; 
}

void BayesianCalculator::SetIntegrationType(const char * type) { 
   // set the integration type (possible type are) : 
   // 1D integration ( used when only one nuisance and when the posterior is scanned):  
   //    adaptive , gauss, nonadaptive
   // multidim: 
   //       ADAPTIVE,   adaptive numerical integration
   //                    The parameter numIters (settable with SetNumIters) is  the max number of function calls. 
   //                    It can be reduced to make teh integration faster but it will be difficult to reach the required tolerance 
   //       VEGAS   MC integration method based on importance sampling - numIters is number of function calls 
   //                Extra Vegas parameter can be set using  IntegratorMultiDimOptions class
   //       MISER    MC integration method based on stratified sampling 
   //                See also http://en.wikipedia.org/wiki/Monte_Carlo_integration for VEGAS and MISER description
   //       PLAIN    simple MC integration method, where the max  number of calls can be specified using SetNumIters(numIters)
   //
   // Extra integration types are: 
   //
   //    TOYMC: 
   //       evaluate posterior by generating toy MC for the nuisance parameters. It is a MC 
   //       integration, where the function is sampled according to the nuisance. It is convenient to use when all 
   //       the nuisance are uncorrelated and it is efficient to generate them
   //       The toy are generated by default for each  poi values 
   //       (this method has been proposed and provided by J.P Chou)
   //    1-TOYMC  : same method as before but in this case the toys are generated only one time and then used for
   //       each poi value. It can be convenient when the generation time is much larger than the evaluation time, 
   //       otherwise it is recoomended to re-generate the toy for each poi scanned point of the posterior function
   //
   //
   //    ROOFIT: 
   //       use roofit default integration methods which will produce a nested integral (not reccomended for more
   //       than 1 nuisance parameters) 
   //
   // if type = 0 use default specified via class IntegratorMultiDimOptions::SetDefaultIntegrator
   fIntegrationType = TString(type); 
   fIntegrationType.ToUpper(); 
}



SimpleInterval* BayesianCalculator::GetInterval() const
{
   // Compute the interval. By Default a central interval is computed 
   // and the result is a SimpleInterval object.
   // Using the method (to be called before SetInterval) SetLeftSideTailFraction the user can choose the type of interval.
   // By default the returned interval is a central interval with the confidence level specified 
   // previously in the constructor ( LeftSideTailFraction = 0.5). 
   //  For lower limit use SetLeftSideTailFraction = 1
   //  For upper limit use SetLeftSideTailFraction = 0
   //  for shortest intervals use SetLeftSideTailFraction = -1 or call the method SetShortestInterval()
   // NOTE: The BayesianCalculator covers only the case with one
   // single parameter of interest
   // NOTE: User takes ownership of the returned object

   if (fValidInterval) 
      coutW(Eval) << "BayesianCalculator::GetInterval - recomputing interval for the same CL and same model" << std::endl;

   RooRealVar* poi = dynamic_cast<RooRealVar*>( fPOI.first() );
   if (!poi) { 
      coutE(Eval) << "BayesianCalculator::GetInterval - no parameter of interest is set " << std::endl;
      return 0; 
   } 



   // get integrated likelihood (posterior function) 
   GetPosteriorFunction();

   //Bool_t silentMode = (RooMsgService::instance().globalKillBelow() >= RooFit::ERROR || RooMsgService::instance().silentMode()) ;
   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors);

   if (fLeftSideFraction < 0 ) { 
      // compute short intervals
      ComputeShortestInterval(); 
   }
   else {
      // compute the other intervals

      double lowerCutOff = fLeftSideFraction * fSize; 
      double upperCutOff = 1. - (1.- fLeftSideFraction) * fSize; 


      if (fNScanBins > 0) { 
         ComputeIntervalFromApproxPosterior(lowerCutOff, upperCutOff);
      }
      
      else { 
         // use integration method if there are nuisance parameters 
         if (fNuisanceParameters.getSize() > 0) { 
            ComputeIntervalFromCdf(lowerCutOff, upperCutOff);      
         }
         else { 
            // case of no nuisance - just use createCdf from roofit
            ComputeIntervalUsingRooFit(lowerCutOff, upperCutOff);      
         }
         // case cdf failed (scan then the posterior)
         if (!fValidInterval) { 
            fNScanBins = 100;
            coutW(Eval) << "BayesianCalculator::GetInterval - computing integral from cdf failed - do a scan in "
                        << fNScanBins << " nbins " << std::endl;
            ComputeIntervalFromApproxPosterior(lowerCutOff, upperCutOff);
         }
      }
   }


   // reset the counts and default mode
   if (RooAbsReal::numEvalErrors() > 0) 
      coutW(Eval) << "BayesianCalculator::GetInterval : " << RooAbsReal::numEvalErrors() << " errors reported in evaluating log-likelihood function "
                   << std::endl;         
   
   RooAbsReal::clearEvalErrorLog();
   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors);

   if (!fValidInterval) { 
      fLower = 1; fUpper = 0;
      coutE(Eval) << "BayesianCalculator::GetInterval - cannot compute a valid interval - return a dummy [1,0] interval"
      <<  std::endl;
   }
   else {
      coutI(Eval) << "BayesianCalculator::GetInterval - found a valid interval : [" << fLower << " , " 
                << fUpper << " ]" << std::endl;
   }
   
   TString interval_name = TString("BayesianInterval_a") + TString(this->GetName());
   SimpleInterval * interval = new SimpleInterval(interval_name,*poi,fLower,fUpper,ConfidenceLevel());
   interval->SetTitle("SimpleInterval from BayesianCalculator");
   
   return interval;
}

double BayesianCalculator::GetMode() const { 
   /// Returns the value of the parameter for the point in
   /// parameter-space that is the most likely.
   ///  How do we do if there are points that are equi-probable?
   /// use approximate posterior
   /// t.b.d use real function to find the mode

   ApproximatePosterior(); 
   TH1 * h = fApproxPosterior->GetHistogram();
   return h->GetBinCenter(h->GetMaximumBin() );
   //return  fApproxPosterior->GetMaximumX();
}

void BayesianCalculator::ComputeIntervalUsingRooFit(double lowerCutOff, double upperCutOff ) const { 
   // internal function compute the interval using RooFit

   coutI(Eval) <<  "BayesianCalculator: Compute interval using RooFit:  posteriorPdf + createCdf + RooBrentRootFinder " << std::endl;

   RooRealVar* poi = dynamic_cast<RooRealVar*>( fPOI.first() ); 
   assert(poi);

   fValidInterval = false;
   if (!fPosteriorPdf) fPosteriorPdf = (RooAbsPdf*) GetPosteriorPdf();
   if (!fPosteriorPdf) return;
         
   RooAbsReal* cdf = fPosteriorPdf->createCdf(fPOI,RooFit::ScanNoCdf());
   if (!cdf) return;
         
   RooAbsFunc* cdf_bind = cdf->bindVars(fPOI,&fPOI);
   if (!cdf_bind) return; 

   RooBrentRootFinder brf(*cdf_bind);
   brf.setTol(fBrfPrecision); // set the brf precision
   
   double tmpVal = poi->getVal();  // patch used because findRoot changes the value of poi
   
   bool ret = true; 
   if (lowerCutOff > 0) { 
      double y = lowerCutOff;
      ret &= brf.findRoot(fLower,poi->getMin(),poi->getMax(),y);
   } 
   else 
      fLower = poi->getMin(); 

   if (upperCutOff < 1.0) { 
      double y=upperCutOff;
      ret &= brf.findRoot(fUpper,poi->getMin(),poi->getMax(),y);
   }
   else 
      fUpper = poi->getMax();
   if (!ret)  coutE(Eval) << "BayesianCalculator::GetInterval "
                           << "Error returned from Root finder, estimated interval is not fully correct" 
                           << std::endl;
   else 
      fValidInterval = true; 
   

   poi->setVal(tmpVal); // patch: restore the original value of poi
   
   delete cdf_bind;
   delete cdf;
}

void BayesianCalculator::ComputeIntervalFromCdf(double lowerCutOff, double upperCutOff ) const { 
   // internal function compute the interval using Cdf integration

   fValidInterval = false; 

   coutI(InputArguments) <<  "BayesianCalculator:GetInterval Compute the interval from the posterior cdf " << std::endl;

   RooRealVar* poi = dynamic_cast<RooRealVar*>( fPOI.first() ); 
   assert(poi);
   if (GetPosteriorFunction() == 0) { 
      coutE(InputArguments) <<  "BayesianCalculator::GetInterval() cannot make posterior Function " << std::endl;      
      return;
   }

   // need to remove the constant parameters
   RooArgList bindParams; 
   bindParams.add(fPOI);
   bindParams.add(fNuisanceParameters);
   
   // this code could be put inside the PosteriorCdfFunction
         
   //bindParams.Print("V");
         
   PosteriorCdfFunction cdf(*fLogLike, bindParams, fPriorPdf, fIntegrationType, fNLLMin ); 
   if( cdf.HasError() ) { 
      coutE(Eval) <<  "BayesianCalculator: Numerical error computing CDF integral - try a different method " << std::endl;      
      return;
   }
                  
   //find the roots 

   ROOT::Math::RootFinder rf(kRootFinderType); 
         
   ccoutD(Eval) << "BayesianCalculator::GetInterval - finding roots of posterior using RF " << rf.Name() 
                << " with precision = " << fBrfPrecision;
         
   if (lowerCutOff > 0) { 
      cdf.SetOffset(lowerCutOff);
      ccoutD(NumIntegration) << "Integrating posterior to get cdf and search lower limit at p =" << lowerCutOff << std::endl;
      bool ok = rf.Solve(cdf, poi->getMin(),poi->getMax() , 200,fBrfPrecision, fBrfPrecision); 
      if( cdf.HasError() ) 
         coutW(Eval) <<  "BayesianCalculator: Numerical error integrating the  CDF   " << std::endl;      
      if (!ok) {
         coutE(NumIntegration) << "BayesianCalculator::GetInterval - Error from root finder when searching lower limit !" << std::endl;
         return;
      }
      fLower = rf.Root(); 
   }
   else { 
      fLower = poi->getMin(); 
   }
   if (upperCutOff < 1.0) {  
      cdf.SetOffset(upperCutOff);
      ccoutD(NumIntegration) << "Integrating posterior to get cdf and search upper interval limit at p =" << upperCutOff << std::endl;
      bool ok = rf.Solve(cdf, fLower,poi->getMax() , 200, fBrfPrecision, fBrfPrecision); 
      if( cdf.HasError() ) 
         coutW(Eval) <<  "BayesianCalculator: Numerical error integrating the  CDF   " << std::endl;            
      if (!ok)  {
         coutE(NumIntegration) << "BayesianCalculator::GetInterval - Error from root finder when searching upper limit !" << std::endl;
         return;
      }   
      fUpper = rf.Root(); 
   }
   else { 
      fUpper = poi->getMax(); 
   }

   fValidInterval = true; 
}

void BayesianCalculator::ApproximatePosterior() const { 
   // approximate posterior in nbins using a TF1 
   // scan the poi values and evaluate the posterior at each point 
   // and save the result in a cloned TF1 
   // For each point the posterior is evaluated by integrating the nuisance 
   // parameters 

   if (fApproxPosterior) { 
      // if number of bins of existing function is >= requested one - no need to redo the scan
      if (fApproxPosterior->GetNpx() >= fNScanBins) return;  
      // otherwise redo the scan
      delete fApproxPosterior; 
      fApproxPosterior = 0;
   }      


   RooAbsReal * posterior = GetPosteriorFunction();
   if (!posterior) return; 


   TF1 * tmp = posterior->asTF(fPOI); 
   assert(tmp != 0);
   // binned the function in nbins and evaluate at thos points
   if (fNScanBins > 0)  tmp->SetNpx(fNScanBins);  // if not use default of TF1 (which is 100)

   coutI(Eval) << "BayesianCalculator - scan posterior function in nbins = " << tmp->GetNpx() << std::endl;

   fApproxPosterior = (TF1*) tmp->Clone();
   // save this function for future reuse 
   // I can delete now original posterior and use this approximated copy
   delete tmp;
   TString name = posterior->GetName() + TString("_approx");
   TString title = posterior->GetTitle() + TString("_approx");
   RooAbsReal * posterior2 = new RooTFnBinding(name,title,fApproxPosterior,fPOI);
   if (posterior == fIntegratedLikelihood) { 
      delete fIntegratedLikelihood;
      fIntegratedLikelihood = posterior2; 
   }
   else if (posterior == fLikelihood) { 
      delete fLikelihood; 
      fLikelihood = posterior2;
   }
   else {
      assert(1); // should never happen this case
   }
}

void BayesianCalculator::ComputeIntervalFromApproxPosterior(double lowerCutOff, double upperCutOff ) const { 
   // compute the interval using the approximate posterior function

   ccoutD(Eval) <<  "BayesianCalculator: Compute interval from the approximate posterior " << std::endl;

   ApproximatePosterior();
   if (!fApproxPosterior) return;

   double prob[2]; 
   double limits[2] = {0,0};
   prob[0] = lowerCutOff;
   prob[1] = upperCutOff; 
   fApproxPosterior->GetQuantiles(2,limits,prob);
   fLower = limits[0]; 
   fUpper = limits[1];
   fValidInterval = true; 
}

void BayesianCalculator::ComputeShortestInterval( ) const { 
   // compute the shortest interval
   coutI(Eval) << "BayesianCalculator - computing shortest interval with CL = " << 1.-fSize << std::endl;

   // compute via the approx posterior function
   ApproximatePosterior(); 
   if (!fApproxPosterior) return;

   TH1D * h1 = dynamic_cast<TH1D*>(fApproxPosterior->GetHistogram() );
   assert(h1 != 0);
   h1->SetName(fApproxPosterior->GetName());
   // get bins and sort them 
   double * bins = h1->GetArray(); 
   int n = h1->GetSize()-2; // exclude under/overflow bins
   std::vector<int> index(n);
   TMath::Sort(n, bins, &index[0]); 
   // find cut off as test size
   double sum = 0; 
   double actualCL = 0;
   double upper =  h1->GetXaxis()->GetXmin();
   double lower =  h1->GetXaxis()->GetXmax();
   double norm = h1->GetSumOfWeights();

   for (int i = 0; i < n; ++i)  {
      int idx = index[i]; 
      double p = bins[ idx] / norm;
      sum += p;
      if (sum > 1.-fSize ) { 
         actualCL = sum - p;
         break; 
      }

      if ( h1->GetBinLowEdge(idx) < lower) 
         lower = h1->GetBinLowEdge(idx);
      if ( h1->GetXaxis()->GetBinUpEdge(idx) > upper) 
         upper = h1->GetXaxis()->GetBinUpEdge(idx);
   }

   ccoutD(Eval) << "BayesianCalculator::ComputeShortestInterval - actual interval CL = " 
                << actualCL << " difference from requested is " << (actualCL-(1.-fSize))/fSize*100. << "%  "
                << " limits are [ " << lower << " , " << " upper ] " << std::endl;


   if (lower < upper) { 
      fLower = lower;
      fUpper = upper; 



      if ( std::abs(actualCL-(1.-fSize)) > 0.1*(1.-fSize) )
         coutW(Eval) << "BayesianCalculator::ComputeShortestInterval - actual interval CL = " 
                     << actualCL << " differs more than 10% from desired CL value - must increase nbins " 
                     << n << " to an higher value " << std::endl;
   }
   else 
      coutE(Eval) << "BayesianCalculator::ComputeShortestInterval " << n << " bins are not sufficient " << std::endl;

   fValidInterval = true; 

}




} // end namespace RooStats