CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Classes | Public Types | Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes
HybridNew Class Reference

#include <HybridNew.h>

Inheritance diagram for HybridNew:
LimitAlgo

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)
 
- Public Member Functions inherited from LimitAlgo
 LimitAlgo ()
 
 LimitAlgo (const char *desc)
 
const
boost::program_options::options_description & 
options () const
 
virtual void setNToys (const int)
 
virtual void setToyNumber (const int)
 
virtual ~LimitAlgo ()
 

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 More...
 
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) More...
 
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 More...
 
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) More...
 
std::pair< double, double > eval (RooStats::HybridCalculator &hc, const RooAbsCollection &rVals, bool adaptive=false, double clsTarget=-1)
 
std::pair< double, double > eval (const RooStats::HypoTestResult &hcres, const RooAbsCollection &rVals)
 
std::pair< double, double > eval (const RooStats::HypoTestResult &hcres, double rVal)
 
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

std::map< double,
RooStats::HypoTestResult * > 
grid_
 
std::auto_ptr< TGraphErrors > limitPlot_
 
float mass_
 
unsigned int perf_totalToysRun_
 

Static Private Attributes

static std::string algo_ = "logSecant"
 
static bool CLs_ = false
 
static double clsAccuracy_ = 0.005
 
static bool clsQuantiles_ = true
 
static bool expectedFromGrid_ = false
 
static bool fitNuisances_ = false
 
static unsigned int fork_ = 1
 
static bool fullBToys_ = false
 
static bool fullGrid_ = false
 
static bool genGlobalObs_ = false
 
static bool genNuisances_ = true
 
static std::string gridFile_ = ""
 
static bool importanceSamplingAlt_ = false
 
static bool importanceSamplingNull_ = false
 
static double interpAccuracy_ = 0.2
 
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 double rAbsAccuracy_ = 0.1
 
static bool readHybridResults_ = false
 
static bool reportPVal_ = false
 
static bool rMaxSet_ = false
 
static bool rMinSet_ = false
 
static double rRelAccuracy_ = 0.05
 
static std::string rule_ = "CLs"
 
static std::string rValue_ = "1.0"
 
static RooArgSet rValues_
 
static bool saveGrid_ = false
 
static bool saveHybridResult_ = false
 
static std::string testStat_ = "LEP"
 
static WorkingMode workingMode_ = MakeLimit
 

Additional Inherited Members

- Protected Attributes inherited from LimitAlgo
boost::program_options::options_description options_
 

Detailed Description

Module to compute limits by tossing toys (CLs, CLs+b, Feldman-Cousins), and related significances

Author
Luca Lista (INFN), Giovanni Petrucciani (UCSD)

Definition at line 21 of file HybridNew.h.

Member Enumeration Documentation

Enumerator
MakeLimit 
MakeSignificance 
MakePValues 
MakeTestStatistics 
MakeSignificanceTestStatistics 

Definition at line 38 of file HybridNew.h.

Constructor & Destructor Documentation

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_.

89  :
90 LimitAlgo("HybridNew specific options") {
91  options_.add_options()
92  ("rule", boost::program_options::value<std::string>(&rule_)->default_value(rule_), "Rule to use: CLs, CLsplusb")
93  ("testStat",boost::program_options::value<std::string>(&testStat_)->default_value(testStat_), "Test statistics: LEP, TEV, LHC (previously known as Atlas), Profile.")
94  ("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,...'")
95  ("onlyTestStat", "Just compute test statistics, not actual p-values (works only with --singlePoint)")
96  ("generateNuisances", boost::program_options::value<bool>(&genNuisances_)->default_value(genNuisances_), "Generate nuisance parameters for each toy")
97  ("generateExternalMeasurements", boost::program_options::value<bool>(&genGlobalObs_)->default_value(genGlobalObs_), "Generate external measurements for each toy, taken from the GlobalObservables of the ModelConfig")
98  ("fitNuisances", boost::program_options::value<bool>(&fitNuisances_)->default_value(fitNuisances_), "Fit the nuisances, when not generating them.")
99  ("searchAlgo", boost::program_options::value<std::string>(&algo_)->default_value(algo_), "Algorithm to use to search for the limit (bisection, logSecant)")
100  ("toysH,T", boost::program_options::value<unsigned int>(&nToys_)->default_value(nToys_), "Number of Toy MC extractions to compute CLs+b, CLb and CLs")
101  ("clsAcc", boost::program_options::value<double>(&clsAccuracy_ )->default_value(clsAccuracy_), "Absolute accuracy on CLs to reach to terminate the scan")
102  ("rAbsAcc", boost::program_options::value<double>(&rAbsAccuracy_)->default_value(rAbsAccuracy_), "Absolute accuracy on r to reach to terminate the scan")
103  ("rRelAcc", boost::program_options::value<double>(&rRelAccuracy_)->default_value(rRelAccuracy_), "Relative accuracy on r to reach to terminate the scan")
104  ("interpAcc", boost::program_options::value<double>(&interpAccuracy_)->default_value(interpAccuracy_), "Minimum uncertainty from interpolation delta(x)/(max(x)-min(x))")
105  ("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)")
106  ("fork", boost::program_options::value<unsigned int>(&fork_)->default_value(fork_), "Fork to N processes before running the toys (set to 0 for debugging)")
107  ("nCPU", boost::program_options::value<unsigned int>(&nCpu_)->default_value(nCpu_), "Use N CPUs with PROOF Lite (experimental!)")
108  ("saveHybridResult", "Save result in the output file")
109  ("readHybridResults", "Read and merge results from file (requires 'toysFile' or 'grid')")
110  ("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.")
111  ("expectedFromGrid", boost::program_options::value<float>(&quantileForExpectedFromGrid_)->default_value(0.5), "Use the grid to compute the expected limit for this quantile")
112  ("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)")
113  ("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")
114  //("importanceSamplingNull", boost::program_options::value<bool>(&importanceSamplingNull_)->default_value(importanceSamplingNull_),
115  // "Enable importance sampling for null hypothesis (background only)")
116  //("importanceSamplingAlt", boost::program_options::value<bool>(&importanceSamplingAlt_)->default_value(importanceSamplingAlt_),
117  // "Enable importance sampling for alternative hypothesis (signal plus background)")
118  ("optimizeTestStatistics", boost::program_options::value<bool>(&optimizeTestStatistics_)->default_value(optimizeTestStatistics_),
119  "Use optimized test statistics if the likelihood is not extended (works for LEP and TEV test statistics).")
120  ("optimizeProductPdf", boost::program_options::value<bool>(&optimizeProductPdf_)->default_value(optimizeProductPdf_),
121  "Optimize the code factorizing pdfs")
122  ("minimizerAlgo", boost::program_options::value<std::string>(&minimizerAlgo_)->default_value(minimizerAlgo_), "Choice of minimizer used for profiling (Minuit vs Minuit2)")
123  ("minimizerTolerance", boost::program_options::value<float>(&minimizerTolerance_)->default_value(minimizerTolerance_), "Tolerance for minimizer used for profiling")
124  ("plot", boost::program_options::value<std::string>(&plot_), "Save a plot of the result (test statistics distributions or limit scan)")
125  ("frequentist", "Shortcut to switch to Frequentist mode (--generateNuisances=0 --generateExternalMeasurements=1 --fitNuisances=1)")
126  ("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.")
127  ("fullGrid", "Evaluate p-values at all grid points, without optimitations")
128  ("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)")
129  ("noUpdateGrid", "Do not update test statistics at grid points")
130  ("fullBToys", "Run as many B toys as S ones (default is to run 1/4 of b-only toys)")
131  ("pvalue", "Report p-value instead of significance (when running with --significance)")
132  ;
133 }
static std::string rule_
Definition: HybridNew.h:44
static unsigned int iterations_
Definition: HybridNew.h:49
static unsigned int fork_
Definition: HybridNew.h:58
static std::string rValue_
Definition: HybridNew.h:47
static std::string algo_
Definition: HybridNew.h:60
static float minimizerTolerance_
Definition: HybridNew.h:63
static double interpAccuracy_
Definition: HybridNew.h:42
static std::string minimizerAlgo_
Definition: HybridNew.h:62
static double rAbsAccuracy_
Definition: HybridNew.h:42
static float quantileForExpectedFromGrid_
Definition: HybridNew.h:53
LimitAlgo()
Definition: LimitAlgo.h:19
static unsigned int nCpu_
Definition: HybridNew.h:58
static unsigned int nToys_
Definition: HybridNew.h:41
static bool genGlobalObs_
Definition: HybridNew.h:46
static std::string testStat_
Definition: HybridNew.h:44
static bool optimizeProductPdf_
Definition: HybridNew.h:65
static std::string gridFile_
Definition: HybridNew.h:51
static bool optimizeTestStatistics_
Definition: HybridNew.h:66
static bool genNuisances_
Definition: HybridNew.h:46
static std::string plot_
Definition: HybridNew.h:61
static double rRelAccuracy_
Definition: HybridNew.h:42
static bool newToyMCSampler_
Definition: HybridNew.h:67
boost::program_options::options_description options_
Definition: LimitAlgo.h:32
static bool clsQuantiles_
Definition: HybridNew.h:52
static bool fitNuisances_
Definition: HybridNew.h:46
static double clsAccuracy_
Definition: HybridNew.h:42

Member Function Documentation

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().

1084  {
1085  RooStats::SamplingDistribution * bDistribution = hcres.GetNullDistribution(), * sDistribution = hcres.GetAltDistribution();
1086  const std::vector<Double_t> & bdist = bDistribution->GetSamplingDistribution();
1087  const std::vector<Double_t> & bweight = bDistribution->GetSampleWeights();
1088  const std::vector<Double_t> & sdist = sDistribution->GetSamplingDistribution();
1089  const std::vector<Double_t> & sweight = sDistribution->GetSampleWeights();
1090  TStopwatch timer;
1091 
1093  timer.Start();
1094  std::vector<std::pair<double,double> > bcumul; bcumul.reserve(bdist.size());
1095  std::vector<std::pair<double,double> > scumul; scumul.reserve(sdist.size());
1096  double btot = 0, stot = 0;
1097  for (std::vector<Double_t>::const_iterator it = bdist.begin(), ed = bdist.end(), itw = bweight.begin(); it != ed; ++it, ++itw) {
1098  bcumul.push_back(std::pair<double,double>(*it, *itw));
1099  btot += *itw;
1100  }
1101  for (std::vector<Double_t>::const_iterator it = sdist.begin(), ed = sdist.end(), itw = sweight.begin(); it != ed; ++it, ++itw) {
1102  scumul.push_back(std::pair<double,double>(*it, *itw));
1103  stot += *itw;
1104  }
1105  double sinv = 1.0/stot, binv = 1.0/btot, runningSum;
1106  // now compute integral distribution of Q(s+b data) so that we can quickly compute the CL_{s+b} for all test stats.
1107  std::sort(scumul.begin(), scumul.end());
1108  runningSum = 0;
1109  for (std::vector<std::pair<double,double> >::reverse_iterator it = scumul.rbegin(), ed = scumul.rend(); it != ed; ++it) {
1110  runningSum += it->second;
1111  it->second = runningSum * sinv;
1112  }
1113  std::sort(bcumul.begin(), bcumul.end());
1114  std::vector<std::pair<double,std::pair<double,double> > > xcumul; xcumul.reserve(bdist.size());
1115  runningSum = 0;
1116  std::vector<std::pair<double,double> >::const_iterator sbegin = scumul.begin(), send = scumul.end();
1117  //int k = 0;
1118  for (std::vector<std::pair<double,double> >::const_reverse_iterator it = bcumul.rbegin(), ed = bcumul.rend(); it != ed; ++it) {
1119  runningSum += it->second;
1120  std::vector<std::pair<double,double> >::const_iterator match = std::upper_bound(sbegin, send, std::pair<double,double>(it->first, 0));
1121  if (match == send) {
1122  //std::cout << "Did not find match, for it->first == " << it->first << ", as back = ( " << scumul.back().first << " , " << scumul.back().second << " ) " << std::endl;
1123  double clsb = (scumul.back().second > 0.5 ? 1.0 : 0.0), clb = runningSum*binv, cls = clsb / clb;
1124  xcumul.push_back(std::make_pair(CLs_ ? cls : clsb, *it));
1125  } else {
1126  double clsb = match->second, clb = runningSum*binv, cls = clsb / clb;
1127  //if ((++k) % 100 == 0) printf("At %+8.5f CLb = %6.4f, CLsplusb = %6.4f, CLs =%7.4f\n", it->first, clb, clsb, cls);
1128  xcumul.push_back(std::make_pair(CLs_ ? cls : clsb, *it));
1129  }
1130  }
1131  // sort
1132  std::sort(xcumul.begin(), xcumul.end());
1133  // get quantile
1134  runningSum = 0; double cut = quantileForExpectedFromGrid_ * btot;
1135  for (std::vector<std::pair<double,std::pair<double,double> > >::const_iterator it = xcumul.begin(), ed = xcumul.end(); it != ed; ++it) {
1136  runningSum += it->second.second;
1137  if (runningSum >= cut) {
1138  hcres.SetTestStatisticData(it->second.first);
1139  //std::cout << "CLs quantile = " << it->first << " for test stat = " << it->second.first << std::endl;
1140  break;
1141  }
1142  }
1143  //std::cout << "CLs quantile = " << (CLs_ ? hcres.CLs() : hcres.CLsplusb()) << std::endl;
1144  //std::cout << "Computed quantiles in " << timer.RealTime() << " s" << std::endl;
1145 #if 0
1146 
1147  timer.Start();
1148  std::vector<std::pair<double, double> > values(bdist.size());
1149  for (int i = 0, n = bdist.size(); i < n; ++i) {
1150  hcres.SetTestStatisticData( bdist[i] );
1151  values[i] = std::pair<double, double>(CLs_ ? hcres.CLs() : hcres.CLsplusb(), bdist[i]);
1152  }
1153  std::sort(values.begin(), values.end());
1154  int index = std::min<int>(floor((1.-quantileForExpectedFromGrid_) * values.size()+0.5), values.size());
1155  std::cout << "CLs quantile = " << values[index].first << " for test stat = " << values[index].second << std::endl;
1156  hcres.SetTestStatisticData(values[index].second);
1157  std::cout << "CLs quantile = " << (CLs_ ? hcres.CLs() : hcres.CLsplusb()) << " for test stat = " << values[index].second << std::endl;
1158  std::cout << "Computed quantiles in " << timer.RealTime() << " s" << std::endl;
1159 #endif
1160 }
int i
Definition: DBlmapReader.cc:9
U second(std::pair< T, U > const &p)
static float quantileForExpectedFromGrid_
Definition: HybridNew.h:53
static bool CLs_
Definition: HybridNew.h:43
tuple cout
Definition: gather_cfg.py:121
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:6
void HybridNew::applyDefaultOptions ( )
virtual

Reimplemented from LimitAlgo.

Definition at line 180 of file HybridNew.cc.

References MakeLimit, validateOptions(), and workingMode_.

180  {
182  validateOptions();
183 }
void validateOptions()
Definition: HybridNew.cc:185
static WorkingMode workingMode_
Definition: HybridNew.h:40
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().

1068  {
1069  if (expectedFromGrid_) {
1071  applySignalQuantile(hcres);
1072  } else if (clsQuantiles_) {
1073  applyClsQuantile(hcres);
1074  } else {
1075  std::vector<Double_t> btoys = hcres.GetNullDistribution()->GetSamplingDistribution();
1076  std::sort(btoys.begin(), btoys.end());
1077  Double_t testStat = btoys[std::min<int>(floor((1.-quantileForExpectedFromGrid_) * btoys.size()+0.5), btoys.size())];
1078  if (verbose > 0) std::cout << "Text statistics for " << quantileForExpectedFromGrid_ << " quantile: " << testStat << std::endl;
1079  hcres.SetTestStatisticData(testStat);
1080  //std::cout << "CLs quantile = " << (CLs_ ? hcres.CLs() : hcres.CLsplusb()) << " for test stat = " << testStat << std::endl;
1081  }
1082  }
1083 }
static WorkingMode workingMode_
Definition: HybridNew.h:40
void applySignalQuantile(RooStats::HypoTestResult &hcres)
Definition: HybridNew.cc:1162
static float quantileForExpectedFromGrid_
Definition: HybridNew.h:53
tuple cout
Definition: gather_cfg.py:121
static bool expectedFromGrid_
Definition: HybridNew.h:52
void applyClsQuantile(RooStats::HypoTestResult &hcres)
Definition: HybridNew.cc:1084
static bool clsQuantiles_
Definition: HybridNew.h:52
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_.

135  {
136  rMinSet_ = vm.count("rMin")>0; rMaxSet_ = vm.count("rMax")>0;
137  if (vm.count("expectedFromGrid") && !vm["expectedFromGrid"].defaulted()) {
138  //if (!vm.count("grid")) throw std::invalid_argument("HybridNew: Can't use --expectedFromGrid without --grid!");
139  if (quantileForExpectedFromGrid_ <= 0 || quantileForExpectedFromGrid_ >= 1.0) throw std::invalid_argument("HybridNew: the quantile for the expected limit must be between 0 and 1");
140  expectedFromGrid_ = true;
142  }
143  if (vm.count("frequentist")) {
145  if (vm["testStat"].defaulted()) testStat_ = "LHC";
146  if (vm["toys"].as<int>() > 0 and vm.count("toysFrequentist")) {
147  if (vm["fitNuisances"].defaulted() && withSystematics) {
148  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;
149  fitNuisances_ = false;
150  }
151  }
152  }
153  if (genGlobalObs_ && genNuisances_) {
154  std::cerr << "ALERT: generating both global observables and nuisance parameters at the same time is not validated." << std::endl;
155  }
156  if (!vm["singlePoint"].defaulted()) {
157  if (doSignificance_) throw std::invalid_argument("HybridNew: Can't use --significance and --singlePoint at the same time");
158  workingMode_ = ( vm.count("onlyTestStat") ? MakeTestStatistics : MakePValues );
159  } else if (vm.count("onlyTestStat")) {
161  else throw std::invalid_argument("HybridNew: --onlyTestStat works only with --singlePoint or --significance");
162  } else if (doSignificance_) {
164  rValue_ = vm["signalForSignificance"].as<std::string>();
165  } else {
167  }
168  saveHybridResult_ = vm.count("saveHybridResult");
169  readHybridResults_ = vm.count("readHybridResults") || vm.count("grid");
170  if (readHybridResults_ && !(vm.count("toysFile") || vm.count("grid"))) throw std::invalid_argument("HybridNew: must have 'toysFile' or 'grid' option to have 'readHybridResults'\n");
171  mass_ = vm["mass"].as<float>();
172  fullGrid_ = vm.count("fullGrid");
173  saveGrid_ = vm.count("saveGrid");
174  fullBToys_ = vm.count("fullBToys");
175  noUpdateGrid_ = vm.count("noUpdateGrid");
176  reportPVal_ = vm.count("pvalue");
177  validateOptions();
178 }
bool doSignificance_
Definition: Combine.cc:69
void validateOptions()
Definition: HybridNew.cc:185
bool withSystematics
Definition: Combine.cc:68
float mass_
Definition: HybridNew.h:70
static std::string rValue_
Definition: HybridNew.h:47
static WorkingMode workingMode_
Definition: HybridNew.h:40
static bool saveHybridResult_
Definition: HybridNew.h:50
static float quantileForExpectedFromGrid_
Definition: HybridNew.h:53
static bool reportPVal_
Definition: HybridNew.h:45
Float_t g_quantileExpected_
Definition: Combine.cc:63
static bool saveGrid_
Definition: HybridNew.h:56
static bool fullGrid_
Definition: HybridNew.h:55
static bool rMaxSet_
Definition: HybridNew.h:69
static bool genGlobalObs_
Definition: HybridNew.h:46
static std::string testStat_
Definition: HybridNew.h:44
static bool noUpdateGrid_
Definition: HybridNew.h:57
static bool genNuisances_
Definition: HybridNew.h:46
static bool rMinSet_
Definition: HybridNew.h:68
tuple cout
Definition: gather_cfg.py:121
static bool expectedFromGrid_
Definition: HybridNew.h:52
static bool fitNuisances_
Definition: HybridNew.h:46
static bool readHybridResults_
Definition: HybridNew.h:50
static bool fullBToys_
Definition: HybridNew.h:54
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().

1162  {
1163  std::vector<Double_t> stoys = hcres.GetAltDistribution()->GetSamplingDistribution();
1164  std::sort(stoys.begin(), stoys.end());
1165  Double_t testStat = stoys[std::min<int>(floor(quantileForExpectedFromGrid_ * stoys.size()+0.5), stoys.size())];
1166  if (verbose > 0) std::cout << "Text statistics for " << quantileForExpectedFromGrid_ << " quantile: " << testStat << std::endl;
1167  hcres.SetTestStatisticData(testStat);
1168 }
static float quantileForExpectedFromGrid_
Definition: HybridNew.h:53
tuple cout
Definition: gather_cfg.py:121
void HybridNew::clearGrid ( )
private

Definition at line 1509 of file HybridNew.cc.

References grid_.

Referenced by readGrid().

1509  {
1510  for (std::map<double, RooStats::HypoTestResult *>::iterator it = grid_.begin(), ed = grid_.end(); it != ed; ++it) {
1511  delete it->second;
1512  }
1513  grid_.clear();
1514 }
std::map< double, RooStats::HypoTestResult * > grid_
Definition: HybridNew.h:119
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().

639  {
640  using namespace RooStats;
641 
642  w->loadSnapshot("clean");
643  // realData_ = &data;
644 
645  RooArgSet poi(*mc_s->GetParametersOfInterest()), params(poi);
646  RooRealVar *r = dynamic_cast<RooRealVar *>(poi.first());
647 
648  if (poi.getSize() == 1) { // here things are a bit more convoluted, although they could probably be cleaned up
649  double rVal = ((RooAbsReal*)rVals.first())->getVal();
650  if (testStat_ != "MLZ") r->setMax(rVal);
651  r->setVal(rVal);
652  if (testStat_ == "LHC" || testStat_ == "Profile") {
653  r->setConstant(false); r->setMin(0);
655  r->setVal(0);
656  r->removeMax();
657  }
658  } else {
659  r->setConstant(true);
660  }
661  } else {
662  if (testStat_ == "Profile") utils::setAllConstant(poi, false);
663  else if (testStat_ == "LEP" || testStat_ == "TEV") utils::setAllConstant(poi, true);
664  poi.assignValueOnly(rVals);
665  }
666 
667  //set value of "globalConstrained" nuisances to the value of the corresponding global observable and set them constant
668  //also fill the RooArgLists which can be passed to the test statistic in order to produce the correct behaviour
669  //this is only supported currently for the optimized version of the LHC-type test statistic
670  RooArgList gobsParams;
671  RooArgList gobs;
673  RooArgList allnuis(*mc_s->GetNuisanceParameters());
674  RooArgList allgobs(*mc_s->GetGlobalObservables());
675  for (int i=0; i<allnuis.getSize(); ++i) {
676  RooRealVar *nuis = (RooRealVar*)allnuis.at(i);
677  if (nuis->getAttribute("globalConstrained")) {
678  RooRealVar *glob = (RooRealVar*)allgobs.find(TString::Format("%s_In",nuis->GetName()));
679  if (glob) {
680  nuis->setVal(glob->getVal());
681  nuis->setConstant();
682  gobsParams.add(*nuis);
683  gobs.add(*glob);
684  }
685  }
686  }
687  }
688 
689  std::auto_ptr<RooFitResult> fitMu, fitZero;
690  if (fitNuisances_ && mc_s->GetNuisanceParameters() && withSystematics) {
691  TStopwatch timer;
692  bool isExt = mc_s->GetPdf()->canBeExtended();
693  utils::setAllConstant(poi, true);
695  RooAbsPdf *pdfB = mc_b->GetPdf();
696  if (poi.getSize() == 1) {
697  r->setVal(0); pdfB = mc_s->GetPdf();
698  }
699  timer.Start();
700  {
701  CloseCoutSentry sentry(verbose < 3);
702  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())));
703  }
704  if (verbose > 1) { std::cout << "Zero signal fit" << std::endl; fitZero->Print("V"); }
705  if (verbose > 1) { std::cout << "Fitting of the background hypothesis done in " << timer.RealTime() << " s" << std::endl; }
706  poi.assignValueOnly(rVals);
707  timer.Start();
708  {
709  CloseCoutSentry sentry(verbose < 3);
710  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())));
711  }
712  if (verbose > 1) { std::cout << "Reference signal fit" << std::endl; fitMu->Print("V"); }
713  if (verbose > 1) { std::cout << "Fitting of the signal-plus-background hypothesis done in " << timer.RealTime() << " s" << std::endl; }
714  } else { fitNuisances_ = false; }
715 
716  // since ModelConfig cannot allow re-setting sets, we have to re-make everything
717  setup.modelConfig = ModelConfig("HybridNew_mc_s", w);
718  setup.modelConfig.SetPdf(*mc_s->GetPdf());
719  setup.modelConfig.SetObservables(*mc_s->GetObservables());
720  setup.modelConfig.SetParametersOfInterest(*mc_s->GetParametersOfInterest());
721  if (withSystematics) {
722  if (genNuisances_ && mc_s->GetNuisanceParameters()) setup.modelConfig.SetNuisanceParameters(*mc_s->GetNuisanceParameters());
723  if (genGlobalObs_ && mc_s->GetGlobalObservables()) setup.modelConfig.SetGlobalObservables(*mc_s->GetGlobalObservables());
724  // if (genGlobalObs_ && mc_s->GetGlobalObservables()) snapGlobalObs_.reset(mc_s->GetGlobalObservables()->snapshot());
725  }
726 
727  setup.modelConfig_bonly = ModelConfig("HybridNew_mc_b", w);
728  setup.modelConfig_bonly.SetPdf(*mc_b->GetPdf());
729  setup.modelConfig_bonly.SetObservables(*mc_b->GetObservables());
730  setup.modelConfig_bonly.SetParametersOfInterest(*mc_b->GetParametersOfInterest());
731  if (withSystematics) {
732  if (genNuisances_ && mc_b->GetNuisanceParameters()) setup.modelConfig_bonly.SetNuisanceParameters(*mc_b->GetNuisanceParameters());
733  if (genGlobalObs_ && mc_b->GetGlobalObservables()) setup.modelConfig_bonly.SetGlobalObservables(*mc_b->GetGlobalObservables());
734  }
735 
736  if (withSystematics && !genNuisances_) {
737  // The pdf will contain non-const parameters which are not observables
738  // and the HybridCalculator will assume they're nuisances and try to generate them
739  // to avoid this, we force him to generate a fake nuisance instead
740  if (w->var("__HybridNew_fake_nuis__") == 0) {
741  w->factory("__HybridNew_fake_nuis__[0.5,0,1]");
742  w->factory("Uniform::__HybridNew_fake_nuisPdf__(__HybridNew_fake_nuis__)");
743  }
744  setup.modelConfig.SetNuisanceParameters(RooArgSet(*w->var("__HybridNew_fake_nuis__")));
745  setup.modelConfig_bonly.SetNuisanceParameters(RooArgSet(*w->var("__HybridNew_fake_nuis__")));
746  }
747 
748 
749  // create snapshots
750  RooArgSet paramsZero;
751  if (poi.getSize() == 1) { // in the trivial 1D case, the background has POI=0.
752  paramsZero.addClone(*rVals.first());
753  paramsZero.setRealValue(rVals.first()->GetName(), 0);
754  }
755  if (fitNuisances_) params.add(fitMu->floatParsFinal());
756  if (fitNuisances_) paramsZero.addClone(fitZero->floatParsFinal());
757  setup.modelConfig.SetSnapshot(params);
758  setup.modelConfig_bonly.SetSnapshot(paramsZero);
759  TString paramsSnapName = TString::Format("%s_%s_snapshot", setup.modelConfig.GetName(), params.GetName());
760  TString paramsZSnapName = TString::Format("%s__snapshot", setup.modelConfig_bonly.GetName());
761  RooFit::MsgLevel level = RooMsgService::instance().globalKillBelow();
762  RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
763  w->defineSet(paramsSnapName, params , true);
764  w->defineSet(paramsZSnapName, paramsZero ,true);
765  RooMsgService::instance().setGlobalKillBelow(level);
766 
767  // Create pdfs without nusiances terms, can be used for LEP tests statistics and
768  // for generating toys when not generating global observables
769  RooAbsPdf *factorizedPdf_s = setup.modelConfig.GetPdf(), *factorizedPdf_b = setup.modelConfig_bonly.GetPdf();
771  RooArgList constraints;
772  RooAbsPdf *factorizedPdf_s = utils::factorizePdf(*mc_s->GetObservables(), *mc_s->GetPdf(), constraints);
773  RooAbsPdf *factorizedPdf_b = utils::factorizePdf(*mc_b->GetObservables(), *mc_b->GetPdf(), constraints);
774  if (factorizedPdf_s != mc_s->GetPdf()) setup.cleanupList.addOwned(*factorizedPdf_s);
775  if (factorizedPdf_b != mc_b->GetPdf()) setup.cleanupList.addOwned(*factorizedPdf_b);
776  setup.modelConfig.SetPdf(*factorizedPdf_s);
777  setup.modelConfig_bonly.SetPdf(*factorizedPdf_b);
778  }
779 
780  if (testStat_ == "LEP") {
781  //SLR is evaluated using the central value of the nuisance parameters, so we have to put them in the parameter sets
782  if (withSystematics) {
783  if (!fitNuisances_) {
784  params.add(*mc_s->GetNuisanceParameters(), true);
785  paramsZero.addClone(*mc_b->GetNuisanceParameters(), true);
786  } else {
787  std::cerr << "ALERT: using LEP test statistics with --fitNuisances is not validated and most likely broken" << std::endl;
788  params.assignValueOnly(*mc_s->GetNuisanceParameters());
789  paramsZero.assignValueOnly(*mc_s->GetNuisanceParameters());
790  }
791  }
792  RooAbsPdf *pdfB = factorizedPdf_b;
793  if (poi.getSize() == 1) pdfB = factorizedPdf_s; // in this case we can remove the arbitrary constant from the test statistics.
795  if (!mc_s->GetPdf()->canBeExtended()) {
796  setup.qvar.reset(new SimplerLikelihoodRatioTestStat(*pdfB,*factorizedPdf_s, paramsZero, params));
797  } else {
798  setup.qvar.reset(new SimplerLikelihoodRatioTestStatOpt(*mc_s->GetObservables(), *pdfB, *factorizedPdf_s, paramsZero, params, withSystematics));
799  }
800  } else {
801  std::cerr << "ALERT: LEP test statistics without optimization not validated." << std::endl;
802  RooArgSet paramsSnap; params.snapshot(paramsSnap); // needs a snapshot
803  setup.qvar.reset(new SimpleLikelihoodRatioTestStat(*pdfB,*factorizedPdf_s));
804  ((SimpleLikelihoodRatioTestStat&)*setup.qvar).SetNullParameters(paramsZero); // Null is B
805  ((SimpleLikelihoodRatioTestStat&)*setup.qvar).SetAltParameters(paramsSnap);
806  }
807  } else if (testStat_ == "TEV") {
808  std::cerr << "ALERT: Tevatron test statistics not yet validated." << std::endl;
809  RooAbsPdf *pdfB = factorizedPdf_b;
810  if (poi.getSize() == 1) pdfB = factorizedPdf_s; // in this case we can remove the arbitrary constant from the test statistics.
812  setup.qvar.reset(new ProfiledLikelihoodRatioTestStatOpt(*mc_s->GetObservables(), *pdfB, *mc_s->GetPdf(), mc_s->GetNuisanceParameters(), paramsZero, params));
813  ((ProfiledLikelihoodRatioTestStatOpt&)*setup.qvar).setPrintLevel(verbose);
814  } else {
815  setup.qvar.reset(new RatioOfProfiledLikelihoodsTestStat(*mc_s->GetPdf(), *pdfB, setup.modelConfig.GetSnapshot()));
816  ((RatioOfProfiledLikelihoodsTestStat&)*setup.qvar).SetSubtractMLE(false);
817  }
818  } else if (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile") {
819  if (poi.getSize() != 1 && testStat_ != "Profile") {
820  throw std::logic_error("ERROR: modified profile likelihood definitions (LHC, LHCFC) do not make sense in more than one dimension");
821  }
825  if (testStat_ == "Profile") side = ProfiledLikelihoodTestStatOpt::twoSidedDef;
826  if (workingMode_ == MakeSignificance) r->setVal(0.0);
827  setup.qvar.reset(new ProfiledLikelihoodTestStatOpt(*mc_s->GetObservables(), *mc_s->GetPdf(), mc_s->GetNuisanceParameters(), params, poi, gobsParams,gobs, verbose, side));
828  } else {
829  std::cerr << "ALERT: LHC test statistics without optimization not validated." << std::endl;
830  setup.qvar.reset(new ProfileLikelihoodTestStat(*mc_s->GetPdf()));
831  if (testStat_ == "LHC") {
832  ((ProfileLikelihoodTestStat&)*setup.qvar).SetOneSided(true);
833  } else if (testStat_ == "LHCFC") {
834  throw std::invalid_argument("Test statistics LHCFC is not supported without optimization");
835  }
836  }
837  } else if (testStat_ == "MLZ") {
838  if (workingMode_ == MakeSignificance) r->setVal(0.0);
839  setup.qvar.reset(new BestFitSigmaTestStat(*mc_s->GetObservables(), *mc_s->GetPdf(), mc_s->GetNuisanceParameters(), params, verbose));
840  }
841 
842  RooAbsPdf *nuisancePdf = 0;
844  nuisancePdf = utils::makeNuisancePdf(*mc_s);
845  if (nuisancePdf) setup.cleanupList.addOwned(*nuisancePdf);
846  }
847  if (newToyMCSampler_) {
848  setup.toymcsampler.reset(new ToyMCSamplerOpt(*setup.qvar, nToys_, nuisancePdf, genNuisances_));
849  } else {
850  std::cerr << "ALERT: running with newToyMC=0 not validated." << std::endl;
851  setup.toymcsampler.reset(new ToyMCSampler(*setup.qvar, nToys_));
852  }
853 
854  if (!mc_b->GetPdf()->canBeExtended()) setup.toymcsampler->SetNEventsPerToy(1);
855 
856  if (nCpu_ > 0) {
857  std::cerr << "ALERT: running with proof not validated." << std::endl;
858  if (verbose > 1) std::cout << " Will use " << nCpu_ << " CPUs." << std::endl;
859  setup.pc.reset(new ProofConfig(*w, nCpu_, "", kFALSE));
860  setup.toymcsampler->SetProofConfig(setup.pc.get());
861  }
862 
863 
864  std::auto_ptr<HybridCalculator> hc(new HybridCalculator(data, setup.modelConfig, setup.modelConfig_bonly, setup.toymcsampler.get()));
865  if (genNuisances_ || !genGlobalObs_) {
866  if (withSystematics) {
867  setup.toymcsampler->SetGlobalObservables(*setup.modelConfig.GetNuisanceParameters());
868  (static_cast<HybridCalculator&>(*hc)).ForcePriorNuisanceNull(*nuisancePdf);
869  (static_cast<HybridCalculator&>(*hc)).ForcePriorNuisanceAlt(*nuisancePdf);
870  }
871  } else if (genGlobalObs_ && !genNuisances_) {
872  setup.toymcsampler->SetGlobalObservables(*setup.modelConfig.GetGlobalObservables());
873  hc->ForcePriorNuisanceNull(*w->pdf("__HybridNew_fake_nuisPdf__"));
874  hc->ForcePriorNuisanceAlt(*w->pdf("__HybridNew_fake_nuisPdf__"));
875  }
876 
877  // we need less B toys than S toys
879  // need only B toys. just keep a few S+B ones to avoid possible divide-by-zero errors somewhere
880  hc->SetToys(nToys_, int(0.01*nToys_)+1);
881  if (fullBToys_) {
882  hc->SetToys(nToys_, nToys_);
883  }
884  } else if (!CLs_) {
885  // we need only S+B toys to compute CLs+b
886  hc->SetToys(fullBToys_ ? nToys_ : int(0.01*nToys_)+1, nToys_);
887  //for two sigma bands need an equal number of B
888  if (expectedFromGrid_ && (fabs(0.5-quantileForExpectedFromGrid_)>=0.4) ) {
889  hc->SetToys(nToys_, nToys_);
890  }
891  } else {
892  // need both, but more S+B than B
893  hc->SetToys(fullBToys_ ? nToys_ : int(0.25*nToys_), nToys_);
894  //for two sigma bands need an equal number of B
895  if (expectedFromGrid_ && (fabs(0.5-quantileForExpectedFromGrid_)>=0.4) ) {
896  hc->SetToys(nToys_, nToys_);
897  }
898  }
899 
900  static const char * istr = "__HybridNew__importanceSamplingDensity";
902  std::cerr << "ALERT: running with importance sampling not validated (and probably not working)." << std::endl;
903  if(verbose > 1) std::cout << ">>> Enabling importance sampling for null hyp." << std::endl;
904  if(!withSystematics) {
905  throw std::invalid_argument("Importance sampling is not available without systematics");
906  }
907  RooArgSet importanceSnapshot;
908  importanceSnapshot.addClone(poi);
909  importanceSnapshot.addClone(*mc_s->GetNuisanceParameters());
910  if (verbose > 2) importanceSnapshot.Print("V");
911  hc->SetNullImportanceDensity(mc_b->GetPdf(), &importanceSnapshot);
912  }
914  std::cerr << "ALERT: running with importance sampling not validated (and probably not working)." << std::endl;
915  if(verbose > 1) std::cout << ">>> Enabling importance sampling for alt. hyp." << std::endl;
916  if(!withSystematics) {
917  throw std::invalid_argument("Importance sampling is not available without systematics");
918  }
919  if (w->pdf(istr) == 0) {
920  w->factory("__oneHalf__[0.5]");
921  RooAddPdf *sum = new RooAddPdf(istr, "fifty-fifty", *mc_s->GetPdf(), *mc_b->GetPdf(), *w->var("__oneHalf__"));
922  w->import(*sum);
923  }
924  RooArgSet importanceSnapshot;
925  importanceSnapshot.addClone(poi);
926  importanceSnapshot.addClone(*mc_s->GetNuisanceParameters());
927  if (verbose > 2) importanceSnapshot.Print("V");
928  hc->SetAltImportanceDensity(w->pdf(istr), &importanceSnapshot);
929  }
930 
931  return hc;
932 }
int i
Definition: DBlmapReader.cc:9
std::auto_ptr< RooStats::ProofConfig > pc
Definition: HybridNew.h:86
static bool importanceSamplingNull_
Definition: HybridNew.h:59
std::auto_ptr< RooStats::ToyMCSampler > toymcsampler
Definition: HybridNew.h:85
bool setAllConstant(const RooAbsCollection &coll, bool constant=true)
set all RooRealVars to constants. return true if at least one changed status
Definition: utils.cc:248
RooStats::ModelConfig modelConfig
Definition: HybridNew.h:83
static PFTauRenderPlugin instance
bool withSystematics
Definition: Combine.cc:68
RooAbsPdf * makeNuisancePdf(RooStats::ModelConfig &model, const char *name="nuisancePdf")
Definition: utils.cc:200
RooArgSet cleanupList
Definition: HybridNew.h:87
static WorkingMode workingMode_
Definition: HybridNew.h:40
static float quantileForExpectedFromGrid_
Definition: HybridNew.h:53
std::auto_ptr< RooStats::TestStatistic > qvar
Definition: HybridNew.h:84
RooAbsPdf * factorizePdf(const RooArgSet &observables, RooAbsPdf &pdf, RooArgList &constraints)
Definition: utils.cc:80
static bool CLs_
Definition: HybridNew.h:43
static unsigned int nCpu_
Definition: HybridNew.h:58
static unsigned int nToys_
Definition: HybridNew.h:41
static bool genGlobalObs_
Definition: HybridNew.h:46
static std::string testStat_
Definition: HybridNew.h:44
static bool optimizeProductPdf_
Definition: HybridNew.h:65
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
static bool optimizeTestStatistics_
Definition: HybridNew.h:66
static bool genNuisances_
Definition: HybridNew.h:46
RooStats::ModelConfig modelConfig_bonly
Definition: HybridNew.h:83
tuple cout
Definition: gather_cfg.py:121
static bool expectedFromGrid_
Definition: HybridNew.h:52
tuple level
Definition: testEve_cfg.py:34
static bool newToyMCSampler_
Definition: HybridNew.h:67
static bool importanceSamplingAlt_
Definition: HybridNew.h:59
static bool fitNuisances_
Definition: HybridNew.h:46
static bool fullBToys_
Definition: HybridNew.h:54
static const int ERROR
T w() const
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.

630  {
631  RooArgSet rValues;
632  mc_s->GetParametersOfInterest()->snapshot(rValues);
633  RooRealVar *r = dynamic_cast<RooRealVar*>(rValues.first());
634  if (rVal > r->getMax()) r->setMax(2*rVal);
635  r->setVal(rVal);
636  return create(w,mc_s,mc_b,data,rValues,setup);
637 }
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
T w() const
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
Definition: HybridNew.cc:639
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().

587  {
588  if (readHybridResults_) {
589  bool isProfile = (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile");
590  std::auto_ptr<RooStats::HypoTestResult> result(readToysFromFile(rVals));
591  std::pair<double, double> ret(-1,-1);
592  assert(result.get() != 0 && "Null result in HybridNew::eval(Workspace,...) after readToysFromFile");
593  if (expectedFromGrid_) {
595  result->SetTestStatisticData(result->GetTestStatisticData() + (isProfile ? -EPS : EPS));
596  } else if (!noUpdateGrid_) {
597  Setup setup;
598  std::auto_ptr<RooStats::HybridCalculator> hc = create(w, mc_s, mc_b, data, rVals, setup);
599  RooArgSet nullPOI(*setup.modelConfig_bonly.GetSnapshot());
600  if (isProfile) nullPOI.assignValueOnly(rVals);
601  double testStat = setup.qvar->Evaluate(data, nullPOI);
602  result->SetTestStatisticData(testStat + (isProfile ? -EPS : EPS));
603  }
604  ret = eval(*result, rVals);
605  return ret;
606  }
607 
609  RooLinkedListIter it = rVals.iterator();
610  for (RooRealVar *rIn = (RooRealVar*) it.Next(); rIn != 0; rIn = (RooRealVar*) it.Next()) {
611  RooRealVar *r = dynamic_cast<RooRealVar *>(mc_s->GetParametersOfInterest()->find(rIn->GetName()));
612  r->setVal(rIn->getVal());
613  if (verbose) std::cout << " " << r->GetName() << " = " << rIn->getVal() << " +/- " << r->getError() << std::endl;
614  }
615  std::auto_ptr<RooStats::HybridCalculator> hc(create(w, mc_s, mc_b, data, rVals, setup));
616  std::pair<double, double> ret = eval(*hc, rVals, adaptive, clsTarget);
617 
618  // add to plot
619  if (limitPlot_.get()) {
620  limitPlot_->Set(limitPlot_->GetN()+1);
621  limitPlot_->SetPoint(limitPlot_->GetN()-1, ((RooAbsReal*)rVals.first())->getVal(), ret.first);
622  limitPlot_->SetPointError(limitPlot_->GetN()-1, 0, ret.second);
623  }
624 
625  return ret;
626 }
#define EPS
Definition: HybridNew.cc:87
void applyExpectedQuantile(RooStats::HypoTestResult &hcres)
Definition: HybridNew.cc:1068
tuple result
Definition: query.py:137
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
Definition: HybridNew.cc:587
static std::string testStat_
Definition: HybridNew.h:44
static bool noUpdateGrid_
Definition: HybridNew.h:57
RooStats::HypoTestResult * readToysFromFile(const RooAbsCollection &rVals)
Definition: HybridNew.cc:1293
std::auto_ptr< TGraphErrors > limitPlot_
Definition: HybridNew.h:73
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
tuple cout
Definition: gather_cfg.py:121
static bool expectedFromGrid_
Definition: HybridNew.h:52
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
static bool readHybridResults_
Definition: HybridNew.h:50
T w() const
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
Definition: HybridNew.cc:639
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.

578  {
579  RooArgSet rValues;
580  mc_s->GetParametersOfInterest()->snapshot(rValues);
581  RooRealVar *r = dynamic_cast<RooRealVar*>(rValues.first());
582  if (rVal > r->getMax()) r->setMax(2*rVal);
583  r->setVal(rVal);
584  return eval(w,mc_s,mc_b,data,rValues,adaptive,clsTarget);
585 }
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
Definition: HybridNew.cc:587
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
T w() const
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.

935  {
936  std::auto_ptr<HypoTestResult> hcResult(evalGeneric(hc));
938  if (hcResult.get() == 0) {
939  std::cerr << "Hypotest failed" << std::endl;
940  return std::pair<double, double>(-1,-1);
941  }
942  if (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile") {
943  // I need to flip the P-values
944  hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()-EPS); // issue with < vs <= in discrete models
945  } else {
946  hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()+EPS); // issue with < vs <= in discrete models
947  hcResult->SetPValueIsRightTail(!hcResult->GetPValueIsRightTail());
948  }
949  std::pair<double,double> cls = eval(*hcResult, rVals);
950  if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tCLsplusb = ") << cls.first << " +/- " << cls.second << std::endl;
951  if (adaptive) {
952  if (CLs_) {
953  hc.SetToys(int(0.25*nToys_ + 1), nToys_);
954  }
955  else {
956  hc.SetToys(1, nToys_);
957  }
958  //for two sigma bands need an equal number of B
959  if (expectedFromGrid_ && (fabs(0.5-quantileForExpectedFromGrid_)>=0.4) ) {
960  hc.SetToys(nToys_, nToys_);
961  }
962  while (cls.second >= clsAccuracy_ && (clsTarget == -1 || fabs(cls.first-clsTarget) < 3*cls.second) ) {
963  std::auto_ptr<HypoTestResult> more(evalGeneric(hc));
964  more->SetBackgroundAsAlt(false);
965  if (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile") more->SetPValueIsRightTail(!more->GetPValueIsRightTail());
966  hcResult->Append(more.get());
968  cls = eval(*hcResult, rVals);
969  if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tCLsplusb = ") << cls.first << " +/- " << cls.second << std::endl;
970  }
971  } else if (iterations_ > 1) {
972  for (unsigned int i = 1; i < iterations_; ++i) {
973  std::auto_ptr<HypoTestResult> more(evalGeneric(hc));
974  more->SetBackgroundAsAlt(false);
975  if (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile") more->SetPValueIsRightTail(!more->GetPValueIsRightTail());
976  hcResult->Append(more.get());
978  cls = eval(*hcResult, rVals);
979  if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tCLsplusb = ") << cls.first << " +/- " << cls.second << std::endl;
980  }
981  }
982 
983  if (verbose > 0) {
984  std::cout <<
985  "\tCLs = " << hcResult->CLs() << " +/- " << hcResult->CLsError() << "\n" <<
986  "\tCLb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << "\n" <<
987  "\tCLsplusb = " << hcResult->CLsplusb() << " +/- " << hcResult->CLsplusbError() << "\n" <<
988  std::endl;
989  }
990 
991  perf_totalToysRun_ += (hcResult->GetAltDistribution()->GetSize() + hcResult->GetNullDistribution()->GetSize());
992 
993  if (!plot_.empty() && workingMode_ != MakeLimit) {
994  HypoTestPlot plot(*hcResult, 30);
995  TCanvas *c1 = new TCanvas("c1","c1");
996  plot.Draw();
997  c1->Print(plot_.c_str());
998  delete c1;
999  }
1000  if (saveHybridResult_) {
1001  TString name = TString::Format("HypoTestResult_mh%g",mass_);
1002  RooLinkedListIter it = rVals.iterator();
1003  for (RooRealVar *rIn = (RooRealVar*) it.Next(); rIn != 0; rIn = (RooRealVar*) it.Next()) {
1004  name += Form("_%s%g", rIn->GetName(), rIn->getVal());
1005  }
1006  name += Form("_%u", RooRandom::integer(std::numeric_limits<UInt_t>::max() - 1));
1007  writeToysHere->WriteTObject(new HypoTestResult(*hcResult), name);
1008  if (verbose) std::cout << "Hybrid result saved as " << name << " in " << writeToysHere->GetFile()->GetName() << " : " << writeToysHere->GetPath() << std::endl;
1009  }
1010 
1011  return cls;
1012 }
int i
Definition: DBlmapReader.cc:9
static unsigned int iterations_
Definition: HybridNew.h:49
RooStats::HypoTestResult * evalGeneric(RooStats::HybridCalculator &hc, bool forceNoFork=false)
Definition: HybridNew.cc:1170
long int integer
Definition: mlp_lapack.h:12
float mass_
Definition: HybridNew.h:70
static WorkingMode workingMode_
Definition: HybridNew.h:40
TDirectory * writeToysHere
Definition: Combine.cc:65
static bool saveHybridResult_
Definition: HybridNew.h:50
virtual const std::string & name() const
Definition: HybridNew.h:32
#define EPS
Definition: HybridNew.cc:87
static float quantileForExpectedFromGrid_
Definition: HybridNew.h:53
void applyExpectedQuantile(RooStats::HypoTestResult &hcres)
Definition: HybridNew.cc:1068
const T & max(const T &a, const T &b)
static bool CLs_
Definition: HybridNew.h:43
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
Definition: HybridNew.cc:587
static unsigned int nToys_
Definition: HybridNew.h:41
static std::string testStat_
Definition: HybridNew.h:44
unsigned int perf_totalToysRun_
Definition: HybridNew.h:76
static std::string plot_
Definition: HybridNew.h:61
perl if(1 lt scalar(@::datatypes))
Definition: edlooper.cc:31
tuple cout
Definition: gather_cfg.py:121
static bool expectedFromGrid_
Definition: HybridNew.h:52
static double clsAccuracy_
Definition: HybridNew.h:42
std::pair< double, double > HybridNew::eval ( const RooStats::HypoTestResult &  hcres,
const RooAbsCollection &  rVals 
)
private

Definition at line 1014 of file HybridNew.cc.

References eval().

1015 {
1016  double rVal = ((RooAbsReal*)rVals.first())->getVal();
1017  return eval(hcres,rVal);
1018 }
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
Definition: HybridNew.cc:587
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_.

1021 {
1022  if (testStat_ == "LHCFC") {
1023  RooStats::SamplingDistribution * bDistribution = hcres.GetNullDistribution(), * sDistribution = hcres.GetAltDistribution();
1024  const std::vector<Double_t> & bdist = bDistribution->GetSamplingDistribution();
1025  const std::vector<Double_t> & bweight = bDistribution->GetSampleWeights();
1026  const std::vector<Double_t> & sdist = sDistribution->GetSamplingDistribution();
1027  const std::vector<Double_t> & sweight = sDistribution->GetSampleWeights();
1028  Double_t data = hcres.GetTestStatisticData();
1029  std::vector<Double_t> absbdist(bdist.size()), abssdist(sdist.size());
1030  std::vector<Double_t> absbweight(bweight), abssweight(sweight);
1031  Double_t absdata;
1032  if (rule_ == "FC") {
1033  for (int i = 0, n = absbdist.size(); i < n; ++i) absbdist[i] = fabs(bdist[i]);
1034  for (int i = 0, n = abssdist.size(); i < n; ++i) abssdist[i] = fabs(sdist[i]);
1035  absdata = fabs(data)-2*minimizerTolerance_; // needed here since zeros are not exact by more than tolerance
1036  } else {
1037  for (int i = 0, n = absbdist.size(); i < n; ++i) absbdist[i] = max(0., bdist[i]);
1038  for (int i = 0, n = abssdist.size(); i < n; ++i) abssdist[i] = max(0., sdist[i]);
1039  absdata = max(0., data) - EPS;
1040  }
1041  if (rVal == 0) { // S+B toys are equal to B ones!
1042  abssdist.reserve(absbdist.size() + abssdist.size());
1043  abssdist.insert(abssdist.end(), absbdist.begin(), absbdist.end());
1044  abssweight.reserve(absbweight.size() + abssweight.size());
1045  abssweight.insert(abssweight.end(), absbweight.begin(), absbweight.end());
1046  }
1047  RooStats::HypoTestResult result;
1048  RooStats::SamplingDistribution *abssDist = new RooStats::SamplingDistribution("s","s",abssdist,abssweight);
1049  RooStats::SamplingDistribution *absbDist = new RooStats::SamplingDistribution("b","b",absbdist,absbweight);
1050  result.SetNullDistribution(absbDist);
1051  result.SetAltDistribution(abssDist);
1052  result.SetTestStatisticData(absdata);
1053  result.SetPValueIsRightTail(!result.GetPValueIsRightTail());
1054  if (CLs_) {
1055  return std::pair<double,double>(result.CLs(), result.CLsError());
1056  } else {
1057  return std::pair<double,double>(result.CLsplusb(), result.CLsplusbError());
1058  }
1059  } else {
1060  if (CLs_) {
1061  return std::pair<double,double>(hcres.CLs(), hcres.CLsError());
1062  } else {
1063  return std::pair<double,double>(hcres.CLsplusb(), hcres.CLsplusbError());
1064  }
1065  }
1066 }
int i
Definition: DBlmapReader.cc:9
static std::string rule_
Definition: HybridNew.h:44
static float minimizerTolerance_
Definition: HybridNew.h:63
#define EPS
Definition: HybridNew.cc:87
const T & max(const T &a, const T &b)
static bool CLs_
Definition: HybridNew.h:43
tuple result
Definition: query.py:137
static std::string testStat_
Definition: HybridNew.h:44
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
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().

1170  {
1171  if (fork_ && !noFork) return evalWithFork(hc);
1172  else {
1173  TStopwatch timer; timer.Start();
1174  RooStats::HypoTestResult * ret = hc.GetHypoTest();
1175  if (runtimedef::get("HybridNew_Timing")) std::cout << "Evaluated toys in " << timer.RealTime() << " s " << std::endl;
1176  return ret;
1177  }
1178 }
static unsigned int fork_
Definition: HybridNew.h:58
int get(const char *name)
tuple cout
Definition: gather_cfg.py:121
RooStats::HypoTestResult * evalWithFork(RooStats::HybridCalculator &hc)
Definition: HybridNew.cc:1180
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().

1180  {
1181  TStopwatch timer;
1182  std::auto_ptr<RooStats::HypoTestResult> result(0);
1183  char *tmpfile = tempnam(NULL,"rstat");
1184  unsigned int ich = 0;
1185  std::vector<UInt_t> newSeeds(fork_);
1186  for (ich = 0; ich < fork_; ++ich) {
1188  if (!fork()) break; // spawn children (but only in the parent thread)
1189  }
1190  if (ich == fork_) { // if i'm the parent
1191  int cstatus, ret;
1192  do {
1193  do { ret = waitpid(-1, &cstatus, 0); } while (ret == -1 && errno == EINTR);
1194  } while (ret != -1);
1195  if (ret == -1 && errno != ECHILD) throw std::runtime_error("Didn't wait for child");
1196  for (ich = 0; ich < fork_; ++ich) {
1197  TFile *f = TFile::Open(TString::Format("%s.%d.root", tmpfile, ich));
1198  if (f == 0) throw std::runtime_error(TString::Format("Child didn't leave output file %s.%d.root", tmpfile, ich).Data());
1199  RooStats::HypoTestResult *res = dynamic_cast<RooStats::HypoTestResult *>(f->Get("result"));
1200  if (res == 0) throw std::runtime_error(TString::Format("Child output file %s.%d.root is corrupted", tmpfile, ich).Data());
1201  if (result.get()) result->Append(res); else result.reset(dynamic_cast<RooStats::HypoTestResult *>(res->Clone()));
1202  f->Close();
1203  unlink(TString::Format("%s.%d.root", tmpfile, ich).Data());
1204  unlink(TString::Format("%s.%d.out.txt", tmpfile, ich).Data());
1205  unlink(TString::Format("%s.%d.err.txt", tmpfile, ich).Data());
1206  }
1207  } else {
1208  RooRandom::randomGenerator()->SetSeed(newSeeds[ich]);
1209  freopen(TString::Format("%s.%d.out.txt", tmpfile, ich).Data(), "w", stdout);
1210  freopen(TString::Format("%s.%d.err.txt", tmpfile, ich).Data(), "w", stderr);
1211  std::cout << " I'm child " << ich << ", seed " << newSeeds[ich] << std::endl;
1212  RooStats::HypoTestResult *hcResult = evalGeneric(hc, /*noFork=*/true);
1213  TFile *f = TFile::Open(TString::Format("%s.%d.root", tmpfile, ich), "RECREATE");
1214  f->WriteTObject(hcResult, "result");
1215  f->ls();
1216  f->Close();
1217  fflush(stdout); fflush(stderr);
1218  std::cout << "And I'm done" << std::endl;
1219  throw std::runtime_error("done"); // I have to throw instead of exiting, otherwise there's no proper stack unwinding
1220  // and deleting of intermediate objects, and when the statics get deleted it crashes
1221  // in 5.27.06 (but not in 5.28)
1222  }
1223  free(tmpfile);
1224  if (verbose > 1) { std::cout << " Evaluation of p-values done in " << timer.RealTime() << " s" << std::endl; }
1225  return result.release();
1226 }
static unsigned int fork_
Definition: HybridNew.h:58
RooStats::HypoTestResult * evalGeneric(RooStats::HybridCalculator &hc, bool forceNoFork=false)
Definition: HybridNew.cc:1170
long int integer
Definition: mlp_lapack.h:12
#define NULL
Definition: scimark2.h:8
const T & max(const T &a, const T &b)
tuple result
Definition: query.py:137
double f[11][100]
tuple cout
Definition: gather_cfg.py:121
virtual const std::string& HybridNew::name ( void  ) const
inlinevirtual

Implements LimitAlgo.

Definition at line 32 of file HybridNew.h.

Referenced by eval(), readGrid(), and runSignificance().

32  {
33  static const std::string name("HybridNew");
34  return name;
35  }
virtual const std::string & name() const
Definition: HybridNew.h:32
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().

1345  {
1346  if (rValues_.getSize() != 1) throw std::runtime_error("Running limits with grid only works in one dimension for the moment");
1347  clearGrid();
1348 
1349  TIter next(toyDir->GetListOfKeys()); TKey *k; const char *poiName = rValues_.first()->GetName();
1350  while ((k = (TKey *) next()) != 0) {
1351  TString name(k->GetName());
1352  if (name.Index("HypoTestResult_mh") == 0) {
1353  if (name.Index(TString::Format("HypoTestResult_mh%g_%s",mass_,poiName)) != 0 || name.Index("_", name.Index("_")+1) == -1) continue;
1354  name.ReplaceAll(TString::Format("HypoTestResult_mh%g_%s",mass_,poiName),""); // remove the prefix
1355  if (name.Index("_") == -1) continue; // check if it has the _<uniqueId> postfix
1356  name.Remove(name.Index("_"),name.Length()); // remove it before calling atof
1357  } else if (name.Index("HypoTestResult_") == 0) {
1358  // let's put a warning here, since results of this form were supported in the past
1359  std::cout << "HybridNew::readGrid: HypoTestResult with non-conformant name " << name << " will be skipped" << std::endl;
1360  continue;
1361  } else continue;
1362  double rVal = atof(name.Data());
1363  if (rVal < rMin || rVal > rMax) continue;
1364  if (verbose > 2) std::cout << " Do " << k->GetName() << " -> " << name << " --> " << rVal << std::endl;
1365  RooStats::HypoTestResult *toy = dynamic_cast<RooStats::HypoTestResult *>(toyDir->Get(k->GetName()));
1366  RooStats::HypoTestResult *&merge = grid_[rVal];
1367  if (merge == 0) merge = new RooStats::HypoTestResult(*toy);
1368  else merge->Append(toy);
1369  merge->ResetBit(1);
1370  }
1371  if (verbose > 1) {
1372  std::cout << "GRID, as is." << std::endl;
1373  typedef std::map<double, RooStats::HypoTestResult *>::iterator point;
1374  for (point it = grid_.begin(), ed = grid_.end(); it != ed; ++it) {
1375  std::cout << " - " << it->first << " (CLs = " << it->second->CLs() << " +/- " << it->second->CLsError() << ")" << std::endl;
1376  }
1377  }
1378 }
float mass_
Definition: HybridNew.h:70
virtual const std::string & name() const
Definition: HybridNew.h:32
void clearGrid()
Definition: HybridNew.cc:1509
std::map< double, RooStats::HypoTestResult * > grid_
Definition: HybridNew.h:119
int k[5][pyjets_maxn]
tuple cout
Definition: gather_cfg.py:121
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
static RooArgSet rValues_
Definition: HybridNew.h:48
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().

1293  {
1294  if (!readToysFromHere) throw std::logic_error("Cannot use readHypoTestResult: option toysFile not specified, or input file empty");
1295  TDirectory *toyDir = readToysFromHere->GetDirectory("toys");
1296  if (!toyDir) throw std::logic_error("Cannot use readHypoTestResult: empty toy dir in input file empty");
1297  if (verbose) std::cout << "Reading toys for ";
1298  TString prefix1 = TString::Format("HypoTestResult_mh%g",mass_);
1299  TString prefix2 = TString::Format("HypoTestResult");
1300  RooLinkedListIter it = rVals.iterator();
1301  for (RooRealVar *rIn = (RooRealVar*) it.Next(); rIn != 0; rIn = (RooRealVar*) it.Next()) {
1302  if (verbose) std::cout << rIn->GetName() << " = " << rIn->getVal() << " ";
1303  prefix1 += Form("_%s%g", rIn->GetName(), rIn->getVal());
1304  prefix2 += Form("_%s%g", rIn->GetName(), rIn->getVal());
1305  }
1306  if (verbose) std::cout << std::endl;
1307  std::auto_ptr<RooStats::HypoTestResult> ret;
1308  TIter next(toyDir->GetListOfKeys()); TKey *k;
1309  while ((k = (TKey *) next()) != 0) {
1310  if (TString(k->GetName()).Index(prefix1) != 0 && TString(k->GetName()).Index(prefix2) != 0) continue;
1311  RooStats::HypoTestResult *toy = dynamic_cast<RooStats::HypoTestResult *>(toyDir->Get(k->GetName()));
1312  if (toy == 0) continue;
1313  if (verbose > 1) std::cout << " - " << k->GetName() << std::endl;
1314  if (ret.get() == 0) {
1315  ret.reset(new RooStats::HypoTestResult(*toy));
1316  } else {
1317  ret->Append(toy);
1318  }
1319  }
1320 
1321  if (ret.get() == 0) {
1322  std::cout << "ERROR: parameter point not found in input root file.\n";
1323  rVals.Print("V");
1324  if (verbose > 0) toyDir->ls();
1325  std::cout << "ERROR: parameter point not found in input root file" << std::endl;
1326  throw std::invalid_argument("Missing input");
1327  }
1328  if (verbose > 0) {
1329  std::cout <<
1330  "\tCLs = " << ret->CLs() << " +/- " << ret->CLsError() << "\n" <<
1331  "\tCLb = " << ret->CLb() << " +/- " << ret->CLbError() << "\n" <<
1332  "\tCLsplusb = " << ret->CLsplusb() << " +/- " << ret->CLsplusbError() << "\n" <<
1333  std::endl;
1334  if (!plot_.empty() && workingMode_ != MakeLimit) {
1335  HypoTestPlot plot(*ret, 30);
1336  TCanvas *c1 = new TCanvas("c1","c1");
1337  plot.Draw();
1338  c1->Print(plot_.c_str());
1339  delete c1;
1340  }
1341  }
1342  return ret.release();
1343 }
TDirectory * readToysFromHere
Definition: Combine.cc:66
float mass_
Definition: HybridNew.h:70
static WorkingMode workingMode_
Definition: HybridNew.h:40
int k[5][pyjets_maxn]
static std::string plot_
Definition: HybridNew.h:61
tuple cout
Definition: gather_cfg.py:121
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_.

251  {
254  perf_totalToysRun_ = 0; // reset performance counter
255  if (rValues_.getSize() == 0) setupPOI(mc_s);
256  switch (workingMode_) {
257  case MakeLimit: return runLimit(w, mc_s, mc_b, data, limit, limitErr, hint);
258  case MakeSignificance: return runSignificance(w, mc_s, mc_b, data, limit, limitErr, hint);
259  case MakePValues: return runSinglePoint(w, mc_s, mc_b, data, limit, limitErr, hint);
260  case MakeTestStatistics:
262  return runTestStatistics(w, mc_s, mc_b, data, limit, limitErr, hint);
263  }
264  assert("Shouldn't get here" == 0);
265 }
#define DEBUG
virtual bool runSignificance(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint)
Definition: HybridNew.cc:267
virtual bool runLimit(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint)
Definition: HybridNew.cc:318
static const int WARNING
static float minimizerTolerance_
Definition: HybridNew.h:63
static std::string minimizerAlgo_
Definition: HybridNew.h:62
static WorkingMode workingMode_
Definition: HybridNew.h:40
virtual bool runTestStatistics(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint)
Definition: HybridNew.cc:556
void setupPOI(RooStats::ModelConfig *mc_s)
Definition: HybridNew.cc:220
Setup Minimizer configuration on creation, reset the previous one on destruction. ...
unsigned int perf_totalToysRun_
Definition: HybridNew.h:76
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
virtual bool runSinglePoint(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint)
Definition: HybridNew.cc:546
static RooArgSet rValues_
Definition: HybridNew.h:48
T w() const
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(), vdt::x, and detailsBasic3DVector::y.

Referenced by run().

318  {
319  if (mc_s->GetParametersOfInterest()->getSize() != 1) throw std::logic_error("Cannot run limits in more than one dimension, for the moment.");
320  RooRealVar *r = dynamic_cast<RooRealVar *>(mc_s->GetParametersOfInterest()->first()); r->setConstant(true);
321  w->loadSnapshot("clean");
322 
323  if ((hint != 0) && (*hint > r->getMin())) {
324  r->setMax(std::min<double>(3.0 * (*hint), r->getMax()));
325  r->setMin(std::max<double>(0.3 * (*hint), r->getMin()));
326  }
327 
328  typedef std::pair<double,double> CLs_t;
329 
330  double clsTarget = 1 - cl;
331  CLs_t clsMin(1,0), clsMax(0,0), clsMid(0,0);
332  double rMin = r->getMin(), rMax = r->getMax();
333  limit = 0.5*(rMax + rMin);
334  limitErr = 0.5*(rMax - rMin);
335  bool done = false;
336 
337  TF1 expoFit("expoFit","[0]*exp([1]*(x-[2]))", rMin, rMax);
338 
339  if (readHybridResults_) {
340  if (verbose > 0) std::cout << "Search for upper limit using pre-computed grid of p-values" << std::endl;
341 
342  if (!gridFile_.empty()) {
343  if (grid_.empty()) {
344  std::auto_ptr<TFile> gridFile(TFile::Open(gridFile_.c_str()));
345  if (gridFile.get() == 0) throw std::runtime_error(("Can't open grid file "+gridFile_).c_str());
346  TDirectory *toyDir = gridFile->GetDirectory("toys");
347  if (!toyDir) throw std::logic_error("Cannot use readHypoTestResult: empty toy dir in input file empty");
348  readGrid(toyDir, rMinSet_ ? rMin : -99e99, rMaxSet_ ? rMax :+99e99);
349  }
350  if (grid_.size() <= 1) throw std::logic_error("The grid must contain at least 2 points.");
351  if (noUpdateGrid_) {
352  if (testStat_ == "LHCFC" && rule_ == "FC" && (saveGrid_ || lowerLimit_)) {
353  std::cout << "Will have to re-run points for which the test statistics was set to zero" << std::endl;
354  updateGridDataFC(w, mc_s, mc_b, data, !fullGrid_, clsTarget);
355  } else {
356  std::cout << "Will use the test statistics that had already been computed" << std::endl;
357  }
358  } else {
359  updateGridData(w, mc_s, mc_b, data, !fullGrid_, clsTarget);
360  }
361  } else throw std::logic_error("When setting a limit reading results from file, a grid file must be specified with option --grid");
362  if (grid_.size() <= 1) throw std::logic_error("The grid must contain at least 2 points.");
363 
364  useGrid();
365 
366  double minDist=1e3;
367  double nearest = 0;
368  rMin = limitPlot_->GetX()[0];
369  rMax = limitPlot_->GetX()[limitPlot_->GetN()-1];
370  for (int i = 0, n = limitPlot_->GetN(); i < n; ++i) {
371  double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i], ey = limitPlot_->GetErrorY(i);
372  if (verbose > 0) printf(" r %.2f: %s = %6.4f +/- %6.4f\n", x, CLs_ ? "CLs" : "CLsplusb", y, ey);
373  if (saveGrid_) { limit = x; limitErr = ey; Combine::commitPoint(false, y); }
374  if (y-3*max(ey,0.01) >= clsTarget) { rMin = x; clsMin = CLs_t(y,ey); }
375  if (fabs(y-clsTarget) < minDist) { nearest = x; minDist = fabs(y-clsTarget); }
376  rMax = x; clsMax = CLs_t(y,ey);
377  if (y+3*max(ey,0.01) <= clsTarget && !fullGrid_) break;
378  }
379  limit = nearest;
380  if (verbose > 0) std::cout << " after scan x ~ " << limit << ", bounds [ " << rMin << ", " << rMax << "]" << std::endl;
381  limitErr = std::max(limit-rMin, rMax-limit);
382  expoFit.SetRange(rMin,rMax);
383 
384  if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
385  if (verbose > 1) std::cout << " reached accuracy " << limitErr << " below " << std::max(rAbsAccuracy_, rRelAccuracy_ * limit) << std::endl;
386  done = true;
387  }
388  } else {
389  limitPlot_.reset(new TGraphErrors());
390 
391  if (verbose > 0) std::cout << "Search for upper limit to the limit" << std::endl;
392  for (int tries = 0; tries < 6; ++tries) {
393  clsMax = eval(w, mc_s, mc_b, data, rMax);
394  if (lowerLimit_) break; // we can't search for lower limits this way
395  if (clsMax.first == 0 || clsMax.first + 3 * fabs(clsMax.second) < clsTarget ) break;
396  rMax += rMax;
397  if (tries == 5) {
398  std::cerr << "Cannot set higher limit: at " << r->GetName() << " = " << rMax << " still get " << (CLs_ ? "CLs" : "CLsplusb") << " = " << clsMax.first << std::endl;
399  return false;
400  }
401  }
402  if (verbose > 0) std::cout << "Search for lower limit to the limit" << std::endl;
403  clsMin = (CLs_ && rMin == 0 ? CLs_t(1,0) : eval(w, mc_s, mc_b, data, rMin));
404  if (!lowerLimit_ && clsMin.first != 1 && clsMin.first - 3 * fabs(clsMin.second) < clsTarget) {
405  if (CLs_) {
406  rMin = 0;
407  clsMin = CLs_t(1,0); // this is always true for CLs
408  } else {
409  rMin = -rMax / 4;
410  for (int tries = 0; tries < 6; ++tries) {
411  clsMin = eval(w, mc_s, mc_b, data, rMin);
412  if (clsMin.first == 1 || clsMin.first - 3 * fabs(clsMin.second) > clsTarget) break;
413  rMin += rMin;
414  if (tries == 5) {
415  std::cerr << "Cannot set lower limit: at " << r->GetName() << " = " << rMin << " still get " << (CLs_ ? "CLs" : "CLsplusb") << " = " << clsMin.first << std::endl;
416  return false;
417  }
418  }
419  }
420  }
421 
422  if (verbose > 0) std::cout << "Now doing proper bracketing & bisection" << std::endl;
423  do {
424  // determine point by bisection or interpolation
425  limit = 0.5*(rMin+rMax); limitErr = 0.5*(rMax-rMin);
426  if (algo_ == "logSecant" && clsMax.first != 0 && clsMin.first != 0) {
427  double logMin = log(clsMin.first), logMax = log(clsMax.first), logTarget = log(clsTarget);
428  limit = rMin + (rMax-rMin) * (logTarget - logMin)/(logMax - logMin);
429  if (clsMax.second != 0 && clsMin.second != 0) {
430  limitErr = hypot((logTarget-logMax) * (clsMin.second/clsMin.first), (logTarget-logMin) * (clsMax.second/clsMax.first));
431  limitErr *= (rMax-rMin)/((logMax-logMin)*(logMax-logMin));
432  }
433  // disallow "too precise" interpolations
434  if (limitErr < interpAccuracy_ * (rMax-rMin)) limitErr = interpAccuracy_ * (rMax-rMin);
435  }
436  r->setError(limitErr);
437 
438  // exit if reached accuracy on r
439  if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
440  if (verbose > 1) std::cout << " reached accuracy " << limitErr << " below " << std::max(rAbsAccuracy_, rRelAccuracy_ * limit) << std::endl;
441  done = true; break;
442  }
443 
444  // evaluate point
445  clsMid = eval(w, mc_s, mc_b, data, limit, true, clsTarget);
446  if (clsMid.second == -1) {
447  std::cerr << "Hypotest failed" << std::endl;
448  return false;
449  }
450 
451  // if sufficiently far away, drop one of the points
452  if (fabs(clsMid.first-clsTarget) >= 2*clsMid.second) {
453  if ((clsMid.first>clsTarget) == (clsMax.first>clsTarget)) {
454  rMax = limit; clsMax = clsMid;
455  } else {
456  rMin = limit; clsMin = clsMid;
457  }
458  } else {
459  if (verbose > 0) std::cout << "Trying to move the interval edges closer" << std::endl;
460  double rMinBound = rMin, rMaxBound = rMax;
461  // try to reduce the size of the interval
462  while (clsMin.second == 0 || fabs(rMin-limit) > std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
463  rMin = 0.5*(rMin+limit);
464  clsMin = eval(w, mc_s, mc_b, data, rMin, true, clsTarget);
465  if (fabs(clsMin.first-clsTarget) <= 2*clsMin.second) break;
466  rMinBound = rMin;
467  }
468  while (clsMax.second == 0 || fabs(rMax-limit) > std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
469  rMax = 0.5*(rMax+limit);
470  clsMax = eval(w, mc_s, mc_b, data, rMax, true, clsTarget);
471  if (fabs(clsMax.first-clsTarget) <= 2*clsMax.second) break;
472  rMaxBound = rMax;
473  }
474  expoFit.SetRange(rMinBound,rMaxBound);
475  break;
476  }
477  } while (true);
478 
479  }
480 
481  if (!done) { // didn't reach accuracy with scan, now do fit
482  double rMinBound, rMaxBound; expoFit.GetRange(rMinBound, rMaxBound);
483  if (verbose) {
484  std::cout << "\n -- HybridNew, before fit -- \n";
485  std::cout << "Limit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " [" << rMinBound << ", " << rMaxBound << "]\n";
486  }
487 
488  expoFit.FixParameter(0,clsTarget);
489  double clsmaxfirst = clsMax.first;
490  if ( clsmaxfirst == 0.0 ) clsmaxfirst = 0.005;
491  double par1guess = log(clsmaxfirst/clsMin.first)/(rMax-rMin);
492  expoFit.SetParameter(1,par1guess);
493  expoFit.SetParameter(2,limit);
494  limitErr = std::max(fabs(rMinBound-limit), fabs(rMaxBound-limit));
495  int npoints = 0;
496  for (int j = 0; j < limitPlot_->GetN(); ++j) {
497  if (limitPlot_->GetX()[j] >= rMinBound && limitPlot_->GetX()[j] <= rMaxBound) npoints++;
498  }
499  for (int i = 0, imax = (readHybridResults_ ? 0 : 8); i <= imax; ++i, ++npoints) {
500  limitPlot_->Sort();
501  limitPlot_->Fit(&expoFit,(verbose <= 1 ? "QNR EX0" : "NR EX0"));
502  if (verbose) {
503  std::cout << "Fit to " << npoints << " points: " << expoFit.GetParameter(2) << " +/- " << expoFit.GetParError(2) << std::endl;
504  }
505  if ((rMinBound < expoFit.GetParameter(2)) && (expoFit.GetParameter(2) < rMaxBound) && (expoFit.GetParError(2) < 0.5*(rMaxBound-rMinBound))) {
506  // sanity check fit result
507  limit = expoFit.GetParameter(2);
508  limitErr = expoFit.GetParError(2);
509  if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) break;
510  }
511  // add one point in the interval.
512  double rTry = RooRandom::uniform()*(rMaxBound-rMinBound)+rMinBound;
513  if (i != imax) eval(w, mc_s, mc_b, data, rTry, true, clsTarget);
514  }
515  }
516 
517  if (!plot_.empty() && limitPlot_.get()) {
518  TCanvas *c1 = new TCanvas("c1","c1");
519  limitPlot_->Sort();
520  limitPlot_->SetLineWidth(2);
521  double xmin = r->getMin(), xmax = r->getMax();
522  for (int j = 0; j < limitPlot_->GetN(); ++j) {
523  if (limitPlot_->GetY()[j] > 1.4*clsTarget || limitPlot_->GetY()[j] < 0.6*clsTarget) continue;
524  xmin = std::min(limitPlot_->GetX()[j], xmin);
525  xmax = std::max(limitPlot_->GetX()[j], xmax);
526  }
527  limitPlot_->GetXaxis()->SetRangeUser(xmin,xmax);
528  limitPlot_->GetYaxis()->SetRangeUser(0.5*clsTarget, 1.5*clsTarget);
529  limitPlot_->Draw("AP");
530  expoFit.Draw("SAME");
531  TLine line(limitPlot_->GetX()[0], clsTarget, limitPlot_->GetX()[limitPlot_->GetN()-1], clsTarget);
532  line.SetLineColor(kRed); line.SetLineWidth(2); line.Draw();
533  line.DrawLine(limit, 0, limit, limitPlot_->GetY()[0]);
534  line.SetLineWidth(1); line.SetLineStyle(2);
535  line.DrawLine(limit-limitErr, 0, limit-limitErr, limitPlot_->GetY()[0]);
536  line.DrawLine(limit+limitErr, 0, limit+limitErr, limitPlot_->GetY()[0]);
537  c1->Print(plot_.c_str());
538  }
539 
540  std::cout << "\n -- Hybrid New -- \n";
541  std::cout << "Limit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " @ " << cl * 100 << "% CL\n";
542  if (verbose > 1) std::cout << "Total toys: " << perf_totalToysRun_ << std::endl;
543  return true;
544 }
int i
Definition: DBlmapReader.cc:9
static std::string rule_
Definition: HybridNew.h:44
bool lowerLimit_
Definition: Combine.cc:70
void updateGridDataFC(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, bool smart, double clsTarget)
Definition: HybridNew.cc:1441
#define min(a, b)
Definition: mlp_lapack.h:161
static std::string algo_
Definition: HybridNew.h:60
static double interpAccuracy_
Definition: HybridNew.h:42
static double rAbsAccuracy_
Definition: HybridNew.h:42
static void commitPoint(bool expected, float quantile)
Save a point into the output tree. Usually if expected = false, quantile should be set to -1 (except ...
Definition: Combine.cc:557
void updateGridData(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, bool smart, double clsTarget)
Definition: HybridNew.cc:1379
const T & max(const T &a, const T &b)
static bool CLs_
Definition: HybridNew.h:43
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
Definition: HybridNew.cc:587
int j
Definition: DBlmapReader.cc:9
static bool saveGrid_
Definition: HybridNew.h:56
std::map< double, RooStats::HypoTestResult * > grid_
Definition: HybridNew.h:119
static bool fullGrid_
Definition: HybridNew.h:55
static bool rMaxSet_
Definition: HybridNew.h:69
float cl
Definition: Combine.cc:71
static std::string testStat_
Definition: HybridNew.h:44
static bool noUpdateGrid_
Definition: HybridNew.h:57
void useGrid()
Definition: HybridNew.cc:1489
static std::string gridFile_
Definition: HybridNew.h:51
unsigned int perf_totalToysRun_
Definition: HybridNew.h:76
std::auto_ptr< TGraphErrors > limitPlot_
Definition: HybridNew.h:73
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
static std::string plot_
Definition: HybridNew.h:61
static double rRelAccuracy_
Definition: HybridNew.h:42
static bool rMinSet_
Definition: HybridNew.h:68
tuple cout
Definition: gather_cfg.py:121
void readGrid(TDirectory *directory, double rMin, double rMax)
Definition: HybridNew.cc:1345
x
Definition: VDTMath.h:216
static bool readHybridResults_
Definition: HybridNew.h:50
T w() const
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().

267  {
269  std::auto_ptr<RooStats::HybridCalculator> hc(create(w, mc_s, mc_b, data, rValues_, setup));
270  std::auto_ptr<HypoTestResult> hcResult;
271  if (readHybridResults_) {
272  hcResult.reset(readToysFromFile(rValues_));
273  } else {
274  hcResult.reset(evalGeneric(*hc));
275  for (unsigned int i = 1; i < iterations_; ++i) {
276  std::auto_ptr<HypoTestResult> more(evalGeneric(*hc));
277  hcResult->Append(more.get());
278  if (verbose) std::cout << "\t1 - CLb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << std::endl;
279  }
280  }
281  if (hcResult.get() == 0) {
282  std::cerr << "Hypotest failed" << std::endl;
283  return false;
284  }
285  if (saveHybridResult_) {
286  TString name = TString::Format("HypoTestResult_mh%g",mass_);
287  RooLinkedListIter it = rValues_.iterator();
288  for (RooRealVar *rIn = (RooRealVar*) it.Next(); rIn != 0; rIn = (RooRealVar*) it.Next()) {
289  name += Form("_%s%g", rIn->GetName(), rIn->getVal());
290  }
291  name += Form("_%u", RooRandom::integer(std::numeric_limits<UInt_t>::max() - 1));
292  writeToysHere->WriteTObject(new HypoTestResult(*hcResult), name);
293  if (verbose) std::cout << "Hybrid result saved as " << name << " in " << writeToysHere->GetFile()->GetName() << " : " << writeToysHere->GetPath() << std::endl;
294  }
295  if (verbose > 1) {
296  std::cout << "Observed test statistics in data: " << hcResult->GetTestStatisticData() << std::endl;
297  std::cout << "Background-only toys sampled: " << hcResult->GetNullDistribution()->GetSize() << std::endl;
298  }
300  // I don't need to flip the P-values for significances, only for limits
301  hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()+EPS); // issue with < vs <= in discrete models
302  double sig = hcResult->Significance();
303  double sigHi = RooStats::PValueToSignificance( (hcResult->CLb() - hcResult->CLbError()) ) - sig;
304  double sigLo = RooStats::PValueToSignificance( (hcResult->CLb() + hcResult->CLbError()) ) - sig;
305  if (reportPVal_) {
306  limit = hcResult->NullPValue();
307  limitErr = hcResult->NullPValueError();
308  } else {
309  limit = sig;
310  limitErr = std::max(sigHi,-sigLo);
311  }
312  std::cout << "\n -- Hybrid New -- \n";
313  std::cout << "Significance: " << sig << " " << sigLo << "/+" << sigHi << "\n";
314  std::cout << "Null p-value: " << hcResult->NullPValue() << " +/- " << hcResult->NullPValueError() << "\n";
315  return std::isfinite(limit);
316 }
int i
Definition: DBlmapReader.cc:9
static unsigned int iterations_
Definition: HybridNew.h:49
RooStats::HypoTestResult * evalGeneric(RooStats::HybridCalculator &hc, bool forceNoFork=false)
Definition: HybridNew.cc:1170
long int integer
Definition: mlp_lapack.h:12
float mass_
Definition: HybridNew.h:70
TDirectory * writeToysHere
Definition: Combine.cc:65
static bool saveHybridResult_
Definition: HybridNew.h:50
virtual const std::string & name() const
Definition: HybridNew.h:32
#define EPS
Definition: HybridNew.cc:87
void applyExpectedQuantile(RooStats::HypoTestResult &hcres)
Definition: HybridNew.cc:1068
const T & max(const T &a, const T &b)
static bool reportPVal_
Definition: HybridNew.h:45
RooStats::HypoTestResult * readToysFromFile(const RooAbsCollection &rVals)
Definition: HybridNew.cc:1293
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
tuple cout
Definition: gather_cfg.py:121
static bool expectedFromGrid_
Definition: HybridNew.h:52
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
static bool readHybridResults_
Definition: HybridNew.h:50
static RooArgSet rValues_
Definition: HybridNew.h:48
T w() const
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
Definition: HybridNew.cc:639
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().

546  {
547  std::pair<double, double> result = eval(w, mc_s, mc_b, data, rValues_, clsAccuracy_ != 0);
548  std::cout << "\n -- Hybrid New -- \n";
549  std::cout << (CLs_ ? "CLs = " : "CLsplusb = ") << result.first << " +/- " << result.second << std::endl;
550  if (verbose > 1) std::cout << "Total toys: " << perf_totalToysRun_ << std::endl;
551  limit = result.first;
552  limitErr = result.second;
553  return true;
554 }
static bool CLs_
Definition: HybridNew.h:43
tuple result
Definition: query.py:137
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
Definition: HybridNew.cc:587
unsigned int perf_totalToysRun_
Definition: HybridNew.h:76
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
perl if(1 lt scalar(@::datatypes))
Definition: edlooper.cc:31
tuple cout
Definition: gather_cfg.py:121
static RooArgSet rValues_
Definition: HybridNew.h:48
static double clsAccuracy_
Definition: HybridNew.h:42
T w() const
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().

556  {
557  bool isProfile = (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile");
559  std::auto_ptr<RooStats::HypoTestResult> result(readToysFromFile(rValues_));
561  limit = -2 * result->GetTestStatisticData();
562  } else {
564  std::auto_ptr<RooStats::HybridCalculator> hc(create(w, mc_s, mc_b, data, rValues_, setup));
565  RooArgSet nullPOI(*setup.modelConfig_bonly.GetSnapshot());
566  if (isProfile) {
568  nullPOI.assignValueOnly(rValues_);
569  }
570  limit = -2 * setup.qvar->Evaluate(data, nullPOI);
571  }
572  if (isProfile) limit = -limit; // there's a sign flip for these two
573  std::cout << "\n -- Hybrid New -- \n";
574  std::cout << "-2 ln Q_{"<< testStat_<<"} = " << limit << std::endl;
575  return true;
576 }
void applyExpectedQuantile(RooStats::HypoTestResult &hcres)
Definition: HybridNew.cc:1068
std::auto_ptr< RooStats::TestStatistic > qvar
Definition: HybridNew.h:84
tuple result
Definition: query.py:137
static std::string testStat_
Definition: HybridNew.h:44
RooStats::HypoTestResult * readToysFromFile(const RooAbsCollection &rVals)
Definition: HybridNew.cc:1293
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
RooStats::ModelConfig modelConfig_bonly
Definition: HybridNew.h:83
tuple cout
Definition: gather_cfg.py:121
static bool expectedFromGrid_
Definition: HybridNew.h:52
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
static bool readHybridResults_
Definition: HybridNew.h:50
static RooArgSet rValues_
Definition: HybridNew.h:48
T w() const
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
Definition: HybridNew.cc:639
void HybridNew::setupPOI ( RooStats::ModelConfig *  mc_s)
private

Definition at line 220 of file HybridNew.cc.

References NULL, rValue_, and rValues_.

Referenced by run().

220  {
221  RooArgSet POI(*mc_s->GetParametersOfInterest());
222  if (rValue_.find("=") == std::string::npos) {
223  if (POI.getSize() != 1) throw std::invalid_argument("Error: the argument to --singlePoint or --signalForSignificance is a single value, but there are multiple POIs");
224  POI.snapshot(rValues_);
225  errno = 0; // check for errors in str->float conversion
226  ((RooRealVar*)rValues_.first())->setVal(strtod(rValue_.c_str(),NULL));
227  if (errno != 0) std::invalid_argument("Error: the argument to --singlePoint or --signalForSignificance is not a valid number.");
228  } else {
229  std::string::size_type eqidx = 0, colidx = 0, colidx2;
230  do {
231  eqidx = rValue_.find("=", colidx);
232  colidx2 = rValue_.find(",", colidx+1);
233  if (eqidx == std::string::npos || (colidx2 != std::string::npos && colidx2 < eqidx)) {
234  throw std::invalid_argument("Error: the argument to --singlePoint or --signalForSignificance is not in the forms 'value' or 'name1=value1,name2=value2,...'\n");
235  }
236  std::string poiName = rValue_.substr(colidx, eqidx);
237  std::string poiVal = rValue_.substr(eqidx+1, (colidx2 == std::string::npos ? std::string::npos : colidx2 - eqidx - 1));
238  RooAbsArg *poi = POI.find(poiName.c_str());
239  if (poi == 0) throw std::invalid_argument("Error: unknown parameter '"+poiName+"' passed to --singlePoint or --signalForSignificance.");
240  rValues_.addClone(*poi);
241  errno = 0;
242  rValues_.setRealValue(poi->GetName(), strtod(poiVal.c_str(),NULL));
243  if (errno != 0) throw std::invalid_argument("Error: invalid value '"+poiVal+"' for parameter '"+poiName+"' passed to --singlePoint or --signalForSignificance.");
244  } while (colidx2 != std::string::npos);
245  if (rValues_.getSize() != POI.getSize()) {
246  throw std::invalid_argument("Error: not all parameters of interest specified in --singlePoint or --signalForSignificance");
247  }
248  }
249 }
static std::string rValue_
Definition: HybridNew.h:47
#define NULL
Definition: scimark2.h:8
uint16_t size_type
static RooArgSet rValues_
Definition: HybridNew.h:48
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().

1379  {
1380  typedef std::map<double, RooStats::HypoTestResult *>::iterator point;
1381  if (!smart) {
1382  for (point it = grid_.begin(), ed = grid_.end(); it != ed; ++it) {
1383  it->second->ResetBit(1);
1384  updateGridPoint(w, mc_s, mc_b, data, it);
1385  }
1386  } else {
1387  typedef std::pair<double,double> CLs_t;
1388  std::vector<point> points; points.reserve(grid_.size());
1389  std::vector<CLs_t> values; values.reserve(grid_.size());
1390  for (point it = grid_.begin(), ed = grid_.end(); it != ed; ++it) { points.push_back(it); values.push_back(CLs_t(-99, -99)); }
1391  int iMin = 0, iMax = points.size()-1;
1392  while (iMax-iMin > 3) {
1393  if (verbose > 1) std::cout << "Bisecting range [" << iMin << ", " << iMax << "]" << std::endl;
1394  int iMid = (iMin+iMax)/2;
1395  CLs_t clsMid = values[iMid] = updateGridPoint(w, mc_s, mc_b, data, points[iMid]);
1396  if (verbose > 1) std::cout << " Midpoint " << iMid << " value " << clsMid.first << " +/- " << clsMid.second << std::endl;
1397  if (clsMid.first - 3*max(clsMid.second,0.01) > clsTarget_) {
1398  if (verbose > 1) std::cout << " Replacing Min" << std::endl;
1399  iMin = iMid; continue;
1400  } else if (clsMid.first + 3*max(clsMid.second,0.01) < clsTarget_) {
1401  if (verbose > 1) std::cout << " Replacing Max" << std::endl;
1402  iMax = iMid; continue;
1403  } else {
1404  if (verbose > 1) std::cout << " Tightening Range" << std::endl;
1405  while (iMin < iMid-1) {
1406  int iLo = (iMin+iMid)/2;
1407  CLs_t clsLo = values[iLo] = updateGridPoint(w, mc_s, mc_b, data, points[iLo]);
1408  if (verbose > 1) std::cout << " Lowpoint " << iLo << " value " << clsLo.first << " +/- " << clsLo.second << std::endl;
1409  if (clsLo.first - 3*max(clsLo.second,0.01) > clsTarget_) iMin = iLo;
1410  else break;
1411  }
1412  while (iMax > iMid+1) {
1413  int iHi = (iMax+iMid)/2;
1414  CLs_t clsHi = values[iHi] = updateGridPoint(w, mc_s, mc_b, data, points[iHi]);
1415  if (verbose > 1) std::cout << " Highpoint " << iHi << " value " << clsHi.first << " +/- " << clsHi.second << std::endl;
1416  if (clsHi.first + 3*max(clsHi.second,0.01) < clsTarget_) iMax = iHi;
1417  else break;
1418  }
1419  break;
1420  }
1421  }
1422  if (verbose > 1) std::cout << "Final range [" << iMin << ", " << iMax << "]" << std::endl;
1423  for (int i = 0; i < iMin; ++i) {
1424  points[i]->second->SetBit(1);
1425  if (verbose > 1) std::cout << " Will not use point " << i << " (r " << points[i]->first << ")" << std::endl;
1426  }
1427  for (int i = iMin; i <= iMax; ++i) {
1428  points[i]->second->ResetBit(1);
1429  if (values[i].first < -2) {
1430  if (verbose > 1) std::cout << " Updaing point " << i << " (r " << points[i]->first << ")" << std::endl;
1431  updateGridPoint(w, mc_s, mc_b, data, points[i]);
1432  }
1433  else if (verbose > 1) std::cout << " Point " << i << " (r " << points[i]->first << ") was already updated during search." << std::endl;
1434  }
1435  for (int i = iMax+1, n = points.size(); i < n; ++i) {
1436  points[i]->second->SetBit(1);
1437  if (verbose > 1) std::cout << " Will not use point " << i << " (r " << points[i]->first << ")" << std::endl;
1438  }
1439  }
1440 }
int i
Definition: DBlmapReader.cc:9
const T & max(const T &a, const T &b)
std::map< double, RooStats::HypoTestResult * > grid_
Definition: HybridNew.h:119
bool first
Definition: L1TdeRCT.cc:94
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
std::pair< double, double > updateGridPoint(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, std::map< double, RooStats::HypoTestResult * >::iterator point)
Definition: HybridNew.cc:1462
tuple cout
Definition: gather_cfg.py:121
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
T w() const
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().

1441  {
1442  typedef std::map<double, RooStats::HypoTestResult *>::iterator point;
1443  std::vector<Double_t> rToUpdate; std::vector<point> pointToUpdate;
1444  for (point it = grid_.begin(), ed = grid_.end(); it != ed; ++it) {
1445  it->second->ResetBit(1);
1446  if (it->first == 0 || fabs(it->second->GetTestStatisticData()) <= 2*minimizerTolerance_) {
1447  rToUpdate.push_back(it->first);
1448  pointToUpdate.push_back(it);
1449  }
1450  }
1451  if (verbose > 0) std::cout << "A total of " << rToUpdate.size() << " points will be updated." << std::endl;
1452  Setup setup;
1453  std::auto_ptr<RooStats::HybridCalculator> hc = create(w, mc_s, mc_b, data, rToUpdate.back(), setup);
1454  RooArgSet nullPOI(*setup.modelConfig_bonly.GetSnapshot());
1455  std::vector<Double_t> qVals = ((ProfiledLikelihoodTestStatOpt&)(*setup.qvar)).Evaluate(data, nullPOI, rToUpdate);
1456  for (int i = 0, n = rToUpdate.size(); i < n; ++i) {
1457  pointToUpdate[i]->second->SetTestStatisticData(qVals[i] - EPS);
1458  }
1459  if (verbose > 0) std::cout << "All points have been updated." << std::endl;
1460 }
int i
Definition: DBlmapReader.cc:9
static float minimizerTolerance_
Definition: HybridNew.h:63
#define EPS
Definition: HybridNew.cc:87
std::map< double, RooStats::HypoTestResult * > grid_
Definition: HybridNew.h:119
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
tuple cout
Definition: gather_cfg.py:121
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
T w() const
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
Definition: HybridNew.cc:639
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().

1462  {
1463  typedef std::pair<double,double> CLs_t;
1464  bool isProfile = (testStat_ == "LHC" || testStat_ == "LHCFC" || testStat_ == "Profile");
1465  if (point->first == 0 && CLs_) return std::pair<double,double>(1,0);
1466  RooArgSet poi(*mc_s->GetParametersOfInterest());
1467  RooRealVar *r = dynamic_cast<RooRealVar *>(poi.first());
1468  if (expectedFromGrid_) {
1469  applyExpectedQuantile(*point->second);
1470  point->second->SetTestStatisticData(point->second->GetTestStatisticData() + (isProfile ? EPS : EPS));
1471  } else {
1472  Setup setup;
1473  std::auto_ptr<RooStats::HybridCalculator> hc = create(w, mc_s, mc_b, data, point->first, setup);
1474  RooArgSet nullPOI(*setup.modelConfig_bonly.GetSnapshot());
1475  if (isProfile) nullPOI.setRealValue(r->GetName(), point->first);
1476  double testStat = setup.qvar->Evaluate(data, nullPOI);
1477  point->second->SetTestStatisticData(testStat + (isProfile ? EPS : EPS));
1478  }
1479  if (verbose > 1) {
1480  std::cout << "At " << r->GetName() << " = " << point->first << ":\n" <<
1481  "\tCLs = " << point->second->CLs() << " +/- " << point->second->CLsError() << "\n" <<
1482  "\tCLb = " << point->second->CLb() << " +/- " << point->second->CLbError() << "\n" <<
1483  "\tCLsplusb = " << point->second->CLsplusb() << " +/- " << point->second->CLsplusbError() << "\n" <<
1484  std::endl;
1485  }
1486 
1487  return eval(*point->second, point->first);
1488 }
#define EPS
Definition: HybridNew.cc:87
void applyExpectedQuantile(RooStats::HypoTestResult &hcres)
Definition: HybridNew.cc:1068
static bool CLs_
Definition: HybridNew.h:43
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
Definition: HybridNew.cc:587
static std::string testStat_
Definition: HybridNew.h:44
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
tuple cout
Definition: gather_cfg.py:121
static bool expectedFromGrid_
Definition: HybridNew.h:52
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
T w() const
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
Definition: HybridNew.cc:639
void HybridNew::useGrid ( )
private

Definition at line 1489 of file HybridNew.cc.

References CLs_, eval(), grid_, i, limitPlot_, and n.

Referenced by runLimit().

1489  {
1490  typedef std::pair<double,double> CLs_t;
1491  int i = 0, n = 0;
1492  limitPlot_.reset(new TGraphErrors(1));
1493  for (std::map<double, RooStats::HypoTestResult *>::iterator itg = grid_.begin(), edg = grid_.end(); itg != edg; ++itg, ++i) {
1494  if (itg->second->TestBit(1)) continue;
1495  CLs_t val(1,0);
1496  if (CLs_) {
1497  if (itg->first > 0) val = eval(*itg->second, itg->first);
1498  } else {
1499  val = eval(*itg->second, itg->first);
1500  }
1501  if (val.first == -1) continue;
1502  if (val.second == 0 && (val.first != 1 && val.first != 0)) continue;
1503  limitPlot_->Set(n+1);
1504  limitPlot_->SetPoint( n, itg->first, val.first);
1505  limitPlot_->SetPointError(n, 0, val.second);
1506  n++;
1507  }
1508 }
int i
Definition: DBlmapReader.cc:9
static bool CLs_
Definition: HybridNew.h:43
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
Definition: HybridNew.cc:587
std::map< double, RooStats::HypoTestResult * > grid_
Definition: HybridNew.h:119
std::auto_ptr< TGraphErrors > limitPlot_
Definition: HybridNew.h:73
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().

185  {
186  if (fork_ > 1) nToys_ /= fork_; // makes more sense
187  if (rule_ == "CLs") {
188  CLs_ = true;
189  } else if (rule_ == "CLsplusb") {
190  CLs_ = false;
191  } else if (rule_ == "FC") {
192  CLs_ = false;
193  } else {
194  throw std::invalid_argument("HybridNew: Rule must be either 'CLs' or 'CLsplusb'");
195  }
196  if (testStat_ == "SimpleLikelihoodRatio" || testStat_ == "SLR" ) { testStat_ = "LEP"; }
197  if (testStat_ == "RatioOfProfiledLikelihoods" || testStat_ == "ROPL") { testStat_ = "TEV"; }
198  if (testStat_ == "ProfileLikelihood" || testStat_ == "PL") { testStat_ = "Profile"; }
199  if (testStat_ == "ModifiedProfileLikelihood" || testStat_ == "MPL") { testStat_ = "LHC"; }
200  if (testStat_ == "SignFlipProfileLikelihood" || testStat_ == "SFPL") { testStat_ = "LHCFC"; }
201  if (testStat_ == "Atlas") { testStat_ = "LHC"; std::cout << "Note: the Atlas test statistics is now known as LHC test statistics.\n" << std::endl; }
202  if (testStat_ != "LEP" && testStat_ != "TEV" && testStat_ != "LHC" && testStat_ != "LHCFC" && testStat_ != "Profile" && testStat_ != "MLZ") {
203  throw std::invalid_argument("HybridNew: Test statistics should be one of 'LEP' or 'TEV' or 'LHC' (previously known as 'Atlas') or 'Profile'");
204  }
205  if (verbose) {
206  if (testStat_ == "LEP") std::cout << ">>> using the Simple Likelihood Ratio test statistics (Q_LEP)" << std::endl;
207  if (testStat_ == "TEV") std::cout << ">>> using the Ratio of Profiled Likelihoods test statistics (Q_TEV)" << std::endl;
208  if (testStat_ == "LHC") std::cout << ">>> using the Profile Likelihood test statistics modified for upper limits (Q_LHC)" << std::endl;
209  if (testStat_ == "LHCFC") std::cout << ">>> using the Profile Likelihood test statistics modified for upper limits and Feldman-Cousins (Q_LHCFC)" << std::endl;
210  if (testStat_ == "Profile") std::cout << ">>> using the Profile Likelihood test statistics not modified for upper limits (Q_Profile)" << std::endl;
211  if (testStat_ == "MLZ") std::cout << ">>> using the Maximum likelihood estimator of the signal strength as test statistics" << std::endl;
212  }
214  // If not generating toys, don't need to fit nuisance parameters
215  fitNuisances_ = false;
216  }
217  if (reportPVal_ && workingMode_ != MakeSignificance) throw std::invalid_argument("HybridNew: option --pvalue must go together with --significance");
218 }
static std::string rule_
Definition: HybridNew.h:44
static unsigned int fork_
Definition: HybridNew.h:58
static WorkingMode workingMode_
Definition: HybridNew.h:40
static bool reportPVal_
Definition: HybridNew.h:45
static bool CLs_
Definition: HybridNew.h:43
static unsigned int nToys_
Definition: HybridNew.h:41
static std::string testStat_
Definition: HybridNew.h:44
tuple cout
Definition: gather_cfg.py:121
static bool fitNuisances_
Definition: HybridNew.h:46
static bool readHybridResults_
Definition: HybridNew.h:50

Member Data Documentation

std::string HybridNew::algo_ = "logSecant"
staticprivate

Definition at line 60 of file HybridNew.h.

Referenced by HybridNew(), and runLimit().

bool HybridNew::CLs_ = false
staticprivate
double HybridNew::clsAccuracy_ = 0.005
staticprivate

Definition at line 42 of file HybridNew.h.

Referenced by eval(), HybridNew(), and runSinglePoint().

bool HybridNew::clsQuantiles_ = true
staticprivate

Definition at line 52 of file HybridNew.h.

Referenced by applyExpectedQuantile(), and HybridNew().

bool HybridNew::expectedFromGrid_ = false
staticprivate
bool HybridNew::fitNuisances_ = false
staticprivate

Definition at line 46 of file HybridNew.h.

Referenced by applyOptions(), create(), HybridNew(), and validateOptions().

unsigned int HybridNew::fork_ = 1
staticprivate

Definition at line 58 of file HybridNew.h.

Referenced by evalGeneric(), evalWithFork(), HybridNew(), and validateOptions().

bool HybridNew::fullBToys_ = false
staticprivate

Definition at line 54 of file HybridNew.h.

Referenced by applyOptions(), and create().

bool HybridNew::fullGrid_ = false
staticprivate

Definition at line 55 of file HybridNew.h.

Referenced by applyOptions(), and runLimit().

bool HybridNew::genGlobalObs_ = false
staticprivate

Definition at line 46 of file HybridNew.h.

Referenced by applyOptions(), create(), and HybridNew().

bool HybridNew::genNuisances_ = true
staticprivate

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_ = ""
staticprivate

Definition at line 51 of file HybridNew.h.

Referenced by HybridNew(), and runLimit().

bool HybridNew::importanceSamplingAlt_ = false
staticprivate

Definition at line 59 of file HybridNew.h.

Referenced by create().

bool HybridNew::importanceSamplingNull_ = false
staticprivate

Definition at line 59 of file HybridNew.h.

Referenced by create().

double HybridNew::interpAccuracy_ = 0.2
staticprivate

Definition at line 42 of file HybridNew.h.

Referenced by HybridNew(), and runLimit().

unsigned int HybridNew::iterations_ = 1
staticprivate

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"
staticprivate

Definition at line 62 of file HybridNew.h.

Referenced by HybridNew(), and run().

float HybridNew::minimizerTolerance_ = 1e-2
staticprivate

Definition at line 63 of file HybridNew.h.

Referenced by eval(), HybridNew(), run(), and updateGridDataFC().

unsigned int HybridNew::nCpu_ = 0
staticprivate

Definition at line 58 of file HybridNew.h.

Referenced by create(), and HybridNew().

bool HybridNew::newToyMCSampler_ = true
staticprivate

Definition at line 67 of file HybridNew.h.

Referenced by create(), and HybridNew().

bool HybridNew::noUpdateGrid_ = false
staticprivate

Definition at line 57 of file HybridNew.h.

Referenced by applyOptions(), eval(), and runLimit().

unsigned int HybridNew::nToys_ = 500
staticprivate

Definition at line 41 of file HybridNew.h.

Referenced by create(), eval(), HybridNew(), and validateOptions().

bool HybridNew::optimizeProductPdf_ = true
staticprivate

Definition at line 65 of file HybridNew.h.

Referenced by create(), and HybridNew().

bool HybridNew::optimizeTestStatistics_ = true
staticprivate

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_
staticprivate

Definition at line 61 of file HybridNew.h.

Referenced by eval(), HybridNew(), readToysFromFile(), and runLimit().

float HybridNew::quantileForExpectedFromGrid_ = 0.5
staticprivate
double HybridNew::rAbsAccuracy_ = 0.1
staticprivate

Definition at line 42 of file HybridNew.h.

Referenced by HybridNew(), and runLimit().

bool HybridNew::readHybridResults_ = false
staticprivate
bool HybridNew::reportPVal_ = false
staticprivate

Definition at line 45 of file HybridNew.h.

Referenced by applyOptions(), runSignificance(), and validateOptions().

bool HybridNew::rMaxSet_ = false
staticprivate

Definition at line 69 of file HybridNew.h.

Referenced by applyOptions(), and runLimit().

bool HybridNew::rMinSet_ = false
staticprivate

Definition at line 68 of file HybridNew.h.

Referenced by applyOptions(), and runLimit().

double HybridNew::rRelAccuracy_ = 0.05
staticprivate

Definition at line 42 of file HybridNew.h.

Referenced by HybridNew(), and runLimit().

std::string HybridNew::rule_ = "CLs"
staticprivate

Definition at line 44 of file HybridNew.h.

Referenced by eval(), HybridNew(), runLimit(), and validateOptions().

std::string HybridNew::rValue_ = "1.0"
staticprivate

Definition at line 47 of file HybridNew.h.

Referenced by applyOptions(), HybridNew(), and setupPOI().

RooArgSet HybridNew::rValues_
staticprivate

Definition at line 48 of file HybridNew.h.

Referenced by readGrid(), run(), runSignificance(), runSinglePoint(), runTestStatistics(), and setupPOI().

bool HybridNew::saveGrid_ = false
staticprivate

Definition at line 56 of file HybridNew.h.

Referenced by applyOptions(), and runLimit().

bool HybridNew::saveHybridResult_ = false
staticprivate

Definition at line 50 of file HybridNew.h.

Referenced by applyOptions(), eval(), and runSignificance().

std::string HybridNew::testStat_ = "LEP"
staticprivate
HybridNew::WorkingMode HybridNew::workingMode_ = MakeLimit
staticprivate