#include <AsymptoticNew.h>

Inheritance diagram for AsymptoticNew:

Public Member Functions
virtual void	applyDefaultOptions ()
virtual void	applyOptions (const boost::program_options::variables_map &vm)
	AsymptoticNew ()
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)
std::vector< std::pair< float, float > >	runLimit (RooWorkspace w, RooStats::ModelConfig mc_s, RooStats::ModelConfig mc_b, RooAbsData &data, double &limit, double &limitErr, const double hint)
Private Attributes
static double	minrscan_ = 0.
Static Private Attributes
static double	maxrscan_ = 20.
static int	nscanpoints_ = 20
static bool	qtilde_ = true
static double	rValue_
static std::string	what_ = "both"

Detailed Description

new CLs asymptotic limits

Author:: Nicholas Wardle (mostly taken from Asymptotic class)

Definition at line 16 of file AsymptoticNew.h.

Constructor & Destructor Documentation

AsymptoticNew::AsymptoticNew ( )

Definition at line 21 of file AsymptoticNew.cc.

References maxrscan_, minrscan_, nscanpoints_, LimitAlgo::options_, qtilde_, and what_.

                             :
LimitAlgo("AsymptoticNew specific options"){
    options_.add_options()
        ("run", boost::program_options::value<std::string>(&what_)->default_value(what_), "What to run: both (default), observed, expected.")
        ("nPoints", boost::program_options::value<int>(&nscanpoints_)->default_value(nscanpoints_), "Number of points in scan for CLs")
        ("qtilde", boost::program_options::value<bool>(&qtilde_)->default_value(qtilde_),  "Allow only non-negative signal strengths (default is true).")
        ("scanMin", boost::program_options::value<double>(&minrscan_)->default_value(minrscan_),  "Minimum value for scan in r.")
        ("scanMax", boost::program_options::value<double>(&maxrscan_)->default_value(maxrscan_),  "Maximum value for scan in r.")
  ;
}

Member Function Documentation

void AsymptoticNew::applyDefaultOptions ( ) [virtual]

Reimplemented from LimitAlgo.

Definition at line 38 of file AsymptoticNew.cc.

References maxrscan_, minrscan_, nscanpoints_, qtilde_, and what_.

                                        { 
    what_ = "observed";
    nscanpoints_ = 20;
    qtilde_ = true;
    minrscan_=0.;
    maxrscan_=20.;
}

void AsymptoticNew::applyOptions ( const boost::program_options::variables_map & vm ) [virtual]

Reimplemented from LimitAlgo.

Definition at line 32 of file AsymptoticNew.cc.

References what_.

                                                                            {
        if (what_ != "observed" && what_ != "expected" && what_ != "both") 
            throw std::invalid_argument("Asymptotic: option 'run' can only be 'observed', 'expected' or 'both' (the default)");

}

virtual const std::string& AsymptoticNew::name ( void ) const [inline, virtual]

Implements LimitAlgo.

Definition at line 23 of file AsymptoticNew.h.

{ static std::string name_ = "AsymptoticNew"; return name_; }

bool AsymptoticNew::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 46 of file AsymptoticNew.cc.

References gather_cfg::cout, DEBUG, runLimit(), dqm::qstatus::WARNING, and what_.

                                                                                                                                                                    {
    RooFitGlobalKillSentry silence(verbose <= 1 ? RooFit::WARNING : RooFit::DEBUG);
    if (verbose > 0) std::cout << "Will compute " << what_ << " limit(s) " << std::endl;
  
    std::vector<std::pair<float,float> > expected;
    expected = runLimit(w, mc_s, mc_b, data, limit, limitErr, hint);

    if (verbose >= 0) {
        const char *rname = mc_s->GetParametersOfInterest()->first()->GetName();
        std::cout << "\n -- AsymptoticNew -- " << "\n";
        if (what_ != "expected") {
            printf("Observed Limit: %s < %6.4f\n", rname, limit);
        }
        for (std::vector<std::pair<float,float> >::const_iterator it = expected.begin(), ed = expected.end(); it != ed; ++it) {
            printf("Expected %4.1f%%: %s < %6.4f\n", it->first*100, rname, it->second);
        }
        std::cout << std::endl;
    }

    return true;
}

std::vector< std::pair< float, float > > AsymptoticNew::runLimit	(	RooWorkspace *	w,
		RooStats::ModelConfig *	mc_s,
		RooStats::ModelConfig *	mc_b,
		RooAbsData &	data,
		double &	limit,
		double &	limitErr,
		const double *	hint
	)

Definition at line 68 of file AsymptoticNew.cc.

References Combine::commitPoint(), maxrscan_, minrscan_, nscanpoints_, qtilde_, alignCSCRings::r, and what_.

Referenced by run().

                                                                                                                                                                                                     {
 
  RooRealVar *poi = (RooRealVar*)mc_s->GetParametersOfInterest()->first();
  mc_s->SetSnapshot(*poi); 
  double oldval = poi->getVal(); 
  poi->setVal(0); mc_b->SetSnapshot(*poi); 
  poi->setVal(oldval);

  // Set up asymptotic calculator
  AsymptoticCalculator * ac = new RooStats::AsymptoticCalculator(data, *mc_b, *mc_s);
  AsymptoticCalculator::SetPrintLevel(verbose>2);
  ac->SetOneSided(true); 
  ac->SetQTilde(qtilde_);
  HypoTestInverter calc(*ac);
  calc.SetVerbose(verbose>2);
  calc.UseCLs(true);
  
  calc.SetFixedScan(nscanpoints_,minrscan_,maxrscan_); 
  RooStats::HypoTestInverterResult * r = calc.GetInterval();

  // Expected Median + bands
  // bands will be based on Z-values
  std::vector<std::pair<float,float> > expected;
  const double quantiles[5] = { 0.025, 0.16, 0.50, 0.84, 0.975 };
  const double zvals[5]     = { -2, -1, 0, 1, 2 };

  for (int iq = 0; iq < 5; ++iq) {
        limit = r->GetExpectedUpperLimit(zvals[iq]);
        limitErr = 0;
        Combine::commitPoint(true, quantiles[iq]);
        expected.push_back(std::pair<float,float>(quantiles[iq], limit));
  }

  
  // Observed Limit
  if (what_!="expected")   {
    limit = r->UpperLimit();
  } else {
    limit = 0 ;
  }

  return expected;
   
}