|
|
|
#include <HybridNew.h>
Classes | |
| struct | Setup |
Public Types | |
| enum | WorkingMode { MakeLimit, MakeSignificance, MakePValues, MakeTestStatistics, MakeSignificanceTestStatistics } |
Public Member Functions | |
| virtual void | applyDefaultOptions () |
| virtual void | applyOptions (const boost::program_options::variables_map &vm) |
| HybridNew () | |
| 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) |
| virtual bool | runLimit (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) |
| virtual bool | runSignificance (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) |
| virtual bool | runSinglePoint (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) |
| virtual bool | runTestStatistics (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) |
Private Member Functions | |
| void | applyClsQuantile (RooStats::HypoTestResult &hcres) |
| void | applyExpectedQuantile (RooStats::HypoTestResult &hcres) |
| void | applySignalQuantile (RooStats::HypoTestResult &hcres) |
| void | clearGrid () |
| std::auto_ptr < RooStats::HybridCalculator > | create (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, const RooAbsCollection &rVals, Setup &setup) |
| generic version | |
| std::auto_ptr < RooStats::HybridCalculator > | create (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double rVal, Setup &setup) |
| specific version used in unidimensional searches (calls the generic one) | |
| std::pair< double, double > | eval (const RooStats::HypoTestResult &hcres, const RooAbsCollection &rVals) |
| std::pair< double, double > | eval (const RooStats::HypoTestResult &hcres, double rVal) |
| std::pair< double, double > | eval (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double rVal, bool adaptive=false, double clsTarget=-1) |
| specific version used in unidimensional searches (calls the generic one) | |
| std::pair< double, double > | eval (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, const RooAbsCollection &rVals, bool adaptive=false, double clsTarget=-1) |
| generic version | |
| std::pair< double, double > | eval (RooStats::HybridCalculator &hc, const RooAbsCollection &rVals, bool adaptive=false, double clsTarget=-1) |
| RooStats::HypoTestResult * | evalGeneric (RooStats::HybridCalculator &hc, bool forceNoFork=false) |
| RooStats::HypoTestResult * | evalWithFork (RooStats::HybridCalculator &hc) |
| void | readGrid (TDirectory *directory, double rMin, double rMax) |
| RooStats::HypoTestResult * | readToysFromFile (const RooAbsCollection &rVals) |
| void | setupPOI (RooStats::ModelConfig *mc_s) |
| void | updateGridData (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, bool smart, double clsTarget) |
| void | updateGridDataFC (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, bool smart, double clsTarget) |
| std::pair< double, double > | updateGridPoint (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, std::map< double, RooStats::HypoTestResult * >::iterator point) |
| void | useGrid () |
| void | validateOptions () |
Private Attributes | |
| static bool | clsQuantiles_ = true |
| static bool | fitNuisances_ = false |
| static unsigned int | fork_ = 1 |
| static bool | genGlobalObs_ = false |
| std::map< double, RooStats::HypoTestResult * > | grid_ |
| static bool | importanceSamplingAlt_ = false |
| static double | interpAccuracy_ = 0.2 |
| std::auto_ptr< TGraphErrors > | limitPlot_ |
| float | mass_ |
| unsigned int | perf_totalToysRun_ |
| static double | rAbsAccuracy_ = 0.1 |
| static bool | readHybridResults_ = false |
| static double | rRelAccuracy_ = 0.05 |
| static std::string | testStat_ = "LEP" |
Static Private Attributes | |
| static std::string | algo_ = "logSecant" |
| static bool | CLs_ = false |
| static double | clsAccuracy_ = 0.005 |
| static bool | expectedFromGrid_ = false |
| static bool | fullBToys_ = false |
| static bool | fullGrid_ = false |
| static bool | genNuisances_ = true |
| static std::string | gridFile_ = "" |
| static bool | importanceSamplingNull_ = false |
| static unsigned int | iterations_ = 1 |
| static std::string | minimizerAlgo_ = "Minuit2" |
| static float | minimizerTolerance_ = 1e-2 |
| static unsigned int | nCpu_ = 0 |
| static bool | newToyMCSampler_ = true |
| static bool | noUpdateGrid_ = false |
| static unsigned int | nToys_ = 500 |
| static bool | optimizeProductPdf_ = true |
| static bool | optimizeTestStatistics_ = true |
| static std::string | plot_ |
| static float | quantileForExpectedFromGrid_ = 0.5 |
| static bool | reportPVal_ = false |
| static bool | rMaxSet_ = false |
| static bool | rMinSet_ = false |
| static std::string | rule_ = "CLs" |
| static std::string | rValue_ = "1.0" |
| static RooArgSet | rValues_ |
| static bool | saveGrid_ = false |
| static bool | saveHybridResult_ = false |
| static WorkingMode | workingMode_ = MakeLimit |
Module to compute limits by tossing toys (CLs, CLs+b, Feldman-Cousins), and related significances
Definition at line 21 of file HybridNew.h.
| MakeLimit | |
| MakeSignificance | |
| MakePValues | |
| MakeTestStatistics | |
| MakeSignificanceTestStatistics |
Definition at line 38 of file HybridNew.h.
| HybridNew::HybridNew | ( | ) |
Definition at line 89 of file HybridNew.cc.
References algo_, clsAccuracy_, clsQuantiles_, fitNuisances_, fork_, genGlobalObs_, genNuisances_, gridFile_, interpAccuracy_, iterations_, minimizerAlgo_, minimizerTolerance_, nCpu_, newToyMCSampler_, nToys_, optimizeProductPdf_, optimizeTestStatistics_, LimitAlgo::options_, plot_, quantileForExpectedFromGrid_, rAbsAccuracy_, rRelAccuracy_, rule_, rValue_, and testStat_.
: LimitAlgo("HybridNew specific options") { options_.add_options() ("rule", boost::program_options::value<std::string>(&rule_)->default_value(rule_), "Rule to use: CLs, CLsplusb") ("testStat",boost::program_options::value<std::string>(&testStat_)->default_value(testStat_), "Test statistics: LEP, TEV, LHC (previously known as Atlas), Profile.") ("singlePoint", boost::program_options::value<std::string>(&rValue_)->default_value(rValue_), "Just compute CLs for the given value of the parameter of interest. In case of multiple parameters, use a syntax 'name=value,name2=value2,...'") ("onlyTestStat", "Just compute test statistics, not actual p-values (works only with --singlePoint)") ("generateNuisances", boost::program_options::value<bool>(&genNuisances_)->default_value(genNuisances_), "Generate nuisance parameters for each toy") ("generateExternalMeasurements", boost::program_options::value<bool>(&genGlobalObs_)->default_value(genGlobalObs_), "Generate external measurements for each toy, taken from the GlobalObservables of the ModelConfig") ("fitNuisances", boost::program_options::value<bool>(&fitNuisances_)->default_value(fitNuisances_), "Fit the nuisances, when not generating them.") ("searchAlgo", boost::program_options::value<std::string>(&algo_)->default_value(algo_), "Algorithm to use to search for the limit (bisection, logSecant)") ("toysH,T", boost::program_options::value<unsigned int>(&nToys_)->default_value(nToys_), "Number of Toy MC extractions to compute CLs+b, CLb and CLs") ("clsAcc", boost::program_options::value<double>(&clsAccuracy_ )->default_value(clsAccuracy_), "Absolute accuracy on CLs to reach to terminate the scan") ("rAbsAcc", boost::program_options::value<double>(&rAbsAccuracy_)->default_value(rAbsAccuracy_), "Absolute accuracy on r to reach to terminate the scan") ("rRelAcc", boost::program_options::value<double>(&rRelAccuracy_)->default_value(rRelAccuracy_), "Relative accuracy on r to reach to terminate the scan") ("interpAcc", boost::program_options::value<double>(&interpAccuracy_)->default_value(interpAccuracy_), "Minimum uncertainty from interpolation delta(x)/(max(x)-min(x))") ("iterations,i", boost::program_options::value<unsigned int>(&iterations_)->default_value(iterations_), "Number of times to throw 'toysH' toys to compute the p-values (for --singlePoint if clsAcc is set to zero disabling adaptive generation)") ("fork", boost::program_options::value<unsigned int>(&fork_)->default_value(fork_), "Fork to N processes before running the toys (set to 0 for debugging)") ("nCPU", boost::program_options::value<unsigned int>(&nCpu_)->default_value(nCpu_), "Use N CPUs with PROOF Lite (experimental!)") ("saveHybridResult", "Save result in the output file") ("readHybridResults", "Read and merge results from file (requires 'toysFile' or 'grid')") ("grid", boost::program_options::value<std::string>(&gridFile_), "Use the specified file containing a grid of SamplingDistributions for the limit (implies readHybridResults).\n For --singlePoint or --signif use --toysFile=x.root --readHybridResult instead of this.") ("expectedFromGrid", boost::program_options::value<float>(&quantileForExpectedFromGrid_)->default_value(0.5), "Use the grid to compute the expected limit for this quantile") ("signalForSignificance", boost::program_options::value<std::string>()->default_value("1"), "Use this value of the parameter of interest when generating signal toys for expected significance (same syntax as --singlePoint)") ("clsQuantiles", boost::program_options::value<bool>(&clsQuantiles_)->default_value(clsQuantiles_), "Compute correct quantiles of CLs or CLsplusb instead of assuming they're the same as CLb ones") //("importanceSamplingNull", boost::program_options::value<bool>(&importanceSamplingNull_)->default_value(importanceSamplingNull_), // "Enable importance sampling for null hypothesis (background only)") //("importanceSamplingAlt", boost::program_options::value<bool>(&importanceSamplingAlt_)->default_value(importanceSamplingAlt_), // "Enable importance sampling for alternative hypothesis (signal plus background)") ("optimizeTestStatistics", boost::program_options::value<bool>(&optimizeTestStatistics_)->default_value(optimizeTestStatistics_), "Use optimized test statistics if the likelihood is not extended (works for LEP and TEV test statistics).") ("optimizeProductPdf", boost::program_options::value<bool>(&optimizeProductPdf_)->default_value(optimizeProductPdf_), "Optimize the code factorizing pdfs") ("minimizerAlgo", boost::program_options::value<std::string>(&minimizerAlgo_)->default_value(minimizerAlgo_), "Choice of minimizer used for profiling (Minuit vs Minuit2)") ("minimizerTolerance", boost::program_options::value<float>(&minimizerTolerance_)->default_value(minimizerTolerance_), "Tolerance for minimizer used for profiling") ("plot", boost::program_options::value<std::string>(&plot_), "Save a plot of the result (test statistics distributions or limit scan)") ("frequentist", "Shortcut to switch to Frequentist mode (--generateNuisances=0 --generateExternalMeasurements=1 --fitNuisances=1)") ("newToyMCSampler", boost::program_options::value<bool>(&newToyMCSampler_)->default_value(newToyMCSampler_), "Use new ToyMC sampler with support for mixed binned-unbinned generation. On by default, you can turn it off if it doesn't work for your workspace.") ("fullGrid", "Evaluate p-values at all grid points, without optimitations") ("saveGrid", "Save CLs or (or FC p-value) at all grid points in the output tree. The value of 'r' is saved in the 'limit' branch, while the CLs or p-value in the 'quantileExpected' branch and the uncertainty on 'limitErr' (since there's no quantileExpectedErr)") ("noUpdateGrid", "Do not update test statistics at grid points") ("fullBToys", "Run as many B toys as S ones (default is to run 1/4 of b-only toys)") ("pvalue", "Report p-value instead of significance (when running with --significance)") ; }
| void HybridNew::applyClsQuantile | ( | RooStats::HypoTestResult & | hcres | ) | [private] |
New test implementation, scales as N*log(N)
Definition at line 1084 of file HybridNew.cc.
References CLs_, gather_cfg::cout, GOODCOLL_filter_cfg::cut, i, getHLTprescales::index, match(), n, quantileForExpectedFromGrid_, edm::second(), python::multivaluedict::sort(), and makeHLTPrescaleTable::values.
Referenced by applyExpectedQuantile().
{
RooStats::SamplingDistribution * bDistribution = hcres.GetNullDistribution(), * sDistribution = hcres.GetAltDistribution();
const std::vector<Double_t> & bdist = bDistribution->GetSamplingDistribution();
const std::vector<Double_t> & bweight = bDistribution->GetSampleWeights();
const std::vector<Double_t> & sdist = sDistribution->GetSamplingDistribution();
const std::vector<Double_t> & sweight = sDistribution->GetSampleWeights();
TStopwatch timer;
timer.Start();
std::vector<std::pair<double,double> > bcumul; bcumul.reserve(bdist.size());
std::vector<std::pair<double,double> > scumul; scumul.reserve(sdist.size());
double btot = 0, stot = 0;
for (std::vector<Double_t>::const_iterator it = bdist.begin(), ed = bdist.end(), itw = bweight.begin(); it != ed; ++it, ++itw) {
bcumul.push_back(std::pair<double,double>(*it, *itw));
btot += *itw;
}
for (std::vector<Double_t>::const_iterator it = sdist.begin(), ed = sdist.end(), itw = sweight.begin(); it != ed; ++it, ++itw) {
scumul.push_back(std::pair<double,double>(*it, *itw));
stot += *itw;
}
double sinv = 1.0/stot, binv = 1.0/btot, runningSum;
// now compute integral distribution of Q(s+b data) so that we can quickly compute the CL_{s+b} for all test stats.
std::sort(scumul.begin(), scumul.end());
runningSum = 0;
for (std::vector<std::pair<double,double> >::reverse_iterator it = scumul.rbegin(), ed = scumul.rend(); it != ed; ++it) {
runningSum += it->second;
it->second = runningSum * sinv;
}
std::sort(bcumul.begin(), bcumul.end());
std::vector<std::pair<double,std::pair<double,double> > > xcumul; xcumul.reserve(bdist.size());
runningSum = 0;
std::vector<std::pair<double,double> >::const_iterator sbegin = scumul.begin(), send = scumul.end();
//int k = 0;
for (std::vector<std::pair<double,double> >::const_reverse_iterator it = bcumul.rbegin(), ed = bcumul.rend(); it != ed; ++it) {
runningSum += it->second;
std::vector<std::pair<double,double> >::const_iterator match = std::upper_bound(sbegin, send, std::pair<double,double>(it->first, 0));
if (match == send) {
//std::cout << "Did not find match, for it->first == " << it->first << ", as back = ( " << scumul.back().first << " , " << scumul.back().second << " ) " << std::endl;
double clsb = (scumul.back().second > 0.5 ? 1.0 : 0.0), clb = runningSum*binv, cls = clsb / clb;
xcumul.push_back(std::make_pair(CLs_ ? cls : clsb, *it));
} else {
double clsb = match->second, clb = runningSum*binv, cls = clsb / clb;
//if ((++k) % 100 == 0) printf("At %+8.5f CLb = %6.4f, CLsplusb = %6.4f, CLs =%7.4f\n", it->first, clb, clsb, cls);
xcumul.push_back(std::make_pair(CLs_ ? cls : clsb, *it));
}
}
// sort
std::sort(xcumul.begin(), xcumul.end());
// get quantile
runningSum = 0; double cut = quantileForExpectedFromGrid_ * btot;
for (std::vector<std::pair<double,std::pair<double,double> > >::const_iterator it = xcumul.begin(), ed = xcumul.end(); it != ed; ++it) {
runningSum += it->second.second;
if (runningSum >= cut) {
hcres.SetTestStatisticData(it->second.first);
//std::cout << "CLs quantile = " << it->first << " for test stat = " << it->second.first << std::endl;
break;
}
}
//std::cout << "CLs quantile = " << (CLs_ ? hcres.CLs() : hcres.CLsplusb()) << std::endl;
//std::cout << "Computed quantiles in " << timer.RealTime() << " s" << std::endl;
#if 0
timer.Start();
std::vector<std::pair<double, double> > values(bdist.size());
for (int i = 0, n = bdist.size(); i < n; ++i) {
hcres.SetTestStatisticData( bdist[i] );
values[i] = std::pair<double, double>(CLs_ ? hcres.CLs() : hcres.CLsplusb(), bdist[i]);
}
std::sort(values.begin(), values.end());
int index = std::min<int>(floor((1.-quantileForExpectedFromGrid_) * values.size()+0.5), values.size());
std::cout << "CLs quantile = " << values[index].first << " for test stat = " << values[index].second << std::endl;
hcres.SetTestStatisticData(values[index].second);
std::cout << "CLs quantile = " << (CLs_ ? hcres.CLs() : hcres.CLsplusb()) << " for test stat = " << values[index].second << std::endl;
std::cout << "Computed quantiles in " << timer.RealTime() << " s" << std::endl;
#endif
}
| void HybridNew::applyDefaultOptions | ( | ) | [virtual] |
Reimplemented from LimitAlgo.
Definition at line 180 of file HybridNew.cc.
References MakeLimit, validateOptions(), and workingMode_.
{
workingMode_ = MakeLimit;
validateOptions();
}
| void HybridNew::applyExpectedQuantile | ( | RooStats::HypoTestResult & | hcres | ) | [private] |
Definition at line 1068 of file HybridNew.cc.
References applyClsQuantile(), applySignalQuantile(), clsQuantiles_, gather_cfg::cout, expectedFromGrid_, MakeSignificance, MakeSignificanceTestStatistics, quantileForExpectedFromGrid_, python::multivaluedict::sort(), and workingMode_.
Referenced by eval(), runSignificance(), runTestStatistics(), and updateGridPoint().
{
if (expectedFromGrid_) {
if (workingMode_ == MakeSignificance || workingMode_ == MakeSignificanceTestStatistics) {
applySignalQuantile(hcres);
} else if (clsQuantiles_) {
applyClsQuantile(hcres);
} else {
std::vector<Double_t> btoys = hcres.GetNullDistribution()->GetSamplingDistribution();
std::sort(btoys.begin(), btoys.end());
Double_t testStat = btoys[std::min<int>(floor((1.-quantileForExpectedFromGrid_) * btoys.size()+0.5), btoys.size())];
if (verbose > 0) std::cout << "Text statistics for " << quantileForExpectedFromGrid_ << " quantile: " << testStat << std::endl;
hcres.SetTestStatisticData(testStat);
//std::cout << "CLs quantile = " << (CLs_ ? hcres.CLs() : hcres.CLsplusb()) << " for test stat = " << testStat << std::endl;
}
}
}
| void HybridNew::applyOptions | ( | const boost::program_options::variables_map & | vm | ) | [virtual] |
Reimplemented from LimitAlgo.
Definition at line 135 of file HybridNew.cc.
References dtNoiseDBValidation_cfg::cerr, gather_cfg::cout, doSignificance_, expectedFromGrid_, fitNuisances_, fullBToys_, fullGrid_, g_quantileExpected_, genGlobalObs_, genNuisances_, MakeLimit, MakePValues, MakeSignificance, MakeSignificanceTestStatistics, MakeTestStatistics, mass_, noUpdateGrid_, quantileForExpectedFromGrid_, readHybridResults_, reportPVal_, rMaxSet_, rMinSet_, rValue_, saveGrid_, saveHybridResult_, testStat_, validateOptions(), withSystematics, and workingMode_.
{
rMinSet_ = vm.count("rMin")>0; rMaxSet_ = vm.count("rMax")>0;
if (vm.count("expectedFromGrid") && !vm["expectedFromGrid"].defaulted()) {
//if (!vm.count("grid")) throw std::invalid_argument("HybridNew: Can't use --expectedFromGrid without --grid!");
if (quantileForExpectedFromGrid_ <= 0 || quantileForExpectedFromGrid_ >= 1.0) throw std::invalid_argument("HybridNew: the quantile for the expected limit must be between 0 and 1");
expectedFromGrid_ = true;
g_quantileExpected_ = quantileForExpectedFromGrid_;
}
if (vm.count("frequentist")) {
genNuisances_ = 0; genGlobalObs_ = withSystematics; fitNuisances_ = withSystematics;
if (vm["testStat"].defaulted()) testStat_ = "LHC";
if (vm["toys"].as<int>() > 0 and vm.count("toysFrequentist")) {
if (vm["fitNuisances"].defaulted() && withSystematics) {
std::cout << "When tossing frequenst toys outside the HybridNew, the nuisances will not be refitted for each toy by default. This can be changed by specifying explicitly the fitNuisances option" << std::endl;
fitNuisances_ = false;
}
}
}
if (genGlobalObs_ && genNuisances_) {
std::cerr << "ALERT: generating both global observables and nuisance parameters at the same time is not validated." << std::endl;
}
if (!vm["singlePoint"].defaulted()) {
if (doSignificance_) throw std::invalid_argument("HybridNew: Can't use --significance and --singlePoint at the same time");
workingMode_ = ( vm.count("onlyTestStat") ? MakeTestStatistics : MakePValues );
} else if (vm.count("onlyTestStat")) {
if (doSignificance_) workingMode_ = MakeSignificanceTestStatistics;
else throw std::invalid_argument("HybridNew: --onlyTestStat works only with --singlePoint or --significance");
} else if (doSignificance_) {
workingMode_ = MakeSignificance;
rValue_ = vm["signalForSignificance"].as<std::string>();
} else {
workingMode_ = MakeLimit;
}
saveHybridResult_ = vm.count("saveHybridResult");
readHybridResults_ = vm.count("readHybridResults") || vm.count("grid");
if (readHybridResults_ && !(vm.count("toysFile") || vm.count("grid"))) throw std::invalid_argument("HybridNew: must have 'toysFile' or 'grid' option to have 'readHybridResults'\n");
mass_ = vm["mass"].as<float>();
fullGrid_ = vm.count("fullGrid");
saveGrid_ = vm.count("saveGrid");
fullBToys_ = vm.count("fullBToys");
noUpdateGrid_ = vm.count("noUpdateGrid");
reportPVal_ = vm.count("pvalue");
validateOptions();
}
| void HybridNew::applySignalQuantile | ( | RooStats::HypoTestResult & | hcres | ) | [private] |
Definition at line 1162 of file HybridNew.cc.
References gather_cfg::cout, quantileForExpectedFromGrid_, and python::multivaluedict::sort().
Referenced by applyExpectedQuantile().
{
std::vector<Double_t> stoys = hcres.GetAltDistribution()->GetSamplingDistribution();
std::sort(stoys.begin(), stoys.end());
Double_t testStat = stoys[std::min<int>(floor(quantileForExpectedFromGrid_ * stoys.size()+0.5), stoys.size())];
if (verbose > 0) std::cout << "Text statistics for " << quantileForExpectedFromGrid_ << " quantile: " << testStat << std::endl;
hcres.SetTestStatisticData(testStat);
}
| void HybridNew::clearGrid | ( | ) | [private] |
| std::auto_ptr< RooStats::HybridCalculator > HybridNew::create | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| const RooAbsCollection & | rVals, | ||
| HybridNew::Setup & | setup | ||
| ) | [private] |
generic version
Background-only fit. In 1D case, can use model_s with POI set to zero, but in nD must use model_b.
Definition at line 639 of file HybridNew.cc.
References dtNoiseDBValidation_cfg::cerr, HybridNew::Setup::cleanupList, CLs_, createBeamHaloJobs::constraints, gather_cfg::cout, dqm::qstatus::ERROR, expectedFromGrid_, utils::factorizePdf(), fitNuisances_, fullBToys_, genGlobalObs_, genNuisances_, i, importanceSamplingAlt_, importanceSamplingNull_, instance, testEve_cfg::level, utils::makeNuisancePdf(), MakeSignificance, MakeSignificanceTestStatistics, HybridNew::Setup::modelConfig, HybridNew::Setup::modelConfig_bonly, nCpu_, newToyMCSampler_, nToys_, ProfiledLikelihoodTestStatOpt::oneSidedDef, optimizeProductPdf_, optimizeTestStatistics_, HybridNew::Setup::pc, quantileForExpectedFromGrid_, HybridNew::Setup::qvar, alignCSCRings::r, utils::setAllConstant(), ProfiledLikelihoodTestStatOpt::signFlipDef, testStat_, HybridNew::Setup::toymcsampler, ProfiledLikelihoodTestStatOpt::twoSidedDef, validate_alignment_devdb10_cfg::verbose, withSystematics, and workingMode_.
Referenced by create(), eval(), runSignificance(), runTestStatistics(), updateGridDataFC(), and updateGridPoint().
{
using namespace RooStats;
w->loadSnapshot("clean");
// realData_ = &data;
RooArgSet poi(*mc_s->GetParametersOfInterest()), params(poi);
RooRealVar *r = dynamic_cast<RooRealVar *>(poi.first());
if (poi.getSize() == 1) { // here things are a bit more convoluted, although they could probably be cleaned up
double rVal = ((RooAbsReal*)rVals.first())->getVal();
if (testStat_ != "MLZ") r->setMax(rVal);
r->setVal(rVal);
if (testStat_ == "LHC" || testStat_ == "Profile") {
r->setConstant(false); r->setMin(0);
if (workingMode_ == MakeSignificance || workingMode_ == MakeSignificanceTestStatistics) {
r->setVal(0);
r->removeMax();
}
} else {
r->setConstant(true);
}
} else {
if (testStat_ == "Profile") utils::setAllConstant(poi, false);
else if (testStat_ == "LEP" || testStat_ == "TEV") utils::setAllConstant(poi, true);
poi.assignValueOnly(rVals);
}
//set value of "globalConstrained" nuisances to the value of the corresponding global observable and set them constant
//also fill the RooArgLists which can be passed to the test statistic in order to produce the correct behaviour
//this is only supported currently for the optimized version of the LHC-type test statistic
RooArgList gobsParams;
RooArgList gobs;
if (withSystematics && testStat_ == "LHC" && optimizeTestStatistics_) {
RooArgList allnuis(*mc_s->GetNuisanceParameters());
RooArgList allgobs(*mc_s->GetGlobalObservables());
for (int i=0; i<allnuis.getSize(); ++i) {
RooRealVar *nuis = (RooRealVar*)allnuis.at(i);
if (nuis->getAttribute("globalConstrained")) {
RooRealVar *glob = (RooRealVar*)allgobs.find(TString::Format("%s_In",nuis->GetName()));
if (glob) {
nuis->setVal(glob->getVal());
nuis->setConstant();
gobsParams.add(*nuis);
gobs.add(*glob);
}
}
}
}
std::auto_ptr<RooFitResult> fitMu, fitZero;
if (fitNuisances_ && mc_s->GetNuisanceParameters() && withSystematics) {
TStopwatch timer;
bool isExt = mc_s->GetPdf()->canBeExtended();
utils::setAllConstant(poi, true);
RooAbsPdf *pdfB = mc_b->GetPdf();
if (poi.getSize() == 1) {
r->setVal(0); pdfB = mc_s->GetPdf();
}
timer.Start();
{
CloseCoutSentry sentry(verbose < 3);
fitZero.reset(pdfB->fitTo(data, RooFit::Save(1), RooFit::Minimizer("Minuit2","minimize"), RooFit::Strategy(1), RooFit::Hesse(0), RooFit::Extended(isExt), RooFit::Constrain(*mc_s->GetNuisanceParameters())));
}
if (verbose > 1) { std::cout << "Zero signal fit" << std::endl; fitZero->Print("V"); }
if (verbose > 1) { std::cout << "Fitting of the background hypothesis done in " << timer.RealTime() << " s" << std::endl; }
poi.assignValueOnly(rVals);
timer.Start();
{
CloseCoutSentry sentry(verbose < 3);
fitMu.reset(mc_s->GetPdf()->fitTo(data, RooFit::Save(1), RooFit::Minimizer("Minuit2","minimize"), RooFit::Strategy(1), RooFit::Hesse(0), RooFit::Extended(isExt), RooFit::Constrain(*mc_s->GetNuisanceParameters())));
}
if (verbose > 1) { std::cout << "Reference signal fit" << std::endl; fitMu->Print("V"); }
if (verbose > 1) { std::cout << "Fitting of the signal-plus-background hypothesis done in " << timer.RealTime() << " s" << std::endl; }
} else { fitNuisances_ = false; }
// since ModelConfig cannot allow re-setting sets, we have to re-make everything
setup.modelConfig = ModelConfig("HybridNew_mc_s", w);
setup.modelConfig.SetPdf(*mc_s->GetPdf());
setup.modelConfig.SetObservables(*mc_s->GetObservables());
setup.modelConfig.SetParametersOfInterest(*mc_s->GetParametersOfInterest());
if (withSystematics) {
if (genNuisances_ && mc_s->GetNuisanceParameters()) setup.modelConfig.SetNuisanceParameters(*mc_s->GetNuisanceParameters());
if (genGlobalObs_ && mc_s->GetGlobalObservables()) setup.modelConfig.SetGlobalObservables(*mc_s->GetGlobalObservables());
// if (genGlobalObs_ && mc_s->GetGlobalObservables()) snapGlobalObs_.reset(mc_s->GetGlobalObservables()->snapshot());
}
setup.modelConfig_bonly = ModelConfig("HybridNew_mc_b", w);
setup.modelConfig_bonly.SetPdf(*mc_b->GetPdf());
setup.modelConfig_bonly.SetObservables(*mc_b->GetObservables());
setup.modelConfig_bonly.SetParametersOfInterest(*mc_b->GetParametersOfInterest());
if (withSystematics) {
if (genNuisances_ && mc_b->GetNuisanceParameters()) setup.modelConfig_bonly.SetNuisanceParameters(*mc_b->GetNuisanceParameters());
if (genGlobalObs_ && mc_b->GetGlobalObservables()) setup.modelConfig_bonly.SetGlobalObservables(*mc_b->GetGlobalObservables());
}
if (withSystematics && !genNuisances_) {
// The pdf will contain non-const parameters which are not observables
// and the HybridCalculator will assume they're nuisances and try to generate them
// to avoid this, we force him to generate a fake nuisance instead
if (w->var("__HybridNew_fake_nuis__") == 0) {
w->factory("__HybridNew_fake_nuis__[0.5,0,1]");
w->factory("Uniform::__HybridNew_fake_nuisPdf__(__HybridNew_fake_nuis__)");
}
setup.modelConfig.SetNuisanceParameters(RooArgSet(*w->var("__HybridNew_fake_nuis__")));
setup.modelConfig_bonly.SetNuisanceParameters(RooArgSet(*w->var("__HybridNew_fake_nuis__")));
}
// create snapshots
RooArgSet paramsZero;
if (poi.getSize() == 1) { // in the trivial 1D case, the background has POI=0.
paramsZero.addClone(*rVals.first());
paramsZero.setRealValue(rVals.first()->GetName(), 0);
}
if (fitNuisances_) params.add(fitMu->floatParsFinal());
if (fitNuisances_) paramsZero.addClone(fitZero->floatParsFinal());
setup.modelConfig.SetSnapshot(params);
setup.modelConfig_bonly.SetSnapshot(paramsZero);
TString paramsSnapName = TString::Format("%s_%s_snapshot", setup.modelConfig.GetName(), params.GetName());
TString paramsZSnapName = TString::Format("%s__snapshot", setup.modelConfig_bonly.GetName());
RooFit::MsgLevel level = RooMsgService::instance().globalKillBelow();
RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
w->defineSet(paramsSnapName, params , true);
w->defineSet(paramsZSnapName, paramsZero ,true);
RooMsgService::instance().setGlobalKillBelow(level);
// Create pdfs without nusiances terms, can be used for LEP tests statistics and
// for generating toys when not generating global observables
RooAbsPdf *factorizedPdf_s = setup.modelConfig.GetPdf(), *factorizedPdf_b = setup.modelConfig_bonly.GetPdf();
if (withSystematics && optimizeProductPdf_ && !genGlobalObs_) {
RooArgList constraints;
RooAbsPdf *factorizedPdf_s = utils::factorizePdf(*mc_s->GetObservables(), *mc_s->GetPdf(), constraints);
RooAbsPdf *factorizedPdf_b = utils::factorizePdf(*mc_b->GetObservables(), *mc_b->GetPdf(), constraints);
if (factorizedPdf_s != mc_s->GetPdf()) setup.cleanupList.addOwned(*factorizedPdf_s);
if (factorizedPdf_b != mc_b->GetPdf()) setup.cleanupList.addOwned(*factorizedPdf_b);
setup.modelConfig.SetPdf(*factorizedPdf_s);
setup.modelConfig_bonly.SetPdf(*factorizedPdf_b);
}
if (testStat_ == "LEP") {
//SLR is evaluated using the central value of the nuisance parameters, so we have to put them in the parameter sets
if (withSystematics) {
if (!fitNuisances_) {
params.add(*mc_s->GetNuisanceParameters(), true);
paramsZero.addClone(*mc_b->GetNuisanceParameters(), true);
} else {
std::cerr << "ALERT: using LEP test statistics with --fitNuisances is not validated and most likely broken" << std::endl;
params.assignValueOnly(*mc_s->GetNuisanceParameters());
paramsZero.assignValueOnly(*mc_s->GetNuisanceParameters());
}
}
RooAbsPdf *pdfB = factorizedPdf_b;
if (poi.getSize() == 1) pdfB = factorizedPdf_s; // in this case we can remove the arbitrary constant from the test statistics.
if (optimizeTestStatistics_) {
if (!mc_s->GetPdf()->canBeExtended()) {
setup.qvar.reset(new SimplerLikelihoodRatioTestStat(*pdfB,*factorizedPdf_s, paramsZero, params));
} else {
setup.qvar.reset(new SimplerLikelihoodRatioTestStatOpt(*mc_s->GetObservables(), *pdfB, *factorizedPdf_s, paramsZero, params, withSystematics));
}
} else {
std::cerr << "ALERT: LEP test statistics without optimization not validated." << std::endl;
RooArgSet paramsSnap; params.snapshot(paramsSnap); // needs a snapshot
setup.qvar.reset(new SimpleLikelihoodRatioTestStat(*pdfB,*factorizedPdf_s));
((SimpleLikelihoodRatioTestStat&)*setup.qvar).SetNullParameters(paramsZero); // Null is B
((SimpleLikelihoodRatioTestStat&)*setup.qvar).SetAltParameters(paramsSnap);
}
} else if (testStat_ == "TEV") {
std::cerr << "ALERT: Tevatron test statistics not yet validated." << std::endl;
RooAbsPdf *pdfB = factorizedPdf_b;
if (poi.getSize() == 1) pdfB = factorizedPdf_s; // in this case we can remove the arbitrary constant from the test statistics.
if (optimizeTestStatistics_) {
setup.qvar.reset(new ProfiledLikelihoodRatioTestStatOpt(*mc_s->GetObservables(), *pdfB, *mc_s->GetPdf(), mc_s->GetNuisanceParameters(), paramsZero, params));
((ProfiledLikelihoodRatioTestStatOpt&)*setup.qvar).setPrintLevel(verbose);
} else {
setup.qvar.reset(new RatioOfProfiledLikelihoodsTestStat(*mc_s->GetPdf(), *pdfB, setup.modelConfig.GetSnapshot()));
((RatioOfProfiledLikelihoodsTestStat&)*setup.qvar).SetSubtractMLE(false);
}
} else if (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile") {
if (poi.getSize() != 1 && testStat_ != "Profile") {
throw std::logic_error("ERROR: modified profile likelihood definitions (LHC, LHCFC) do not make sense in more than one dimension");
}
if (optimizeTestStatistics_) {
ProfiledLikelihoodTestStatOpt::OneSidedness side = ProfiledLikelihoodTestStatOpt::oneSidedDef;
if (testStat_ == "LHCFC") side = ProfiledLikelihoodTestStatOpt::signFlipDef;
if (testStat_ == "Profile") side = ProfiledLikelihoodTestStatOpt::twoSidedDef;
if (workingMode_ == MakeSignificance) r->setVal(0.0);
setup.qvar.reset(new ProfiledLikelihoodTestStatOpt(*mc_s->GetObservables(), *mc_s->GetPdf(), mc_s->GetNuisanceParameters(), params, poi, gobsParams,gobs, verbose, side));
} else {
std::cerr << "ALERT: LHC test statistics without optimization not validated." << std::endl;
setup.qvar.reset(new ProfileLikelihoodTestStat(*mc_s->GetPdf()));
if (testStat_ == "LHC") {
((ProfileLikelihoodTestStat&)*setup.qvar).SetOneSided(true);
} else if (testStat_ == "LHCFC") {
throw std::invalid_argument("Test statistics LHCFC is not supported without optimization");
}
}
} else if (testStat_ == "MLZ") {
if (workingMode_ == MakeSignificance) r->setVal(0.0);
setup.qvar.reset(new BestFitSigmaTestStat(*mc_s->GetObservables(), *mc_s->GetPdf(), mc_s->GetNuisanceParameters(), params, verbose));
}
RooAbsPdf *nuisancePdf = 0;
if (withSystematics && (genNuisances_ || (newToyMCSampler_ && genGlobalObs_))) {
nuisancePdf = utils::makeNuisancePdf(*mc_s);
if (nuisancePdf) setup.cleanupList.addOwned(*nuisancePdf);
}
if (newToyMCSampler_) {
setup.toymcsampler.reset(new ToyMCSamplerOpt(*setup.qvar, nToys_, nuisancePdf, genNuisances_));
} else {
std::cerr << "ALERT: running with newToyMC=0 not validated." << std::endl;
setup.toymcsampler.reset(new ToyMCSampler(*setup.qvar, nToys_));
}
if (!mc_b->GetPdf()->canBeExtended()) setup.toymcsampler->SetNEventsPerToy(1);
if (nCpu_ > 0) {
std::cerr << "ALERT: running with proof not validated." << std::endl;
if (verbose > 1) std::cout << " Will use " << nCpu_ << " CPUs." << std::endl;
setup.pc.reset(new ProofConfig(*w, nCpu_, "", kFALSE));
setup.toymcsampler->SetProofConfig(setup.pc.get());
}
std::auto_ptr<HybridCalculator> hc(new HybridCalculator(data, setup.modelConfig, setup.modelConfig_bonly, setup.toymcsampler.get()));
if (genNuisances_ || !genGlobalObs_) {
if (withSystematics) {
setup.toymcsampler->SetGlobalObservables(*setup.modelConfig.GetNuisanceParameters());
(static_cast<HybridCalculator&>(*hc)).ForcePriorNuisanceNull(*nuisancePdf);
(static_cast<HybridCalculator&>(*hc)).ForcePriorNuisanceAlt(*nuisancePdf);
}
} else if (genGlobalObs_ && !genNuisances_) {
setup.toymcsampler->SetGlobalObservables(*setup.modelConfig.GetGlobalObservables());
hc->ForcePriorNuisanceNull(*w->pdf("__HybridNew_fake_nuisPdf__"));
hc->ForcePriorNuisanceAlt(*w->pdf("__HybridNew_fake_nuisPdf__"));
}
// we need less B toys than S toys
if (workingMode_ == MakeSignificance) {
// need only B toys. just keep a few S+B ones to avoid possible divide-by-zero errors somewhere
hc->SetToys(nToys_, int(0.01*nToys_)+1);
if (fullBToys_) {
hc->SetToys(nToys_, nToys_);
}
} else if (!CLs_) {
// we need only S+B toys to compute CLs+b
hc->SetToys(fullBToys_ ? nToys_ : int(0.01*nToys_)+1, nToys_);
//for two sigma bands need an equal number of B
if (expectedFromGrid_ && (fabs(0.5-quantileForExpectedFromGrid_)>=0.4) ) {
hc->SetToys(nToys_, nToys_);
}
} else {
// need both, but more S+B than B
hc->SetToys(fullBToys_ ? nToys_ : int(0.25*nToys_), nToys_);
//for two sigma bands need an equal number of B
if (expectedFromGrid_ && (fabs(0.5-quantileForExpectedFromGrid_)>=0.4) ) {
hc->SetToys(nToys_, nToys_);
}
}
static const char * istr = "__HybridNew__importanceSamplingDensity";
if(importanceSamplingNull_) {
std::cerr << "ALERT: running with importance sampling not validated (and probably not working)." << std::endl;
if(verbose > 1) std::cout << ">>> Enabling importance sampling for null hyp." << std::endl;
if(!withSystematics) {
throw std::invalid_argument("Importance sampling is not available without systematics");
}
RooArgSet importanceSnapshot;
importanceSnapshot.addClone(poi);
importanceSnapshot.addClone(*mc_s->GetNuisanceParameters());
if (verbose > 2) importanceSnapshot.Print("V");
hc->SetNullImportanceDensity(mc_b->GetPdf(), &importanceSnapshot);
}
if(importanceSamplingAlt_) {
std::cerr << "ALERT: running with importance sampling not validated (and probably not working)." << std::endl;
if(verbose > 1) std::cout << ">>> Enabling importance sampling for alt. hyp." << std::endl;
if(!withSystematics) {
throw std::invalid_argument("Importance sampling is not available without systematics");
}
if (w->pdf(istr) == 0) {
w->factory("__oneHalf__[0.5]");
RooAddPdf *sum = new RooAddPdf(istr, "fifty-fifty", *mc_s->GetPdf(), *mc_b->GetPdf(), *w->var("__oneHalf__"));
w->import(*sum);
}
RooArgSet importanceSnapshot;
importanceSnapshot.addClone(poi);
importanceSnapshot.addClone(*mc_s->GetNuisanceParameters());
if (verbose > 2) importanceSnapshot.Print("V");
hc->SetAltImportanceDensity(w->pdf(istr), &importanceSnapshot);
}
return hc;
}
| std::auto_ptr< RooStats::HybridCalculator > HybridNew::create | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| double | rVal, | ||
| HybridNew::Setup & | setup | ||
| ) | [private] |
specific version used in unidimensional searches (calls the generic one)
Definition at line 630 of file HybridNew.cc.
References create(), and alignCSCRings::r.
| std::pair< double, double > HybridNew::eval | ( | const RooStats::HypoTestResult & | hcres, |
| double | rVal | ||
| ) | [private] |
Definition at line 1020 of file HybridNew.cc.
References CLs_, data, EPS, i, max(), minimizerTolerance_, n, query::result, rule_, and testStat_.
{
if (testStat_ == "LHCFC") {
RooStats::SamplingDistribution * bDistribution = hcres.GetNullDistribution(), * sDistribution = hcres.GetAltDistribution();
const std::vector<Double_t> & bdist = bDistribution->GetSamplingDistribution();
const std::vector<Double_t> & bweight = bDistribution->GetSampleWeights();
const std::vector<Double_t> & sdist = sDistribution->GetSamplingDistribution();
const std::vector<Double_t> & sweight = sDistribution->GetSampleWeights();
Double_t data = hcres.GetTestStatisticData();
std::vector<Double_t> absbdist(bdist.size()), abssdist(sdist.size());
std::vector<Double_t> absbweight(bweight), abssweight(sweight);
Double_t absdata;
if (rule_ == "FC") {
for (int i = 0, n = absbdist.size(); i < n; ++i) absbdist[i] = fabs(bdist[i]);
for (int i = 0, n = abssdist.size(); i < n; ++i) abssdist[i] = fabs(sdist[i]);
absdata = fabs(data)-2*minimizerTolerance_; // needed here since zeros are not exact by more than tolerance
} else {
for (int i = 0, n = absbdist.size(); i < n; ++i) absbdist[i] = max(0., bdist[i]);
for (int i = 0, n = abssdist.size(); i < n; ++i) abssdist[i] = max(0., sdist[i]);
absdata = max(0., data) - EPS;
}
if (rVal == 0) { // S+B toys are equal to B ones!
abssdist.reserve(absbdist.size() + abssdist.size());
abssdist.insert(abssdist.end(), absbdist.begin(), absbdist.end());
abssweight.reserve(absbweight.size() + abssweight.size());
abssweight.insert(abssweight.end(), absbweight.begin(), absbweight.end());
}
RooStats::HypoTestResult result;
RooStats::SamplingDistribution *abssDist = new RooStats::SamplingDistribution("s","s",abssdist,abssweight);
RooStats::SamplingDistribution *absbDist = new RooStats::SamplingDistribution("b","b",absbdist,absbweight);
result.SetNullDistribution(absbDist);
result.SetAltDistribution(abssDist);
result.SetTestStatisticData(absdata);
result.SetPValueIsRightTail(!result.GetPValueIsRightTail());
if (CLs_) {
return std::pair<double,double>(result.CLs(), result.CLsError());
} else {
return std::pair<double,double>(result.CLsplusb(), result.CLsplusbError());
}
} else {
if (CLs_) {
return std::pair<double,double>(hcres.CLs(), hcres.CLsError());
} else {
return std::pair<double,double>(hcres.CLsplusb(), hcres.CLsplusbError());
}
}
}
| std::pair< double, double > HybridNew::eval | ( | const RooStats::HypoTestResult & | hcres, |
| const RooAbsCollection & | rVals | ||
| ) | [private] |
Definition at line 1014 of file HybridNew.cc.
References eval().
{
double rVal = ((RooAbsReal*)rVals.first())->getVal();
return eval(hcres,rVal);
}
| std::pair< double, double > HybridNew::eval | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| const RooAbsCollection & | rVals, | ||
| bool | adaptive = false, |
||
| double | clsTarget = -1 |
||
| ) | [private] |
generic version
Definition at line 587 of file HybridNew.cc.
References applyExpectedQuantile(), gather_cfg::cout, create(), EPS, expectedFromGrid_, limitPlot_, HybridNew::Setup::modelConfig_bonly, noUpdateGrid_, HybridNew::Setup::qvar, alignCSCRings::r, readHybridResults_, readToysFromFile(), query::result, run_regression::ret, HcalObjRepresent::setup(), and testStat_.
Referenced by eval(), runLimit(), runSinglePoint(), updateGridPoint(), and useGrid().
{
if (readHybridResults_) {
bool isProfile = (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile");
std::auto_ptr<RooStats::HypoTestResult> result(readToysFromFile(rVals));
std::pair<double, double> ret(-1,-1);
assert(result.get() != 0 && "Null result in HybridNew::eval(Workspace,...) after readToysFromFile");
if (expectedFromGrid_) {
applyExpectedQuantile(*result);
result->SetTestStatisticData(result->GetTestStatisticData() + (isProfile ? -EPS : EPS));
} else if (!noUpdateGrid_) {
Setup setup;
std::auto_ptr<RooStats::HybridCalculator> hc = create(w, mc_s, mc_b, data, rVals, setup);
RooArgSet nullPOI(*setup.modelConfig_bonly.GetSnapshot());
if (isProfile) nullPOI.assignValueOnly(rVals);
double testStat = setup.qvar->Evaluate(data, nullPOI);
result->SetTestStatisticData(testStat + (isProfile ? -EPS : EPS));
}
ret = eval(*result, rVals);
return ret;
}
HybridNew::Setup setup;
RooLinkedListIter it = rVals.iterator();
for (RooRealVar *rIn = (RooRealVar*) it.Next(); rIn != 0; rIn = (RooRealVar*) it.Next()) {
RooRealVar *r = dynamic_cast<RooRealVar *>(mc_s->GetParametersOfInterest()->find(rIn->GetName()));
r->setVal(rIn->getVal());
if (verbose) std::cout << " " << r->GetName() << " = " << rIn->getVal() << " +/- " << r->getError() << std::endl;
}
std::auto_ptr<RooStats::HybridCalculator> hc(create(w, mc_s, mc_b, data, rVals, setup));
std::pair<double, double> ret = eval(*hc, rVals, adaptive, clsTarget);
// add to plot
if (limitPlot_.get()) {
limitPlot_->Set(limitPlot_->GetN()+1);
limitPlot_->SetPoint(limitPlot_->GetN()-1, ((RooAbsReal*)rVals.first())->getVal(), ret.first);
limitPlot_->SetPointError(limitPlot_->GetN()-1, 0, ret.second);
}
return ret;
}
| std::pair< double, double > HybridNew::eval | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| double | rVal, | ||
| bool | adaptive = false, |
||
| double | clsTarget = -1 |
||
| ) | [private] |
specific version used in unidimensional searches (calls the generic one)
Definition at line 578 of file HybridNew.cc.
References eval(), and alignCSCRings::r.
| std::pair< double, double > HybridNew::eval | ( | RooStats::HybridCalculator & | hc, |
| const RooAbsCollection & | rVals, | ||
| bool | adaptive = false, |
||
| double | clsTarget = -1 |
||
| ) | [private] |
Definition at line 935 of file HybridNew.cc.
References applyExpectedQuantile(), alignmentValidation::c1, dtNoiseDBValidation_cfg::cerr, CLs_, clsAccuracy_, gather_cfg::cout, EPS, eval(), evalGeneric(), expectedFromGrid_, i, if(), iterations_, MakeLimit, mass_, max(), name(), nToys_, perf_totalToysRun_, plotResiduals::plot(), plot_, quantileForExpectedFromGrid_, saveHybridResult_, testStat_, workingMode_, and writeToysHere.
{
std::auto_ptr<HypoTestResult> hcResult(evalGeneric(hc));
if (expectedFromGrid_) applyExpectedQuantile(*hcResult);
if (hcResult.get() == 0) {
std::cerr << "Hypotest failed" << std::endl;
return std::pair<double, double>(-1,-1);
}
if (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile") {
// I need to flip the P-values
hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()-EPS); // issue with < vs <= in discrete models
} else {
hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()+EPS); // issue with < vs <= in discrete models
hcResult->SetPValueIsRightTail(!hcResult->GetPValueIsRightTail());
}
std::pair<double,double> cls = eval(*hcResult, rVals);
if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tCLsplusb = ") << cls.first << " +/- " << cls.second << std::endl;
if (adaptive) {
if (CLs_) {
hc.SetToys(int(0.25*nToys_ + 1), nToys_);
}
else {
hc.SetToys(1, nToys_);
}
//for two sigma bands need an equal number of B
if (expectedFromGrid_ && (fabs(0.5-quantileForExpectedFromGrid_)>=0.4) ) {
hc.SetToys(nToys_, nToys_);
}
while (cls.second >= clsAccuracy_ && (clsTarget == -1 || fabs(cls.first-clsTarget) < 3*cls.second) ) {
std::auto_ptr<HypoTestResult> more(evalGeneric(hc));
more->SetBackgroundAsAlt(false);
if (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile") more->SetPValueIsRightTail(!more->GetPValueIsRightTail());
hcResult->Append(more.get());
if (expectedFromGrid_) applyExpectedQuantile(*hcResult);
cls = eval(*hcResult, rVals);
if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tCLsplusb = ") << cls.first << " +/- " << cls.second << std::endl;
}
} else if (iterations_ > 1) {
for (unsigned int i = 1; i < iterations_; ++i) {
std::auto_ptr<HypoTestResult> more(evalGeneric(hc));
more->SetBackgroundAsAlt(false);
if (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile") more->SetPValueIsRightTail(!more->GetPValueIsRightTail());
hcResult->Append(more.get());
if (expectedFromGrid_) applyExpectedQuantile(*hcResult);
cls = eval(*hcResult, rVals);
if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tCLsplusb = ") << cls.first << " +/- " << cls.second << std::endl;
}
}
if (verbose > 0) {
std::cout <<
"\tCLs = " << hcResult->CLs() << " +/- " << hcResult->CLsError() << "\n" <<
"\tCLb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << "\n" <<
"\tCLsplusb = " << hcResult->CLsplusb() << " +/- " << hcResult->CLsplusbError() << "\n" <<
std::endl;
}
perf_totalToysRun_ += (hcResult->GetAltDistribution()->GetSize() + hcResult->GetNullDistribution()->GetSize());
if (!plot_.empty() && workingMode_ != MakeLimit) {
HypoTestPlot plot(*hcResult, 30);
TCanvas *c1 = new TCanvas("c1","c1");
plot.Draw();
c1->Print(plot_.c_str());
delete c1;
}
if (saveHybridResult_) {
TString name = TString::Format("HypoTestResult_mh%g",mass_);
RooLinkedListIter it = rVals.iterator();
for (RooRealVar *rIn = (RooRealVar*) it.Next(); rIn != 0; rIn = (RooRealVar*) it.Next()) {
name += Form("_%s%g", rIn->GetName(), rIn->getVal());
}
name += Form("_%u", RooRandom::integer(std::numeric_limits<UInt_t>::max() - 1));
writeToysHere->WriteTObject(new HypoTestResult(*hcResult), name);
if (verbose) std::cout << "Hybrid result saved as " << name << " in " << writeToysHere->GetFile()->GetName() << " : " << writeToysHere->GetPath() << std::endl;
}
return cls;
}
| RooStats::HypoTestResult * HybridNew::evalGeneric | ( | RooStats::HybridCalculator & | hc, |
| bool | forceNoFork = false |
||
| ) | [private] |
Definition at line 1170 of file HybridNew.cc.
References gather_cfg::cout, evalWithFork(), fork_, runtimedef::get(), and run_regression::ret.
Referenced by eval(), evalWithFork(), and runSignificance().
{
if (fork_ && !noFork) return evalWithFork(hc);
else {
TStopwatch timer; timer.Start();
RooStats::HypoTestResult * ret = hc.GetHypoTest();
if (runtimedef::get("HybridNew_Timing")) std::cout << "Evaluated toys in " << timer.RealTime() << " s " << std::endl;
return ret;
}
}
| RooStats::HypoTestResult * HybridNew::evalWithFork | ( | RooStats::HybridCalculator & | hc | ) | [private] |
Definition at line 1180 of file HybridNew.cc.
References gather_cfg::cout, evalGeneric(), f, fork_, max(), NULL, query::result, run_regression::ret, and createPayload::tmpfile.
Referenced by evalGeneric().
{
TStopwatch timer;
std::auto_ptr<RooStats::HypoTestResult> result(0);
char *tmpfile = tempnam(NULL,"rstat");
unsigned int ich = 0;
std::vector<UInt_t> newSeeds(fork_);
for (ich = 0; ich < fork_; ++ich) {
newSeeds[ich] = RooRandom::integer(std::numeric_limits<UInt_t>::max()-1);
if (!fork()) break; // spawn children (but only in the parent thread)
}
if (ich == fork_) { // if i'm the parent
int cstatus, ret;
do {
do { ret = waitpid(-1, &cstatus, 0); } while (ret == -1 && errno == EINTR);
} while (ret != -1);
if (ret == -1 && errno != ECHILD) throw std::runtime_error("Didn't wait for child");
for (ich = 0; ich < fork_; ++ich) {
TFile *f = TFile::Open(TString::Format("%s.%d.root", tmpfile, ich));
if (f == 0) throw std::runtime_error(TString::Format("Child didn't leave output file %s.%d.root", tmpfile, ich).Data());
RooStats::HypoTestResult *res = dynamic_cast<RooStats::HypoTestResult *>(f->Get("result"));
if (res == 0) throw std::runtime_error(TString::Format("Child output file %s.%d.root is corrupted", tmpfile, ich).Data());
if (result.get()) result->Append(res); else result.reset(dynamic_cast<RooStats::HypoTestResult *>(res->Clone()));
f->Close();
unlink(TString::Format("%s.%d.root", tmpfile, ich).Data());
unlink(TString::Format("%s.%d.out.txt", tmpfile, ich).Data());
unlink(TString::Format("%s.%d.err.txt", tmpfile, ich).Data());
}
} else {
RooRandom::randomGenerator()->SetSeed(newSeeds[ich]);
freopen(TString::Format("%s.%d.out.txt", tmpfile, ich).Data(), "w", stdout);
freopen(TString::Format("%s.%d.err.txt", tmpfile, ich).Data(), "w", stderr);
std::cout << " I'm child " << ich << ", seed " << newSeeds[ich] << std::endl;
RooStats::HypoTestResult *hcResult = evalGeneric(hc, /*noFork=*/true);
TFile *f = TFile::Open(TString::Format("%s.%d.root", tmpfile, ich), "RECREATE");
f->WriteTObject(hcResult, "result");
f->ls();
f->Close();
fflush(stdout); fflush(stderr);
std::cout << "And I'm done" << std::endl;
throw std::runtime_error("done"); // I have to throw instead of exiting, otherwise there's no proper stack unwinding
// and deleting of intermediate objects, and when the statics get deleted it crashes
// in 5.27.06 (but not in 5.28)
}
free(tmpfile);
if (verbose > 1) { std::cout << " Evaluation of p-values done in " << timer.RealTime() << " s" << std::endl; }
return result.release();
}
| virtual const std::string& HybridNew::name | ( | void | ) | const [inline, virtual] |
Implements LimitAlgo.
Definition at line 32 of file HybridNew.h.
Referenced by eval(), readGrid(), and runSignificance().
| void HybridNew::readGrid | ( | TDirectory * | directory, |
| double | rMin, | ||
| double | rMax | ||
| ) | [private] |
Definition at line 1345 of file HybridNew.cc.
References clearGrid(), gather_cfg::cout, grid_, gen::k, mass_, relval_steps::merge(), name(), GetRecoTauVFromDQM_MC_cff::next, point, and rValues_.
Referenced by runLimit().
{
if (rValues_.getSize() != 1) throw std::runtime_error("Running limits with grid only works in one dimension for the moment");
clearGrid();
TIter next(toyDir->GetListOfKeys()); TKey *k; const char *poiName = rValues_.first()->GetName();
while ((k = (TKey *) next()) != 0) {
TString name(k->GetName());
if (name.Index("HypoTestResult_mh") == 0) {
if (name.Index(TString::Format("HypoTestResult_mh%g_%s",mass_,poiName)) != 0 || name.Index("_", name.Index("_")+1) == -1) continue;
name.ReplaceAll(TString::Format("HypoTestResult_mh%g_%s",mass_,poiName),""); // remove the prefix
if (name.Index("_") == -1) continue; // check if it has the _<uniqueId> postfix
name.Remove(name.Index("_"),name.Length()); // remove it before calling atof
} else if (name.Index("HypoTestResult_") == 0) {
// let's put a warning here, since results of this form were supported in the past
std::cout << "HybridNew::readGrid: HypoTestResult with non-conformant name " << name << " will be skipped" << std::endl;
continue;
} else continue;
double rVal = atof(name.Data());
if (rVal < rMin || rVal > rMax) continue;
if (verbose > 2) std::cout << " Do " << k->GetName() << " -> " << name << " --> " << rVal << std::endl;
RooStats::HypoTestResult *toy = dynamic_cast<RooStats::HypoTestResult *>(toyDir->Get(k->GetName()));
RooStats::HypoTestResult *&merge = grid_[rVal];
if (merge == 0) merge = new RooStats::HypoTestResult(*toy);
else merge->Append(toy);
merge->ResetBit(1);
}
if (verbose > 1) {
std::cout << "GRID, as is." << std::endl;
typedef std::map<double, RooStats::HypoTestResult *>::iterator point;
for (point it = grid_.begin(), ed = grid_.end(); it != ed; ++it) {
std::cout << " - " << it->first << " (CLs = " << it->second->CLs() << " +/- " << it->second->CLsError() << ")" << std::endl;
}
}
}
| RooStats::HypoTestResult * HybridNew::readToysFromFile | ( | const RooAbsCollection & | rVals | ) | [private] |
Definition at line 1293 of file HybridNew.cc.
References alignmentValidation::c1, gather_cfg::cout, gen::k, MakeLimit, mass_, GetRecoTauVFromDQM_MC_cff::next, plotResiduals::plot(), plot_, readToysFromHere, run_regression::ret, and workingMode_.
Referenced by eval(), runSignificance(), and runTestStatistics().
{
if (!readToysFromHere) throw std::logic_error("Cannot use readHypoTestResult: option toysFile not specified, or input file empty");
TDirectory *toyDir = readToysFromHere->GetDirectory("toys");
if (!toyDir) throw std::logic_error("Cannot use readHypoTestResult: empty toy dir in input file empty");
if (verbose) std::cout << "Reading toys for ";
TString prefix1 = TString::Format("HypoTestResult_mh%g",mass_);
TString prefix2 = TString::Format("HypoTestResult");
RooLinkedListIter it = rVals.iterator();
for (RooRealVar *rIn = (RooRealVar*) it.Next(); rIn != 0; rIn = (RooRealVar*) it.Next()) {
if (verbose) std::cout << rIn->GetName() << " = " << rIn->getVal() << " ";
prefix1 += Form("_%s%g", rIn->GetName(), rIn->getVal());
prefix2 += Form("_%s%g", rIn->GetName(), rIn->getVal());
}
if (verbose) std::cout << std::endl;
std::auto_ptr<RooStats::HypoTestResult> ret;
TIter next(toyDir->GetListOfKeys()); TKey *k;
while ((k = (TKey *) next()) != 0) {
if (TString(k->GetName()).Index(prefix1) != 0 && TString(k->GetName()).Index(prefix2) != 0) continue;
RooStats::HypoTestResult *toy = dynamic_cast<RooStats::HypoTestResult *>(toyDir->Get(k->GetName()));
if (toy == 0) continue;
if (verbose > 1) std::cout << " - " << k->GetName() << std::endl;
if (ret.get() == 0) {
ret.reset(new RooStats::HypoTestResult(*toy));
} else {
ret->Append(toy);
}
}
if (ret.get() == 0) {
std::cout << "ERROR: parameter point not found in input root file.\n";
rVals.Print("V");
if (verbose > 0) toyDir->ls();
std::cout << "ERROR: parameter point not found in input root file" << std::endl;
throw std::invalid_argument("Missing input");
}
if (verbose > 0) {
std::cout <<
"\tCLs = " << ret->CLs() << " +/- " << ret->CLsError() << "\n" <<
"\tCLb = " << ret->CLb() << " +/- " << ret->CLbError() << "\n" <<
"\tCLsplusb = " << ret->CLsplusb() << " +/- " << ret->CLsplusbError() << "\n" <<
std::endl;
if (!plot_.empty() && workingMode_ != MakeLimit) {
HypoTestPlot plot(*ret, 30);
TCanvas *c1 = new TCanvas("c1","c1");
plot.Draw();
c1->Print(plot_.c_str());
delete c1;
}
}
return ret.release();
}
| bool HybridNew::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 251 of file HybridNew.cc.
References DEBUG, MakeLimit, MakePValues, MakeSignificance, MakeSignificanceTestStatistics, MakeTestStatistics, minimizerAlgo_, minimizerTolerance_, perf_totalToysRun_, runLimit(), runSignificance(), runSinglePoint(), runTestStatistics(), rValues_, setupPOI(), dqm::qstatus::WARNING, and workingMode_.
{
RooFitGlobalKillSentry silence(verbose <= 1 ? RooFit::WARNING : RooFit::DEBUG);
ProfileLikelihood::MinimizerSentry minimizerConfig(minimizerAlgo_, minimizerTolerance_);
perf_totalToysRun_ = 0; // reset performance counter
if (rValues_.getSize() == 0) setupPOI(mc_s);
switch (workingMode_) {
case MakeLimit: return runLimit(w, mc_s, mc_b, data, limit, limitErr, hint);
case MakeSignificance: return runSignificance(w, mc_s, mc_b, data, limit, limitErr, hint);
case MakePValues: return runSinglePoint(w, mc_s, mc_b, data, limit, limitErr, hint);
case MakeTestStatistics:
case MakeSignificanceTestStatistics:
return runTestStatistics(w, mc_s, mc_b, data, limit, limitErr, hint);
}
assert("Shouldn't get here" == 0);
}
| bool HybridNew::runLimit | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| double & | limit, | ||
| double & | limitErr, | ||
| const double * | hint | ||
| ) | [virtual] |
Definition at line 318 of file HybridNew.cc.
References algo_, alignmentValidation::c1, dtNoiseDBValidation_cfg::cerr, cl, CLs_, Combine::commitPoint(), gather_cfg::cout, run_regression::done, eval(), fullGrid_, grid_, gridFile_, i, interpAccuracy_, j, MessageLogger_cff::limit, limitPlot_, geometryCSVtoXML::line, create_public_lumi_plots::log, lowerLimit_, max(), min, n, noUpdateGrid_, perf_totalToysRun_, plot_, alignCSCRings::r, rAbsAccuracy_, readGrid(), readHybridResults_, rMaxSet_, rMinSet_, rRelAccuracy_, rule_, saveGrid_, testStat_, updateGridData(), updateGridDataFC(), useGrid(), x, and detailsBasic3DVector::y.
Referenced by run().
{
if (mc_s->GetParametersOfInterest()->getSize() != 1) throw std::logic_error("Cannot run limits in more than one dimension, for the moment.");
RooRealVar *r = dynamic_cast<RooRealVar *>(mc_s->GetParametersOfInterest()->first()); r->setConstant(true);
w->loadSnapshot("clean");
if ((hint != 0) && (*hint > r->getMin())) {
r->setMax(std::min<double>(3.0 * (*hint), r->getMax()));
r->setMin(std::max<double>(0.3 * (*hint), r->getMin()));
}
typedef std::pair<double,double> CLs_t;
double clsTarget = 1 - cl;
CLs_t clsMin(1,0), clsMax(0,0), clsMid(0,0);
double rMin = r->getMin(), rMax = r->getMax();
limit = 0.5*(rMax + rMin);
limitErr = 0.5*(rMax - rMin);
bool done = false;
TF1 expoFit("expoFit","[0]*exp([1]*(x-[2]))", rMin, rMax);
if (readHybridResults_) {
if (verbose > 0) std::cout << "Search for upper limit using pre-computed grid of p-values" << std::endl;
if (!gridFile_.empty()) {
if (grid_.empty()) {
std::auto_ptr<TFile> gridFile(TFile::Open(gridFile_.c_str()));
if (gridFile.get() == 0) throw std::runtime_error(("Can't open grid file "+gridFile_).c_str());
TDirectory *toyDir = gridFile->GetDirectory("toys");
if (!toyDir) throw std::logic_error("Cannot use readHypoTestResult: empty toy dir in input file empty");
readGrid(toyDir, rMinSet_ ? rMin : -99e99, rMaxSet_ ? rMax :+99e99);
}
if (grid_.size() <= 1) throw std::logic_error("The grid must contain at least 2 points.");
if (noUpdateGrid_) {
if (testStat_ == "LHCFC" && rule_ == "FC" && (saveGrid_ || lowerLimit_)) {
std::cout << "Will have to re-run points for which the test statistics was set to zero" << std::endl;
updateGridDataFC(w, mc_s, mc_b, data, !fullGrid_, clsTarget);
} else {
std::cout << "Will use the test statistics that had already been computed" << std::endl;
}
} else {
updateGridData(w, mc_s, mc_b, data, !fullGrid_, clsTarget);
}
} else throw std::logic_error("When setting a limit reading results from file, a grid file must be specified with option --grid");
if (grid_.size() <= 1) throw std::logic_error("The grid must contain at least 2 points.");
useGrid();
double minDist=1e3;
double nearest = 0;
rMin = limitPlot_->GetX()[0];
rMax = limitPlot_->GetX()[limitPlot_->GetN()-1];
for (int i = 0, n = limitPlot_->GetN(); i < n; ++i) {
double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i], ey = limitPlot_->GetErrorY(i);
if (verbose > 0) printf(" r %.2f: %s = %6.4f +/- %6.4f\n", x, CLs_ ? "CLs" : "CLsplusb", y, ey);
if (saveGrid_) { limit = x; limitErr = ey; Combine::commitPoint(false, y); }
if (y-3*max(ey,0.01) >= clsTarget) { rMin = x; clsMin = CLs_t(y,ey); }
if (fabs(y-clsTarget) < minDist) { nearest = x; minDist = fabs(y-clsTarget); }
rMax = x; clsMax = CLs_t(y,ey);
if (y+3*max(ey,0.01) <= clsTarget && !fullGrid_) break;
}
limit = nearest;
if (verbose > 0) std::cout << " after scan x ~ " << limit << ", bounds [ " << rMin << ", " << rMax << "]" << std::endl;
limitErr = std::max(limit-rMin, rMax-limit);
expoFit.SetRange(rMin,rMax);
if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
if (verbose > 1) std::cout << " reached accuracy " << limitErr << " below " << std::max(rAbsAccuracy_, rRelAccuracy_ * limit) << std::endl;
done = true;
}
} else {
limitPlot_.reset(new TGraphErrors());
if (verbose > 0) std::cout << "Search for upper limit to the limit" << std::endl;
for (int tries = 0; tries < 6; ++tries) {
clsMax = eval(w, mc_s, mc_b, data, rMax);
if (lowerLimit_) break; // we can't search for lower limits this way
if (clsMax.first == 0 || clsMax.first + 3 * fabs(clsMax.second) < clsTarget ) break;
rMax += rMax;
if (tries == 5) {
std::cerr << "Cannot set higher limit: at " << r->GetName() << " = " << rMax << " still get " << (CLs_ ? "CLs" : "CLsplusb") << " = " << clsMax.first << std::endl;
return false;
}
}
if (verbose > 0) std::cout << "Search for lower limit to the limit" << std::endl;
clsMin = (CLs_ && rMin == 0 ? CLs_t(1,0) : eval(w, mc_s, mc_b, data, rMin));
if (!lowerLimit_ && clsMin.first != 1 && clsMin.first - 3 * fabs(clsMin.second) < clsTarget) {
if (CLs_) {
rMin = 0;
clsMin = CLs_t(1,0); // this is always true for CLs
} else {
rMin = -rMax / 4;
for (int tries = 0; tries < 6; ++tries) {
clsMin = eval(w, mc_s, mc_b, data, rMin);
if (clsMin.first == 1 || clsMin.first - 3 * fabs(clsMin.second) > clsTarget) break;
rMin += rMin;
if (tries == 5) {
std::cerr << "Cannot set lower limit: at " << r->GetName() << " = " << rMin << " still get " << (CLs_ ? "CLs" : "CLsplusb") << " = " << clsMin.first << std::endl;
return false;
}
}
}
}
if (verbose > 0) std::cout << "Now doing proper bracketing & bisection" << std::endl;
do {
// determine point by bisection or interpolation
limit = 0.5*(rMin+rMax); limitErr = 0.5*(rMax-rMin);
if (algo_ == "logSecant" && clsMax.first != 0 && clsMin.first != 0) {
double logMin = log(clsMin.first), logMax = log(clsMax.first), logTarget = log(clsTarget);
limit = rMin + (rMax-rMin) * (logTarget - logMin)/(logMax - logMin);
if (clsMax.second != 0 && clsMin.second != 0) {
limitErr = hypot((logTarget-logMax) * (clsMin.second/clsMin.first), (logTarget-logMin) * (clsMax.second/clsMax.first));
limitErr *= (rMax-rMin)/((logMax-logMin)*(logMax-logMin));
}
// disallow "too precise" interpolations
if (limitErr < interpAccuracy_ * (rMax-rMin)) limitErr = interpAccuracy_ * (rMax-rMin);
}
r->setError(limitErr);
// exit if reached accuracy on r
if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
if (verbose > 1) std::cout << " reached accuracy " << limitErr << " below " << std::max(rAbsAccuracy_, rRelAccuracy_ * limit) << std::endl;
done = true; break;
}
// evaluate point
clsMid = eval(w, mc_s, mc_b, data, limit, true, clsTarget);
if (clsMid.second == -1) {
std::cerr << "Hypotest failed" << std::endl;
return false;
}
// if sufficiently far away, drop one of the points
if (fabs(clsMid.first-clsTarget) >= 2*clsMid.second) {
if ((clsMid.first>clsTarget) == (clsMax.first>clsTarget)) {
rMax = limit; clsMax = clsMid;
} else {
rMin = limit; clsMin = clsMid;
}
} else {
if (verbose > 0) std::cout << "Trying to move the interval edges closer" << std::endl;
double rMinBound = rMin, rMaxBound = rMax;
// try to reduce the size of the interval
while (clsMin.second == 0 || fabs(rMin-limit) > std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
rMin = 0.5*(rMin+limit);
clsMin = eval(w, mc_s, mc_b, data, rMin, true, clsTarget);
if (fabs(clsMin.first-clsTarget) <= 2*clsMin.second) break;
rMinBound = rMin;
}
while (clsMax.second == 0 || fabs(rMax-limit) > std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
rMax = 0.5*(rMax+limit);
clsMax = eval(w, mc_s, mc_b, data, rMax, true, clsTarget);
if (fabs(clsMax.first-clsTarget) <= 2*clsMax.second) break;
rMaxBound = rMax;
}
expoFit.SetRange(rMinBound,rMaxBound);
break;
}
} while (true);
}
if (!done) { // didn't reach accuracy with scan, now do fit
double rMinBound, rMaxBound; expoFit.GetRange(rMinBound, rMaxBound);
if (verbose) {
std::cout << "\n -- HybridNew, before fit -- \n";
std::cout << "Limit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " [" << rMinBound << ", " << rMaxBound << "]\n";
}
expoFit.FixParameter(0,clsTarget);
double clsmaxfirst = clsMax.first;
if ( clsmaxfirst == 0.0 ) clsmaxfirst = 0.005;
double par1guess = log(clsmaxfirst/clsMin.first)/(rMax-rMin);
expoFit.SetParameter(1,par1guess);
expoFit.SetParameter(2,limit);
limitErr = std::max(fabs(rMinBound-limit), fabs(rMaxBound-limit));
int npoints = 0;
for (int j = 0; j < limitPlot_->GetN(); ++j) {
if (limitPlot_->GetX()[j] >= rMinBound && limitPlot_->GetX()[j] <= rMaxBound) npoints++;
}
for (int i = 0, imax = (readHybridResults_ ? 0 : 8); i <= imax; ++i, ++npoints) {
limitPlot_->Sort();
limitPlot_->Fit(&expoFit,(verbose <= 1 ? "QNR EX0" : "NR EX0"));
if (verbose) {
std::cout << "Fit to " << npoints << " points: " << expoFit.GetParameter(2) << " +/- " << expoFit.GetParError(2) << std::endl;
}
if ((rMinBound < expoFit.GetParameter(2)) && (expoFit.GetParameter(2) < rMaxBound) && (expoFit.GetParError(2) < 0.5*(rMaxBound-rMinBound))) {
// sanity check fit result
limit = expoFit.GetParameter(2);
limitErr = expoFit.GetParError(2);
if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) break;
}
// add one point in the interval.
double rTry = RooRandom::uniform()*(rMaxBound-rMinBound)+rMinBound;
if (i != imax) eval(w, mc_s, mc_b, data, rTry, true, clsTarget);
}
}
if (!plot_.empty() && limitPlot_.get()) {
TCanvas *c1 = new TCanvas("c1","c1");
limitPlot_->Sort();
limitPlot_->SetLineWidth(2);
double xmin = r->getMin(), xmax = r->getMax();
for (int j = 0; j < limitPlot_->GetN(); ++j) {
if (limitPlot_->GetY()[j] > 1.4*clsTarget || limitPlot_->GetY()[j] < 0.6*clsTarget) continue;
xmin = std::min(limitPlot_->GetX()[j], xmin);
xmax = std::max(limitPlot_->GetX()[j], xmax);
}
limitPlot_->GetXaxis()->SetRangeUser(xmin,xmax);
limitPlot_->GetYaxis()->SetRangeUser(0.5*clsTarget, 1.5*clsTarget);
limitPlot_->Draw("AP");
expoFit.Draw("SAME");
TLine line(limitPlot_->GetX()[0], clsTarget, limitPlot_->GetX()[limitPlot_->GetN()-1], clsTarget);
line.SetLineColor(kRed); line.SetLineWidth(2); line.Draw();
line.DrawLine(limit, 0, limit, limitPlot_->GetY()[0]);
line.SetLineWidth(1); line.SetLineStyle(2);
line.DrawLine(limit-limitErr, 0, limit-limitErr, limitPlot_->GetY()[0]);
line.DrawLine(limit+limitErr, 0, limit+limitErr, limitPlot_->GetY()[0]);
c1->Print(plot_.c_str());
}
std::cout << "\n -- Hybrid New -- \n";
std::cout << "Limit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " @ " << cl * 100 << "% CL\n";
if (verbose > 1) std::cout << "Total toys: " << perf_totalToysRun_ << std::endl;
return true;
}
| bool HybridNew::runSignificance | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| double & | limit, | ||
| double & | limitErr, | ||
| const double * | hint | ||
| ) | [virtual] |
Definition at line 267 of file HybridNew.cc.
References applyExpectedQuantile(), dtNoiseDBValidation_cfg::cerr, gather_cfg::cout, create(), EPS, evalGeneric(), expectedFromGrid_, i, iterations_, mass_, max(), name(), readHybridResults_, readToysFromFile(), reportPVal_, rValues_, saveHybridResult_, HcalObjRepresent::setup(), and writeToysHere.
Referenced by run().
{
HybridNew::Setup setup;
std::auto_ptr<RooStats::HybridCalculator> hc(create(w, mc_s, mc_b, data, rValues_, setup));
std::auto_ptr<HypoTestResult> hcResult;
if (readHybridResults_) {
hcResult.reset(readToysFromFile(rValues_));
} else {
hcResult.reset(evalGeneric(*hc));
for (unsigned int i = 1; i < iterations_; ++i) {
std::auto_ptr<HypoTestResult> more(evalGeneric(*hc));
hcResult->Append(more.get());
if (verbose) std::cout << "\t1 - CLb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << std::endl;
}
}
if (hcResult.get() == 0) {
std::cerr << "Hypotest failed" << std::endl;
return false;
}
if (saveHybridResult_) {
TString name = TString::Format("HypoTestResult_mh%g",mass_);
RooLinkedListIter it = rValues_.iterator();
for (RooRealVar *rIn = (RooRealVar*) it.Next(); rIn != 0; rIn = (RooRealVar*) it.Next()) {
name += Form("_%s%g", rIn->GetName(), rIn->getVal());
}
name += Form("_%u", RooRandom::integer(std::numeric_limits<UInt_t>::max() - 1));
writeToysHere->WriteTObject(new HypoTestResult(*hcResult), name);
if (verbose) std::cout << "Hybrid result saved as " << name << " in " << writeToysHere->GetFile()->GetName() << " : " << writeToysHere->GetPath() << std::endl;
}
if (verbose > 1) {
std::cout << "Observed test statistics in data: " << hcResult->GetTestStatisticData() << std::endl;
std::cout << "Background-only toys sampled: " << hcResult->GetNullDistribution()->GetSize() << std::endl;
}
if (expectedFromGrid_) applyExpectedQuantile(*hcResult);
// I don't need to flip the P-values for significances, only for limits
hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()+EPS); // issue with < vs <= in discrete models
double sig = hcResult->Significance();
double sigHi = RooStats::PValueToSignificance( (hcResult->CLb() - hcResult->CLbError()) ) - sig;
double sigLo = RooStats::PValueToSignificance( (hcResult->CLb() + hcResult->CLbError()) ) - sig;
if (reportPVal_) {
limit = hcResult->NullPValue();
limitErr = hcResult->NullPValueError();
} else {
limit = sig;
limitErr = std::max(sigHi,-sigLo);
}
std::cout << "\n -- Hybrid New -- \n";
std::cout << "Significance: " << sig << " " << sigLo << "/+" << sigHi << "\n";
std::cout << "Null p-value: " << hcResult->NullPValue() << " +/- " << hcResult->NullPValueError() << "\n";
return isfinite(limit);
}
| bool HybridNew::runSinglePoint | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| double & | limit, | ||
| double & | limitErr, | ||
| const double * | hint | ||
| ) | [virtual] |
Definition at line 546 of file HybridNew.cc.
References CLs_, clsAccuracy_, gather_cfg::cout, eval(), if(), perf_totalToysRun_, query::result, and rValues_.
Referenced by run().
{
std::pair<double, double> result = eval(w, mc_s, mc_b, data, rValues_, clsAccuracy_ != 0);
std::cout << "\n -- Hybrid New -- \n";
std::cout << (CLs_ ? "CLs = " : "CLsplusb = ") << result.first << " +/- " << result.second << std::endl;
if (verbose > 1) std::cout << "Total toys: " << perf_totalToysRun_ << std::endl;
limit = result.first;
limitErr = result.second;
return true;
}
| bool HybridNew::runTestStatistics | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| double & | limit, | ||
| double & | limitErr, | ||
| const double * | hint | ||
| ) | [virtual] |
Probably useless, but should double-check before deleting this.
Definition at line 556 of file HybridNew.cc.
References applyExpectedQuantile(), gather_cfg::cout, create(), expectedFromGrid_, MessageLogger_cff::limit, HybridNew::Setup::modelConfig_bonly, HybridNew::Setup::qvar, readHybridResults_, readToysFromFile(), query::result, rValues_, HcalObjRepresent::setup(), and testStat_.
Referenced by run().
{
bool isProfile = (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile");
if (readHybridResults_ && expectedFromGrid_) {
std::auto_ptr<RooStats::HypoTestResult> result(readToysFromFile(rValues_));
applyExpectedQuantile(*result);
limit = -2 * result->GetTestStatisticData();
} else {
HybridNew::Setup setup;
std::auto_ptr<RooStats::HybridCalculator> hc(create(w, mc_s, mc_b, data, rValues_, setup));
RooArgSet nullPOI(*setup.modelConfig_bonly.GetSnapshot());
if (isProfile) {
nullPOI.assignValueOnly(rValues_);
}
limit = -2 * setup.qvar->Evaluate(data, nullPOI);
}
if (isProfile) limit = -limit; // there's a sign flip for these two
std::cout << "\n -- Hybrid New -- \n";
std::cout << "-2 ln Q_{"<< testStat_<<"} = " << limit << std::endl;
return true;
}
| void HybridNew::setupPOI | ( | RooStats::ModelConfig * | mc_s | ) | [private] |
Definition at line 220 of file HybridNew.cc.
References NULL, rValue_, and rValues_.
Referenced by run().
{
RooArgSet POI(*mc_s->GetParametersOfInterest());
if (rValue_.find("=") == std::string::npos) {
if (POI.getSize() != 1) throw std::invalid_argument("Error: the argument to --singlePoint or --signalForSignificance is a single value, but there are multiple POIs");
POI.snapshot(rValues_);
errno = 0; // check for errors in str->float conversion
((RooRealVar*)rValues_.first())->setVal(strtod(rValue_.c_str(),NULL));
if (errno != 0) std::invalid_argument("Error: the argument to --singlePoint or --signalForSignificance is not a valid number.");
} else {
std::string::size_type eqidx = 0, colidx = 0, colidx2;
do {
eqidx = rValue_.find("=", colidx);
colidx2 = rValue_.find(",", colidx+1);
if (eqidx == std::string::npos || (colidx2 != std::string::npos && colidx2 < eqidx)) {
throw std::invalid_argument("Error: the argument to --singlePoint or --signalForSignificance is not in the forms 'value' or 'name1=value1,name2=value2,...'\n");
}
std::string poiName = rValue_.substr(colidx, eqidx);
std::string poiVal = rValue_.substr(eqidx+1, (colidx2 == std::string::npos ? std::string::npos : colidx2 - eqidx - 1));
RooAbsArg *poi = POI.find(poiName.c_str());
if (poi == 0) throw std::invalid_argument("Error: unknown parameter '"+poiName+"' passed to --singlePoint or --signalForSignificance.");
rValues_.addClone(*poi);
errno = 0;
rValues_.setRealValue(poi->GetName(), strtod(poiVal.c_str(),NULL));
if (errno != 0) throw std::invalid_argument("Error: invalid value '"+poiVal+"' for parameter '"+poiName+"' passed to --singlePoint or --signalForSignificance.");
} while (colidx2 != std::string::npos);
if (rValues_.getSize() != POI.getSize()) {
throw std::invalid_argument("Error: not all parameters of interest specified in --singlePoint or --signalForSignificance");
}
}
}
| void HybridNew::updateGridData | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| bool | smart, | ||
| double | clsTarget | ||
| ) | [private] |
Definition at line 1379 of file HybridNew.cc.
References gather_cfg::cout, first, grid_, i, max(), n, point, updateGridPoint(), and makeHLTPrescaleTable::values.
Referenced by runLimit().
{
typedef std::map<double, RooStats::HypoTestResult *>::iterator point;
if (!smart) {
for (point it = grid_.begin(), ed = grid_.end(); it != ed; ++it) {
it->second->ResetBit(1);
updateGridPoint(w, mc_s, mc_b, data, it);
}
} else {
typedef std::pair<double,double> CLs_t;
std::vector<point> points; points.reserve(grid_.size());
std::vector<CLs_t> values; values.reserve(grid_.size());
for (point it = grid_.begin(), ed = grid_.end(); it != ed; ++it) { points.push_back(it); values.push_back(CLs_t(-99, -99)); }
int iMin = 0, iMax = points.size()-1;
while (iMax-iMin > 3) {
if (verbose > 1) std::cout << "Bisecting range [" << iMin << ", " << iMax << "]" << std::endl;
int iMid = (iMin+iMax)/2;
CLs_t clsMid = values[iMid] = updateGridPoint(w, mc_s, mc_b, data, points[iMid]);
if (verbose > 1) std::cout << " Midpoint " << iMid << " value " << clsMid.first << " +/- " << clsMid.second << std::endl;
if (clsMid.first - 3*max(clsMid.second,0.01) > clsTarget_) {
if (verbose > 1) std::cout << " Replacing Min" << std::endl;
iMin = iMid; continue;
} else if (clsMid.first + 3*max(clsMid.second,0.01) < clsTarget_) {
if (verbose > 1) std::cout << " Replacing Max" << std::endl;
iMax = iMid; continue;
} else {
if (verbose > 1) std::cout << " Tightening Range" << std::endl;
while (iMin < iMid-1) {
int iLo = (iMin+iMid)/2;
CLs_t clsLo = values[iLo] = updateGridPoint(w, mc_s, mc_b, data, points[iLo]);
if (verbose > 1) std::cout << " Lowpoint " << iLo << " value " << clsLo.first << " +/- " << clsLo.second << std::endl;
if (clsLo.first - 3*max(clsLo.second,0.01) > clsTarget_) iMin = iLo;
else break;
}
while (iMax > iMid+1) {
int iHi = (iMax+iMid)/2;
CLs_t clsHi = values[iHi] = updateGridPoint(w, mc_s, mc_b, data, points[iHi]);
if (verbose > 1) std::cout << " Highpoint " << iHi << " value " << clsHi.first << " +/- " << clsHi.second << std::endl;
if (clsHi.first + 3*max(clsHi.second,0.01) < clsTarget_) iMax = iHi;
else break;
}
break;
}
}
if (verbose > 1) std::cout << "Final range [" << iMin << ", " << iMax << "]" << std::endl;
for (int i = 0; i < iMin; ++i) {
points[i]->second->SetBit(1);
if (verbose > 1) std::cout << " Will not use point " << i << " (r " << points[i]->first << ")" << std::endl;
}
for (int i = iMin; i <= iMax; ++i) {
points[i]->second->ResetBit(1);
if (values[i].first < -2) {
if (verbose > 1) std::cout << " Updaing point " << i << " (r " << points[i]->first << ")" << std::endl;
updateGridPoint(w, mc_s, mc_b, data, points[i]);
}
else if (verbose > 1) std::cout << " Point " << i << " (r " << points[i]->first << ") was already updated during search." << std::endl;
}
for (int i = iMax+1, n = points.size(); i < n; ++i) {
points[i]->second->SetBit(1);
if (verbose > 1) std::cout << " Will not use point " << i << " (r " << points[i]->first << ")" << std::endl;
}
}
}
| void HybridNew::updateGridDataFC | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| bool | smart, | ||
| double | clsTarget | ||
| ) | [private] |
Definition at line 1441 of file HybridNew.cc.
References gather_cfg::cout, create(), EPS, grid_, i, minimizerTolerance_, HybridNew::Setup::modelConfig_bonly, n, point, HybridNew::Setup::qvar, and HcalObjRepresent::setup().
Referenced by runLimit().
{
typedef std::map<double, RooStats::HypoTestResult *>::iterator point;
std::vector<Double_t> rToUpdate; std::vector<point> pointToUpdate;
for (point it = grid_.begin(), ed = grid_.end(); it != ed; ++it) {
it->second->ResetBit(1);
if (it->first == 0 || fabs(it->second->GetTestStatisticData()) <= 2*minimizerTolerance_) {
rToUpdate.push_back(it->first);
pointToUpdate.push_back(it);
}
}
if (verbose > 0) std::cout << "A total of " << rToUpdate.size() << " points will be updated." << std::endl;
Setup setup;
std::auto_ptr<RooStats::HybridCalculator> hc = create(w, mc_s, mc_b, data, rToUpdate.back(), setup);
RooArgSet nullPOI(*setup.modelConfig_bonly.GetSnapshot());
std::vector<Double_t> qVals = ((ProfiledLikelihoodTestStatOpt&)(*setup.qvar)).Evaluate(data, nullPOI, rToUpdate);
for (int i = 0, n = rToUpdate.size(); i < n; ++i) {
pointToUpdate[i]->second->SetTestStatisticData(qVals[i] - EPS);
}
if (verbose > 0) std::cout << "All points have been updated." << std::endl;
}
| std::pair< double, double > HybridNew::updateGridPoint | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| std::map< double, RooStats::HypoTestResult * >::iterator | point | ||
| ) | [private] |
Definition at line 1462 of file HybridNew.cc.
References applyExpectedQuantile(), CLs_, gather_cfg::cout, create(), EPS, eval(), expectedFromGrid_, HybridNew::Setup::modelConfig_bonly, HybridNew::Setup::qvar, alignCSCRings::r, HcalObjRepresent::setup(), and testStat_.
Referenced by updateGridData().
{
typedef std::pair<double,double> CLs_t;
bool isProfile = (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile");
if (point->first == 0 && CLs_) return std::pair<double,double>(1,0);
RooArgSet poi(*mc_s->GetParametersOfInterest());
RooRealVar *r = dynamic_cast<RooRealVar *>(poi.first());
if (expectedFromGrid_) {
applyExpectedQuantile(*point->second);
point->second->SetTestStatisticData(point->second->GetTestStatisticData() + (isProfile ? EPS : EPS));
} else {
Setup setup;
std::auto_ptr<RooStats::HybridCalculator> hc = create(w, mc_s, mc_b, data, point->first, setup);
RooArgSet nullPOI(*setup.modelConfig_bonly.GetSnapshot());
if (isProfile) nullPOI.setRealValue(r->GetName(), point->first);
double testStat = setup.qvar->Evaluate(data, nullPOI);
point->second->SetTestStatisticData(testStat + (isProfile ? EPS : EPS));
}
if (verbose > 1) {
std::cout << "At " << r->GetName() << " = " << point->first << ":\n" <<
"\tCLs = " << point->second->CLs() << " +/- " << point->second->CLsError() << "\n" <<
"\tCLb = " << point->second->CLb() << " +/- " << point->second->CLbError() << "\n" <<
"\tCLsplusb = " << point->second->CLsplusb() << " +/- " << point->second->CLsplusbError() << "\n" <<
std::endl;
}
return eval(*point->second, point->first);
}
| void HybridNew::useGrid | ( | ) | [private] |
Definition at line 1489 of file HybridNew.cc.
References CLs_, eval(), grid_, i, limitPlot_, and n.
Referenced by runLimit().
{
typedef std::pair<double,double> CLs_t;
int i = 0, n = 0;
limitPlot_.reset(new TGraphErrors(1));
for (std::map<double, RooStats::HypoTestResult *>::iterator itg = grid_.begin(), edg = grid_.end(); itg != edg; ++itg, ++i) {
if (itg->second->TestBit(1)) continue;
CLs_t val(1,0);
if (CLs_) {
if (itg->first > 0) val = eval(*itg->second, itg->first);
} else {
val = eval(*itg->second, itg->first);
}
if (val.first == -1) continue;
if (val.second == 0 && (val.first != 1 && val.first != 0)) continue;
limitPlot_->Set(n+1);
limitPlot_->SetPoint( n, itg->first, val.first);
limitPlot_->SetPointError(n, 0, val.second);
n++;
}
}
| void HybridNew::validateOptions | ( | ) | [private] |
Definition at line 185 of file HybridNew.cc.
References CLs_, gather_cfg::cout, fitNuisances_, fork_, MakeSignificance, MakeSignificanceTestStatistics, MakeTestStatistics, nToys_, readHybridResults_, reportPVal_, rule_, testStat_, and workingMode_.
Referenced by applyDefaultOptions(), and applyOptions().
{
if (fork_ > 1) nToys_ /= fork_; // makes more sense
if (rule_ == "CLs") {
CLs_ = true;
} else if (rule_ == "CLsplusb") {
CLs_ = false;
} else if (rule_ == "FC") {
CLs_ = false;
} else {
throw std::invalid_argument("HybridNew: Rule must be either 'CLs' or 'CLsplusb'");
}
if (testStat_ == "SimpleLikelihoodRatio" || testStat_ == "SLR" ) { testStat_ = "LEP"; }
if (testStat_ == "RatioOfProfiledLikelihoods" || testStat_ == "ROPL") { testStat_ = "TEV"; }
if (testStat_ == "ProfileLikelihood" || testStat_ == "PL") { testStat_ = "Profile"; }
if (testStat_ == "ModifiedProfileLikelihood" || testStat_ == "MPL") { testStat_ = "LHC"; }
if (testStat_ == "SignFlipProfileLikelihood" || testStat_ == "SFPL") { testStat_ = "LHCFC"; }
if (testStat_ == "Atlas") { testStat_ = "LHC"; std::cout << "Note: the Atlas test statistics is now known as LHC test statistics.\n" << std::endl; }
if (testStat_ != "LEP" && testStat_ != "TEV" && testStat_ != "LHC" && testStat_ != "LHCFC" && testStat_ != "Profile" && testStat_ != "MLZ") {
throw std::invalid_argument("HybridNew: Test statistics should be one of 'LEP' or 'TEV' or 'LHC' (previously known as 'Atlas') or 'Profile'");
}
if (verbose) {
if (testStat_ == "LEP") std::cout << ">>> using the Simple Likelihood Ratio test statistics (Q_LEP)" << std::endl;
if (testStat_ == "TEV") std::cout << ">>> using the Ratio of Profiled Likelihoods test statistics (Q_TEV)" << std::endl;
if (testStat_ == "LHC") std::cout << ">>> using the Profile Likelihood test statistics modified for upper limits (Q_LHC)" << std::endl;
if (testStat_ == "LHCFC") std::cout << ">>> using the Profile Likelihood test statistics modified for upper limits and Feldman-Cousins (Q_LHCFC)" << std::endl;
if (testStat_ == "Profile") std::cout << ">>> using the Profile Likelihood test statistics not modified for upper limits (Q_Profile)" << std::endl;
if (testStat_ == "MLZ") std::cout << ">>> using the Maximum likelihood estimator of the signal strength as test statistics" << std::endl;
}
if (readHybridResults_ || workingMode_ == MakeTestStatistics || workingMode_ == MakeSignificanceTestStatistics) {
// If not generating toys, don't need to fit nuisance parameters
fitNuisances_ = false;
}
if (reportPVal_ && workingMode_ != MakeSignificance) throw std::invalid_argument("HybridNew: option --pvalue must go together with --significance");
}
std::string HybridNew::algo_ = "logSecant" [static, private] |
Definition at line 60 of file HybridNew.h.
Referenced by HybridNew(), and runLimit().
bool HybridNew::CLs_ = false [static, private] |
Definition at line 43 of file HybridNew.h.
Referenced by applyClsQuantile(), create(), eval(), runLimit(), runSinglePoint(), updateGridPoint(), useGrid(), and validateOptions().
double HybridNew::clsAccuracy_ = 0.005 [static, private] |
Definition at line 42 of file HybridNew.h.
Referenced by eval(), HybridNew(), and runSinglePoint().
bool HybridNew::clsQuantiles_ = true [private] |
Definition at line 52 of file HybridNew.h.
Referenced by applyExpectedQuantile(), and HybridNew().
bool HybridNew::expectedFromGrid_ = false [static, private] |
Definition at line 52 of file HybridNew.h.
Referenced by applyExpectedQuantile(), applyOptions(), create(), eval(), runSignificance(), runTestStatistics(), and updateGridPoint().
bool HybridNew::fitNuisances_ = false [private] |
Definition at line 46 of file HybridNew.h.
Referenced by applyOptions(), create(), HybridNew(), and validateOptions().
unsigned int HybridNew::fork_ = 1 [private] |
Definition at line 58 of file HybridNew.h.
Referenced by evalGeneric(), evalWithFork(), HybridNew(), and validateOptions().
bool HybridNew::fullBToys_ = false [static, private] |
Definition at line 54 of file HybridNew.h.
Referenced by applyOptions(), and create().
bool HybridNew::fullGrid_ = false [static, private] |
Definition at line 55 of file HybridNew.h.
Referenced by applyOptions(), and runLimit().
bool HybridNew::genGlobalObs_ = false [private] |
Definition at line 46 of file HybridNew.h.
Referenced by applyOptions(), create(), and HybridNew().
bool HybridNew::genNuisances_ = true [static, private] |
Definition at line 46 of file HybridNew.h.
Referenced by applyOptions(), create(), and HybridNew().
std::map<double, RooStats::HypoTestResult *> HybridNew::grid_ [private] |
Definition at line 119 of file HybridNew.h.
Referenced by clearGrid(), readGrid(), runLimit(), updateGridData(), updateGridDataFC(), and useGrid().
std::string HybridNew::gridFile_ = "" [static, private] |
Definition at line 51 of file HybridNew.h.
Referenced by HybridNew(), and runLimit().
bool HybridNew::importanceSamplingAlt_ = false [private] |
Definition at line 59 of file HybridNew.h.
Referenced by create().
bool HybridNew::importanceSamplingNull_ = false [static, private] |
Definition at line 59 of file HybridNew.h.
Referenced by create().
double HybridNew::interpAccuracy_ = 0.2 [private] |
Definition at line 42 of file HybridNew.h.
Referenced by HybridNew(), and runLimit().
unsigned int HybridNew::iterations_ = 1 [static, private] |
Definition at line 49 of file HybridNew.h.
Referenced by eval(), HybridNew(), and runSignificance().
std::auto_ptr<TGraphErrors> HybridNew::limitPlot_ [private] |
Definition at line 73 of file HybridNew.h.
Referenced by eval(), runLimit(), and useGrid().
float HybridNew::mass_ [private] |
Definition at line 70 of file HybridNew.h.
Referenced by applyOptions(), eval(), readGrid(), readToysFromFile(), and runSignificance().
std::string HybridNew::minimizerAlgo_ = "Minuit2" [static, private] |
Definition at line 62 of file HybridNew.h.
Referenced by HybridNew(), and run().
float HybridNew::minimizerTolerance_ = 1e-2 [static, private] |
Definition at line 63 of file HybridNew.h.
Referenced by eval(), HybridNew(), run(), and updateGridDataFC().
unsigned int HybridNew::nCpu_ = 0 [static, private] |
Definition at line 58 of file HybridNew.h.
Referenced by create(), and HybridNew().
bool HybridNew::newToyMCSampler_ = true [static, private] |
Definition at line 67 of file HybridNew.h.
Referenced by create(), and HybridNew().
bool HybridNew::noUpdateGrid_ = false [static, private] |
Definition at line 57 of file HybridNew.h.
Referenced by applyOptions(), eval(), and runLimit().
unsigned int HybridNew::nToys_ = 500 [static, private] |
Definition at line 41 of file HybridNew.h.
Referenced by create(), eval(), HybridNew(), and validateOptions().
bool HybridNew::optimizeProductPdf_ = true [static, private] |
Definition at line 65 of file HybridNew.h.
Referenced by create(), and HybridNew().
bool HybridNew::optimizeTestStatistics_ = true [static, private] |
Definition at line 66 of file HybridNew.h.
Referenced by create(), and HybridNew().
unsigned int HybridNew::perf_totalToysRun_ [private] |
Definition at line 76 of file HybridNew.h.
Referenced by eval(), run(), runLimit(), and runSinglePoint().
std::string HybridNew::plot_ [static, private] |
Definition at line 61 of file HybridNew.h.
Referenced by eval(), HybridNew(), readToysFromFile(), and runLimit().
float HybridNew::quantileForExpectedFromGrid_ = 0.5 [static, private] |
Definition at line 53 of file HybridNew.h.
Referenced by applyClsQuantile(), applyExpectedQuantile(), applyOptions(), applySignalQuantile(), create(), eval(), and HybridNew().
double HybridNew::rAbsAccuracy_ = 0.1 [private] |
Definition at line 42 of file HybridNew.h.
Referenced by HybridNew(), and runLimit().
bool HybridNew::readHybridResults_ = false [private] |
Definition at line 50 of file HybridNew.h.
Referenced by applyOptions(), eval(), runLimit(), runSignificance(), runTestStatistics(), and validateOptions().
bool HybridNew::reportPVal_ = false [static, private] |
Definition at line 45 of file HybridNew.h.
Referenced by applyOptions(), runSignificance(), and validateOptions().
bool HybridNew::rMaxSet_ = false [static, private] |
Definition at line 69 of file HybridNew.h.
Referenced by applyOptions(), and runLimit().
bool HybridNew::rMinSet_ = false [static, private] |
Definition at line 68 of file HybridNew.h.
Referenced by applyOptions(), and runLimit().
double HybridNew::rRelAccuracy_ = 0.05 [private] |
Definition at line 42 of file HybridNew.h.
Referenced by HybridNew(), and runLimit().
std::string HybridNew::rule_ = "CLs" [static, private] |
Definition at line 44 of file HybridNew.h.
Referenced by eval(), HybridNew(), runLimit(), and validateOptions().
std::string HybridNew::rValue_ = "1.0" [static, private] |
Definition at line 47 of file HybridNew.h.
Referenced by applyOptions(), HybridNew(), and setupPOI().
RooArgSet HybridNew::rValues_ [static, private] |
Definition at line 48 of file HybridNew.h.
Referenced by readGrid(), run(), runSignificance(), runSinglePoint(), runTestStatistics(), and setupPOI().
bool HybridNew::saveGrid_ = false [static, private] |
Definition at line 56 of file HybridNew.h.
Referenced by applyOptions(), and runLimit().
bool HybridNew::saveHybridResult_ = false [static, private] |
Definition at line 50 of file HybridNew.h.
Referenced by applyOptions(), eval(), and runSignificance().
std::string HybridNew::testStat_ = "LEP" [private] |
Definition at line 44 of file HybridNew.h.
Referenced by applyOptions(), create(), eval(), HybridNew(), runLimit(), runTestStatistics(), updateGridPoint(), and validateOptions().
HybridNew::WorkingMode HybridNew::workingMode_ = MakeLimit [static, private] |
Definition at line 40 of file HybridNew.h.
Referenced by applyDefaultOptions(), applyExpectedQuantile(), applyOptions(), create(), eval(), readToysFromFile(), run(), and validateOptions().
1.7.3