BayesianToyMC Class Reference

#include <BayesianToyMC.h>

Inheritance diagram for BayesianToyMC:

Public Member Functions
virtual void	applyOptions (const boost::program_options::variables_map &vm)
	BayesianToyMC ()
virtual const std::string &	name () const
virtual bool	run (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint)
Public Member Functions inherited from LimitAlgo
virtual void	applyDefaultOptions ()
	LimitAlgo ()
	LimitAlgo (const char *desc)
const boost::program_options::options_description &	options () const
virtual void	setNToys (const int)
virtual void	setToyNumber (const int)
virtual	~LimitAlgo ()

Static Private Attributes
static float	hintSafetyFactor_ = 5.
	Safety factor for hint (integrate up to this number of times the hinted limit) More...
static std::string	integrationType_ = "toymc"
	numerical integration algorithm More...
static int	numIters_ = 1000
	number of iterations for each toy mc computation More...
static unsigned int	tries_ = 1
	number of toy mc computations to run More...

Additional Inherited Members
Protected Attributes inherited from LimitAlgo
boost::program_options::options_description	options_

Detailed Description

Interface to BayesianCalculator used with toymc algorithms

Author

Giovanni Petrucciani (UCSD)

Definition at line 13 of file BayesianToyMC.h.

Constructor & Destructor Documentation

BayesianToyMC::BayesianToyMC

(

)

Definition at line 21 of file BayesianToyMC.cc.

References hintSafetyFactor_, integrationType_, numIters_, LimitAlgo::options_, and tries_.

                             :
    LimitAlgo("BayesianToyMC specific options")
{
    options_.add_options()
        ("integrationType", boost::program_options::value<std::string>(&integrationType_)->default_value(integrationType_), "Integration algorithm to use")
        ("tries",      boost::program_options::value<unsigned int>(&tries_)->default_value(tries_), "Number of times to run the ToyMC on the same data")
        ("numIters,i", boost::program_options::value<int>(&numIters_)->default_value(numIters_),    "Number of iterations or calls used within iteration (0=ROOT Default)")
        ("hintSafetyFactor",
                boost::program_options::value<float>(&hintSafetyFactor_)->default_value(hintSafetyFactor_),
                "set range of integration equal to this number of times the hinted limit")
        ;
}

Member Function Documentation

void BayesianToyMC::applyOptions ( const boost::program_options::variables_map &  vm )

virtual

Reimplemented from LimitAlgo.

Definition at line 34 of file BayesianToyMC.cc.

References gather_cfg::cout, integrationType_, numIters_, tries_, and withSystematics.

                                                                            {
    if (!withSystematics) {
        std::cout << "BayesianToyMC: when running with no systematics, BayesianToyMC is identical to BayesianSimple." << std::endl;
        tries_    = 1;
        numIters_ = 0;
        integrationType_ = "";
    }
}

virtual const std::string& BayesianToyMC::name ( void  ) const

inlinevirtual

Implements LimitAlgo.

Definition at line 18 of file BayesianToyMC.h.

                                         {
    static const std::string name("BayesianToyMC");
    return name;
  }

bool BayesianToyMC::run ( RooWorkspace *  w, RooStats::ModelConfig *  mc_s, RooStats::ModelConfig *  mc_b, RooAbsData &  data, double &  limit, double &  limitErr, const double *  hint  )

virtual

Implements LimitAlgo.

Definition at line 42 of file BayesianToyMC.cc.

References cl, gather_cfg::cout, hintSafetyFactor_, i, integrationType_, utils::makeNuisancePdf(), numIters_, alignCSCRings::r, mathSSE::sqrt(), tries_, and withSystematics.

                                                                                                                                                                    {
  RooArgSet  poi(*mc_s->GetParametersOfInterest());

  RooRealVar *r = dynamic_cast<RooRealVar *>(poi.first());
  if ((hint != 0) && (*hint > r->getMin())) {
    r->setMax(hintSafetyFactor_*(*hint));
  }
  double rMax = r->getMax();
  
  std::auto_ptr<RooStats::ModelConfig> mc_noNuis(0);
  if (!withSystematics && mc_s->GetNuisanceParameters() != 0) {
    mc_noNuis.reset(new RooStats::ModelConfig(w));
    mc_noNuis->SetPdf(*mc_s->GetPdf());
    mc_noNuis->SetObservables(*mc_s->GetObservables());
    mc_noNuis->SetParametersOfInterest(*mc_s->GetParametersOfInterest());
    mc_noNuis->SetPriorPdf(*mc_s->GetPriorPdf());
    mc_s = mc_noNuis.get();
  }
  std::auto_ptr<RooAbsPdf> nuisancePdf;
  for (;;) {
    limit = 0; limitErr = 0; bool rerun = false;
    for (unsigned int i = 0; i < tries_; ++i) {
        BayesianCalculator bcalc(data, *mc_s);
        bcalc.SetLeftSideTailFraction(0);
        bcalc.SetConfidenceLevel(cl); 
        if (!integrationType_.empty()) bcalc.SetIntegrationType(integrationType_.c_str());
        if (numIters_) bcalc.SetNumIters(numIters_);
        if (integrationType_ == "toymc") {
            if (nuisancePdf.get() == 0)  nuisancePdf.reset(utils::makeNuisancePdf(*mc_s));
            bcalc.ForceNuisancePdf(*nuisancePdf);
        }
        // get the interval
        std::auto_ptr<SimpleInterval> bcInterval(bcalc.GetInterval());
        if (bcInterval.get() == 0) return false;
        double lim = bcInterval->UpperLimit();
        // check against bound
        if (lim >= 0.5*r->getMax()) { 
            std::cout << "Limit " << r->GetName() << " < " << lim << "; " << r->GetName() << " max < " << r->getMax() << std::endl;
            if (r->getMax()/rMax > 20) return false;
            r->setMax(r->getMax()*2); 
            rerun = true; break;
        }
        // add to running sum(x) and sum(x2)
        limit    += lim;
        limitErr += lim*lim;
        if (tries_ > 1 && verbose > 1) std::cout << " - limit from try " << i << ": " << lim << std::endl;
    }
    if (rerun) continue;
    limit /= tries_; 
    limitErr = (tries_ > 1 ? std::sqrt((limitErr/tries_ - limit*limit)/((tries_-1)*tries_)) : 0);
    if (verbose > -1) {
        std::cout << "\n -- BayesianToyMC -- " << "\n";
        if (limitErr > 0) {
            std::cout << "Limit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " @ " << cl * 100 << "% CL" << std::endl;
        } else {
            std::cout << "Limit: " << r->GetName() << " < " << limit << " @ " << cl * 100 << "% CL" << std::endl;
        }
    }
    break;
  }
  return true;
}

Member Data Documentation

float BayesianToyMC::hintSafetyFactor_ = 5.

staticprivate

Safety factor for hint (integrate up to this number of times the hinted limit)

Definition at line 30 of file BayesianToyMC.h.

Referenced by BayesianToyMC(), and run().

std::string BayesianToyMC::integrationType_ = "toymc"

staticprivate

numerical integration algorithm

Definition at line 24 of file BayesianToyMC.h.

Referenced by applyOptions(), BayesianToyMC(), and run().

int BayesianToyMC::numIters_ = 1000

staticprivate

number of iterations for each toy mc computation

Definition at line 26 of file BayesianToyMC.h.

Referenced by applyOptions(), BayesianToyMC(), and run().

unsigned int BayesianToyMC::tries_ = 1

staticprivate

number of toy mc computations to run

Definition at line 28 of file BayesianToyMC.h.

Referenced by applyOptions(), BayesianToyMC(), and run().