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

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 [inline, virtual]

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;
}