Logo ROOT   master
Reference Guide
HypoTestInverter.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
3 // Contributions: Giovanni Petrucciani and Annapaola Decosa
4 /*************************************************************************
5  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class RooStats::HypoTestInverter
13  \ingroup Roostats
14 
15 A class for performing a hypothesis test inversion by scanning
16 the hypothesis test results of a HypoTestCalculator for various values of the
17 parameter of interest. By looking at the confidence level curve of the result, an
18 upper limit can be derived by computing the intersection of the confidence level curve with the desired confidence level.
19 The class implements the RooStats::IntervalCalculator interface, and returns a
20 RooStats::HypoTestInverterResult. The result is a SimpleInterval, which
21 via the method UpperLimit() returns to the user the upper limit value.
22 
23 ## Scanning options
24 The HypoTestInverter implements various options for performing the scan.
25 - HypoTestInverter::RunFixedScan will scan the parameter of interest using a fixed grid.
26 - HypoTestInverter::SetAutoScan will perform an automatic scan to find
27 optimally the curve. It will stop when the desired precision is obtained.
28 - HypoTestInverter::RunOnePoint computes the confidence level at a given point.
29 
30 ### CLs presciption
31 The class can scan the CLs+b values or alternatively CLs. For the latter,
32 call HypoTestInverter::UseCLs().
33 */
34 
35 // include other header files
36 
37 #include "RooAbsData.h"
38 #
39 #include "TMath.h"
40 
41 #include "RooStats/HybridResult.h"
42 
44 #include <cassert>
45 #include <cmath>
46 
47 #include "TF1.h"
48 #include "TFile.h"
49 #include "TH1.h"
50 #include "TLine.h"
51 #include "TCanvas.h"
52 #include "TGraphErrors.h"
53 #include "RooRealVar.h"
54 #include "RooArgSet.h"
55 #include "RooAbsPdf.h"
56 #include "RooRandom.h"
57 #include "RooConstVar.h"
58 #include "RooMsgService.h"
59 #include "RooStats/ModelConfig.h"
66 #include "RooStats/ToyMCSampler.h"
67 #include "RooStats/HypoTestPlot.h"
69 
70 #include "RooStats/ProofConfig.h"
71 
73 
74 using namespace RooStats;
75 using namespace std;
76 
77 // static variable definitions
78 double HypoTestInverter::fgCLAccuracy = 0.005;
79 unsigned int HypoTestInverter::fgNToys = 500;
80 
83 std::string HypoTestInverter::fgAlgo = "logSecant";
84 
86 
87 // helper class to wrap the functionality of the various HypoTestCalculators
88 
89 template<class HypoTestType>
90 struct HypoTestWrapper {
91 
92  static void SetToys(HypoTestType * h, int toyNull, int toyAlt) { h->SetToys(toyNull,toyAlt); }
93 
94 };
95 
96 ////////////////////////////////////////////////////////////////////////////////
97 /// set flag to close proof for every new run
98 
100  fgCloseProof = flag;
101 }
102 
103 ////////////////////////////////////////////////////////////////////////////////
104 /// get the variable to scan
105 /// try first with null model if not go to alternate model
106 
108 
109  RooRealVar * varToScan = 0;
110  const ModelConfig * mc = hc.GetNullModel();
111  if (mc) {
112  const RooArgSet * poi = mc->GetParametersOfInterest();
113  if (poi) varToScan = dynamic_cast<RooRealVar*> (poi->first() );
114  }
115  if (!varToScan) {
116  mc = hc.GetAlternateModel();
117  if (mc) {
118  const RooArgSet * poi = mc->GetParametersOfInterest();
119  if (poi) varToScan = dynamic_cast<RooRealVar*> (poi->first() );
120  }
121  }
122  return varToScan;
123 }
124 
125 ////////////////////////////////////////////////////////////////////////////////
126 /// check the model given the given hypotestcalculator
127 
129  const ModelConfig * modelSB = hc.GetNullModel();
130  const ModelConfig * modelB = hc.GetAlternateModel();
131  if (!modelSB || ! modelB)
132  oocoutF((TObject*)0,InputArguments) << "HypoTestInverter - model are not existing" << std::endl;
133  assert(modelSB && modelB);
134 
135  oocoutI((TObject*)0,InputArguments) << "HypoTestInverter ---- Input models: \n"
136  << "\t\t using as S+B (null) model : "
137  << modelSB->GetName() << "\n"
138  << "\t\t using as B (alternate) model : "
139  << modelB->GetName() << "\n" << std::endl;
140 
141  // check if scanVariable is included in B model pdf
142  RooAbsPdf * bPdf = modelB->GetPdf();
143  const RooArgSet * bObs = modelB->GetObservables();
144  if (!bPdf || !bObs) {
145  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - B model has no pdf or observables defined" << std::endl;
146  return;
147  }
148  RooArgSet * bParams = bPdf->getParameters(*bObs);
149  if (!bParams) {
150  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - pdf of B model has no parameters" << std::endl;
151  return;
152  }
153  if (bParams->find(scanVariable.GetName() ) ) {
154  const RooArgSet * poiB = modelB->GetSnapshot();
155  if (!poiB || !poiB->find(scanVariable.GetName()) ||
156  ( (RooRealVar*) poiB->find(scanVariable.GetName()) )->getVal() != 0 )
157  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter - using a B model with POI "
158  << scanVariable.GetName() << " not equal to zero "
159  << " user must check input model configurations " << endl;
160  if (poiB) delete poiB;
161  }
162  delete bParams;
163 }
164 
165 ////////////////////////////////////////////////////////////////////////////////
166 /// default constructor (doesn't do anything)
167 
169  fTotalToysRun(0),
170  fMaxToys(0),
171  fCalculator0(0),
172  fScannedVariable(0),
173  fResults(0),
174  fUseCLs(false),
175  fScanLog(false),
176  fSize(0),
177  fVerbose(0),
178  fCalcType(kUndefined),
179  fNBins(0), fXmin(1), fXmax(1),
180  fNumErr(0)
181 {
182 }
183 
184 ////////////////////////////////////////////////////////////////////////////////
185 /// Constructor from a HypoTestCalculatorGeneric
186 /// The HypoTest calculator must be a FrequentistCalculator or HybridCalculator type
187 /// Other type of calculators are not supported.
188 /// The calculator must be created before by using the S+B model for the null and
189 /// the B model for the alt
190 /// If no variable to scan are given they are assumed to be the first variable
191 /// from the parameter of interests of the null model
192 
194  RooRealVar* scannedVariable, double size ) :
195  fTotalToysRun(0),
196  fMaxToys(0),
197  fCalculator0(0),
198  fScannedVariable(scannedVariable),
199  fResults(0),
200  fUseCLs(false),
201  fScanLog(false),
202  fSize(size),
203  fVerbose(0),
204  fCalcType(kUndefined),
205  fNBins(0), fXmin(1), fXmax(1),
206  fNumErr(0)
207 {
208 
209  if (!fScannedVariable) {
211  }
212  if (!fScannedVariable)
213  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
214  else
216 
217  HybridCalculator * hybCalc = dynamic_cast<HybridCalculator*>(&hc);
218  if (hybCalc) {
219  fCalcType = kHybrid;
220  fCalculator0 = hybCalc;
221  return;
222  }
223  FrequentistCalculator * freqCalc = dynamic_cast<FrequentistCalculator*>(&hc);
224  if (freqCalc) {
226  fCalculator0 = freqCalc;
227  return;
228  }
229  AsymptoticCalculator * asymCalc = dynamic_cast<AsymptoticCalculator*>(&hc);
230  if (asymCalc) {
232  fCalculator0 = asymCalc;
233  return;
234  }
235  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Type of hypotest calculator is not supported " <<std::endl;
236  fCalculator0 = &hc;
237 }
238 
239 ////////////////////////////////////////////////////////////////////////////////
240 /// Constructor from a reference to a HybridCalculator
241 /// The calculator must be created before by using the S+B model for the null and
242 /// the B model for the alt
243 /// If no variable to scan are given they are assumed to be the first variable
244 /// from the parameter of interests of the null model
245 
247  RooRealVar* scannedVariable, double size ) :
248  fTotalToysRun(0),
249  fMaxToys(0),
250  fCalculator0(&hc),
251  fScannedVariable(scannedVariable),
252  fResults(0),
253  fUseCLs(false),
254  fScanLog(false),
255  fSize(size),
256  fVerbose(0),
257  fCalcType(kHybrid),
258  fNBins(0), fXmin(1), fXmax(1),
259  fNumErr(0)
260 {
261 
262  if (!fScannedVariable) {
264  }
265  if (!fScannedVariable)
266  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
267  else
269 
270 }
271 
272 ////////////////////////////////////////////////////////////////////////////////
273 /// Constructor from a reference to a FrequentistCalculator
274 /// The calculator must be created before by using the S+B model for the null and
275 /// the B model for the alt
276 /// If no variable to scan are given they are assumed to be the first variable
277 /// from the parameter of interests of the null model
278 
280  RooRealVar* scannedVariable, double size ) :
281  fTotalToysRun(0),
282  fMaxToys(0),
283  fCalculator0(&hc),
284  fScannedVariable(scannedVariable),
285  fResults(0),
286  fUseCLs(false),
287  fScanLog(false),
288  fSize(size),
289  fVerbose(0),
290  fCalcType(kFrequentist),
291  fNBins(0), fXmin(1), fXmax(1),
292  fNumErr(0)
293 {
294 
295  if (!fScannedVariable) {
297  }
298  if (!fScannedVariable)
299  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
300  else
302 }
303 
304 ////////////////////////////////////////////////////////////////////////////////
305 /// Constructor from a reference to a AsymptoticCalculator
306 /// The calculator must be created before by using the S+B model for the null and
307 /// the B model for the alt
308 /// If no variable to scan are given they are assumed to be the first variable
309 /// from the parameter of interests of the null model
310 
312  RooRealVar* scannedVariable, double size ) :
313  fTotalToysRun(0),
314  fMaxToys(0),
315  fCalculator0(&hc),
316  fScannedVariable(scannedVariable),
317  fResults(0),
318  fUseCLs(false),
319  fScanLog(false),
320  fSize(size),
321  fVerbose(0),
322  fCalcType(kAsymptotic),
323  fNBins(0), fXmin(1), fXmax(1),
324  fNumErr(0)
325 {
326 
327  if (!fScannedVariable) {
329  }
330  if (!fScannedVariable)
331  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
332  else
334 
335 }
336 
337 ////////////////////////////////////////////////////////////////////////////////
338 /// Constructor from a model for B model and a model for S+B.
339 /// An HypoTestCalculator (Hybrid of Frequentis) will be created using the
340 /// S+B model as the null and the B model as the alternate
341 /// If no variable to scan are given they are assumed to be the first variable
342 /// from the parameter of interests of the null model
343 
345  RooRealVar * scannedVariable, ECalculatorType type, double size) :
346  fTotalToysRun(0),
347  fMaxToys(0),
348  fCalculator0(0),
349  fScannedVariable(scannedVariable),
350  fResults(0),
351  fUseCLs(false),
352  fScanLog(false),
353  fSize(size),
354  fVerbose(0),
355  fCalcType(type),
356  fNBins(0), fXmin(1), fXmax(1),
357  fNumErr(0)
358 {
359  if(fCalcType==kFrequentist) fHC.reset(new FrequentistCalculator(data, bModel, sbModel));
360  if(fCalcType==kHybrid) fHC.reset( new HybridCalculator(data, bModel, sbModel)) ;
361  if(fCalcType==kAsymptotic) fHC.reset( new AsymptoticCalculator(data, bModel, sbModel));
362  fCalculator0 = fHC.get();
363  // get scanned variable
364  if (!fScannedVariable) {
366  }
367  if (!fScannedVariable)
368  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
369  else
371 
372 }
373 
374 ////////////////////////////////////////////////////////////////////////////////
375 /// copy-constructor
376 /// NOTE: this class does not copy the contained result and
377 /// the HypoTestCalculator, but only the pointers
378 /// It requires the original HTI to be alive
379 
382  fTotalToysRun(0),
383  fCalculator0(0), fScannedVariable(0), // add these for Coverity
384  fResults(0)
385 {
386  (*this) = rhs;
387 }
388 
389 ////////////////////////////////////////////////////////////////////////////////
390 /// assignment operator
391 /// NOTE: this class does not copy the contained result and
392 /// the HypoTestCalculator, but only the pointers
393 /// It requires the original HTI to be alive
394 
396  if (this == &rhs) return *this;
397  fTotalToysRun = 0;
398  fMaxToys = rhs.fMaxToys;
401  fUseCLs = rhs.fUseCLs;
402  fScanLog = rhs.fScanLog;
403  fSize = rhs.fSize;
404  fVerbose = rhs.fVerbose;
405  fCalcType = rhs.fCalcType;
406  fNBins = rhs.fNBins;
407  fXmin = rhs.fXmin;
408  fXmax = rhs.fXmax;
409  fNumErr = rhs.fNumErr;
410 
411  return *this;
412 }
413 
414 ////////////////////////////////////////////////////////////////////////////////
415 /// destructor (delete the HypoTestInverterResult)
416 
418 {
419  if (fResults) delete fResults;
420  fCalculator0 = 0;
421 }
422 
423 ////////////////////////////////////////////////////////////////////////////////
424 /// return the test statistic which is or will be used by the class
425 
427 {
430  }
431  else
432  return 0;
433 }
434 
435 ////////////////////////////////////////////////////////////////////////////////
436 /// set the test statistic to use
437 
439 {
442  return true;
443  }
444  else return false;
445 }
446 
447 ////////////////////////////////////////////////////////////////////////////////
448 /// delete contained result and graph
449 
451  if (fResults) delete fResults;
452  fResults = 0;
453  fLimitPlot.reset(nullptr);
454 }
455 
456 ////////////////////////////////////////////////////////////////////////////////
457 /// create a new HypoTestInverterResult to hold all computed results
458 
460  if (fResults == 0) {
461  TString results_name = "result_";
462  results_name += fScannedVariable->GetName();
464  TString title = "HypoTestInverter Result For ";
465  title += fScannedVariable->GetName();
466  fResults->SetTitle(title);
467  }
470  // check if one or two sided scan
471  if (fCalculator0) {
472  // if asymptotic calculator
474  if (ac)
476  else {
477  // in case of the other calculators
479  if (sampler) {
480  ProfileLikelihoodTestStat * pl = dynamic_cast<ProfileLikelihoodTestStat*>(sampler->GetTestStatistic());
481  if (pl && pl->IsTwoSided() ) fResults->fIsTwoSided = true;
482  }
483  }
484  }
485 }
486 
487 ////////////////////////////////////////////////////////////////////////////////
488 /// Run a fixed scan or the automatic scan depending on the configuration.
489 /// Return if needed a copy of the result object which will be managed by the user.
490 
492 
493  // if having a result with at least one point return it
494  if (fResults && fResults->ArraySize() >= 1) {
495  oocoutI((TObject*)0,Eval) << "HypoTestInverter::GetInterval - return an already existing interval " << std::endl;
496  return (HypoTestInverterResult*)(fResults->Clone());
497  }
498 
499  if (fNBins > 0) {
500  oocoutI((TObject*)0,Eval) << "HypoTestInverter::GetInterval - run a fixed scan" << std::endl;
501  bool ret = RunFixedScan(fNBins, fXmin, fXmax, fScanLog);
502  if (!ret)
503  oocoutE((TObject*)0,Eval) << "HypoTestInverter::GetInterval - error running a fixed scan " << std::endl;
504  }
505  else {
506  oocoutI((TObject*)0,Eval) << "HypoTestInverter::GetInterval - run an automatic scan" << std::endl;
507  double limit(0),err(0);
508  bool ret = RunLimit(limit,err);
509  if (!ret)
510  oocoutE((TObject*)0,Eval) << "HypoTestInverter::GetInterval - error running an auto scan " << std::endl;
511  }
512 
514 
515  return (HypoTestInverterResult*) (fResults->Clone());
516 }
517 
518 ////////////////////////////////////////////////////////////////////////////////
519 /// Run the Hypothesis test at a previous configured point
520 /// (internal function called by RunOnePoint)
521 
522 HypoTestResult * HypoTestInverter::Eval(HypoTestCalculatorGeneric &hc, bool adaptive, double clsTarget) const {
523  //for debug
524  // std::cout << ">>>>>>>>>>> " << std::endl;
525  // std::cout << "alternate model " << std::endl;
526  // hc.GetAlternateModel()->GetNuisanceParameters()->Print("V");
527  // hc.GetAlternateModel()->GetParametersOfInterest()->Print("V");
528  // std::cout << "Null model " << std::endl;
529  // hc.GetNullModel()->GetNuisanceParameters()->Print("V");
530  // hc.GetNullModel()->GetParametersOfInterest()->Print("V");
531  // std::cout << "<<<<<<<<<<<<<<< " << std::endl;
532 
533  // run the hypothesis test
534  HypoTestResult * hcResult = hc.GetHypoTest();
535  if (hcResult == 0) {
536  oocoutE((TObject*)0,Eval) << "HypoTestInverter::Eval - HypoTest failed" << std::endl;
537  return hcResult;
538  }
539 
540  // since the b model is the alt need to set the flag
541  hcResult->SetBackgroundAsAlt(true);
542 
543 
544  // bool flipPvalue = false;
545  // if (flipPValues)
546  // hcResult->SetPValueIsRightTail(!hcResult->GetPValueIsRightTail());
547 
548  // adjust for some numerical error in discrete models and == is not anymore
549  if (hcResult->GetPValueIsRightTail() )
550  hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()-fNumErr); // issue with < vs <= in discrete models
551  else
552  hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()+fNumErr); // issue with < vs <= in discrete models
553 
554  double clsMid = (fUseCLs ? hcResult->CLs() : hcResult->CLsplusb());
555  double clsMidErr = (fUseCLs ? hcResult->CLsError() : hcResult->CLsplusbError());
556 
557  //if (fVerbose) std::cout << (fUseCLs ? "\tCLs = " : "\tCLsplusb = ") << clsMid << " +/- " << clsMidErr << std::endl;
558 
559  if (adaptive) {
560 
561  if (fCalcType == kHybrid) HypoTestWrapper<HybridCalculator>::SetToys((HybridCalculator*)&hc, fUseCLs ? fgNToys : 1, 4*fgNToys);
562  if (fCalcType == kFrequentist) HypoTestWrapper<FrequentistCalculator>::SetToys((FrequentistCalculator*)&hc, fUseCLs ? fgNToys : 1, 4*fgNToys);
563 
564  while (clsMidErr >= fgCLAccuracy && (clsTarget == -1 || fabs(clsMid-clsTarget) < 3*clsMidErr) ) {
565  std::unique_ptr<HypoTestResult> more(hc.GetHypoTest());
566 
567  // if (flipPValues)
568  // more->SetPValueIsRightTail(!more->GetPValueIsRightTail());
569 
570  hcResult->Append(more.get());
571  clsMid = (fUseCLs ? hcResult->CLs() : hcResult->CLsplusb());
572  clsMidErr = (fUseCLs ? hcResult->CLsError() : hcResult->CLsplusbError());
573  if (fVerbose) std::cout << (fUseCLs ? "\tCLs = " : "\tCLsplusb = ") << clsMid << " +/- " << clsMidErr << std::endl;
574  }
575 
576  }
577  if (fVerbose ) {
578  oocoutP((TObject*)0,Eval) << "P values for " << fScannedVariable->GetName() << " = " <<
579  fScannedVariable->getVal() << "\n" <<
580  "\tCLs = " << hcResult->CLs() << " +/- " << hcResult->CLsError() << "\n" <<
581  "\tCLb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << "\n" <<
582  "\tCLsplusb = " << hcResult->CLsplusb() << " +/- " << hcResult->CLsplusbError() << "\n" <<
583  std::endl;
584  }
585 
586  if (fCalcType == kFrequentist || fCalcType == kHybrid) {
587  fTotalToysRun += (hcResult->GetAltDistribution()->GetSize() + hcResult->GetNullDistribution()->GetSize());
588 
589  // set sampling distribution name
590  TString nullDistName = TString::Format("%s_%s_%4.2f",hcResult->GetNullDistribution()->GetName(),
592  TString altDistName = TString::Format("%s_%s_%4.2f",hcResult->GetAltDistribution()->GetName(),
594 
595  hcResult->GetNullDistribution()->SetName(nullDistName);
596  hcResult->GetAltDistribution()->SetName(altDistName);
597  }
598 
599  return hcResult;
600 }
601 
602 ////////////////////////////////////////////////////////////////////////////////
603 /// Run a Fixed scan in npoints between min and max
604 
605 bool HypoTestInverter::RunFixedScan( int nBins, double xMin, double xMax, bool scanLog ) const
606 {
607 
608  CreateResults();
609  // interpolate the limits
610  fResults->fFittedLowerLimit = false;
611  fResults->fFittedUpperLimit = false;
612 
613  // safety checks
614  if ( nBins<=0 ) {
615  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - Please provide nBins>0\n";
616  return false;
617  }
618  if ( nBins==1 && xMin!=xMax ) {
619  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - nBins==1 -> I will run for xMin (" << xMin << ")\n";
620  }
621  if ( xMin==xMax && nBins>1 ) {
622  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - xMin==xMax -> I will enforce nBins==1\n";
623  nBins = 1;
624  }
625  if ( xMin>xMax ) {
626  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - Please provide xMin ("
627  << xMin << ") smaller than xMax (" << xMax << ")\n";
628  return false;
629  }
630 
631  if (xMin < fScannedVariable->getMin()) {
632  xMin = fScannedVariable->getMin();
633  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - xMin < lower bound, using xmin = "
634  << xMin << std::endl;
635  }
636  if (xMax > fScannedVariable->getMax()) {
637  xMax = fScannedVariable->getMax();
638  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - xMax > upper bound, using xmax = "
639  << xMax << std::endl;
640  }
641 
642  if (xMin <= 0. && scanLog) {
643  oocoutE((TObject*)nullptr, InputArguments) << "HypoTestInverter::RunFixedScan - cannot go in log steps if xMin <= 0" << std::endl;
644  return false;
645  }
646 
647  double thisX = xMin;
648  for (int i=0; i<nBins; i++) {
649 
650  if (i > 0) { // avoids case of nBins = 1
651  if (scanLog)
652  thisX = exp( log(xMin) + i*(log(xMax)-log(xMin))/(nBins-1) ); // scan in log x
653  else
654  thisX = xMin + i*(xMax-xMin)/(nBins-1); // linear scan in x
655  }
656 
657  bool status = RunOnePoint(thisX);
658 
659  // check if failed status
660  if ( status==false ) {
661  std::cout << "\t\tLoop interrupted because of failed status\n";
662  return false;
663  }
664  }
665 
666  return true;
667 }
668 
669 ////////////////////////////////////////////////////////////////////////////////
670 /// run only one point at the given POI value
671 
672 bool HypoTestInverter::RunOnePoint( double rVal, bool adaptive, double clTarget) const
673 {
674 
675  CreateResults();
676 
677  // check if rVal is in the range specified for fScannedVariable
678  if ( rVal < fScannedVariable->getMin() ) {
679  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunOnePoint - Out of range: using the lower bound "
681  << " on the scanned variable rather than " << rVal<< "\n";
682  rVal = fScannedVariable->getMin();
683  }
684  if ( rVal > fScannedVariable->getMax() ) {
685  // print a message when you have a significative difference since rval is computed
686  if ( rVal > fScannedVariable->getMax()*(1.+1.E-12) )
687  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunOnePoint - Out of range: using the upper bound "
689  << " on the scanned variable rather than " << rVal<< "\n";
690  rVal = fScannedVariable->getMax();
691  }
692 
693  // save old value
694  double oldValue = fScannedVariable->getVal();
695 
696  // evaluate hybrid calculator at a single point
697  fScannedVariable->setVal(rVal);
698  // need to set value of rval in hybridcalculator
699  // assume null model is S+B and alternate is B only
700  const ModelConfig * sbModel = fCalculator0->GetNullModel();
701  RooArgSet poi; poi.add(*sbModel->GetParametersOfInterest());
702  // set poi to right values
703  poi = RooArgSet(*fScannedVariable);
704  const_cast<ModelConfig*>(sbModel)->SetSnapshot(poi);
705 
706  if (fVerbose > 0)
707  oocoutP((TObject*)0,Eval) << "Running for " << fScannedVariable->GetName() << " = " << fScannedVariable->getVal() << endl;
708 
709  // compute the results
710  HypoTestResult* result = Eval(*fCalculator0,adaptive,clTarget);
711  if (!result) {
712  oocoutE((TObject*)0,Eval) << "HypoTestInverter - Error running point " << fScannedVariable->GetName() << " = " <<
713  fScannedVariable->getVal() << endl;
714  return false;
715  }
716  // in case of a dummy result
717  if (TMath::IsNaN(result->NullPValue() ) && TMath::IsNaN(result->AlternatePValue() ) ) {
718  oocoutW((TObject*)0,Eval) << "HypoTestInverter - Skip invalid result for point " << fScannedVariable->GetName() << " = " <<
719  fScannedVariable->getVal() << endl;
720  return true; // need to return true to avoid breaking the scan loop
721  }
722 
723  double lastXtested;
724  if ( fResults->ArraySize()!=0 ) lastXtested = fResults->GetXValue(fResults->ArraySize()-1);
725  else lastXtested = -999;
726 
727  if ( (std::abs(rVal) < 1 && TMath::AreEqualAbs(rVal, lastXtested,1.E-12) ) ||
728  (std::abs(rVal) >= 1 && TMath::AreEqualRel(rVal, lastXtested,1.E-12) ) ) {
729 
730  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunOnePoint - Merge with previous result for "
731  << fScannedVariable->GetName() << " = " << rVal << std::endl;
732  HypoTestResult* prevResult = fResults->GetResult(fResults->ArraySize()-1);
733  if (prevResult && prevResult->GetNullDistribution() && prevResult->GetAltDistribution()) {
734  prevResult->Append(result);
735  delete result; // we can delete the result
736  }
737  else {
738  // if it was empty we re-use it
739  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunOnePoint - replace previous empty result\n";
740  fResults->fYObjects.Remove( prevResult);
741  fResults->fYObjects.Add(result);
742  }
743 
744  } else {
745 
746  // fill the results in the HypoTestInverterResult array
747  fResults->fXValues.push_back(rVal);
748  fResults->fYObjects.Add(result);
749 
750  }
751 
752  // std::cout << "computed value for poi " << rVal << " : " << fResults->GetYValue(fResults->ArraySize()-1)
753  // << " +/- " << fResults->GetYError(fResults->ArraySize()-1) << endl;
754 
755  fScannedVariable->setVal(oldValue);
756 
757  return true;
758 }
759 
760 ////////////////////////////////////////////////////////////////////////////////
761 /// Run an automatic scan until the desired accuracy is reached.
762 /// Start by default from the full interval (min,max) of the POI and then via bisection find the line crossing
763 /// the target line.
764 /// Optionally, a hint can be provided and the scan will be done closer to that value.
765 /// If by bisection the desired accuracy will not be reached, a fit to the points is performed.
766 /// \param[out] limit The limit.
767 /// \param[out] limitErr The error of the limit.
768 /// \param[in] absAccuracy Desired absolute accuracy.
769 /// \param[in] relAccuracy Desired relative accuracy.
770 /// \param[in] hint Hint to start from or nullptr for no hint.
771 
772 bool HypoTestInverter::RunLimit(double &limit, double &limitErr, double absAccuracy, double relAccuracy, const double*hint) const {
773 
774 
775  // routine from G. Petrucciani (from HiggsCombination CMS package)
776 // bool HybridNew::runLimit(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) {
777 
779  //w->loadSnapshot("clean");
780 
781  if ((hint != 0) && (*hint > r->getMin())) {
782  r->setMax(std::min<double>(3.0 * (*hint), r->getMax()));
783  r->setMin(std::max<double>(0.3 * (*hint), r->getMin()));
784  oocoutI((TObject*)0,InputArguments) << "HypoTestInverter::RunLimit - Use hint value " << *hint
785  << " search in interval " << r->getMin() << " , " << r->getMax() << std::endl;
786  }
787 
788  // if not specified use the default values for rel and absolute accuracy
789  if (absAccuracy <= 0) absAccuracy = fgAbsAccuracy;
790  if (relAccuracy <= 0) relAccuracy = fgRelAccuracy;
791 
792  typedef std::pair<double,double> CLs_t;
793  double clsTarget = fSize;
794  CLs_t clsMin(1,0), clsMax(0,0), clsMid(0,0);
795  double rMin = r->getMin(), rMax = r->getMax();
796  limit = 0.5*(rMax + rMin);
797  limitErr = 0.5*(rMax - rMin);
798  bool done = false;
799 
800  TF1 expoFit("expoFit","[0]*exp([1]*(x-[2]))", rMin, rMax);
801 
802  // if (readHybridResults_) {
803  // if (verbose > 0) std::cout << "Search for upper limit using pre-computed grid of p-values" << std::endl;
804 
805  // readAllToysFromFile();
806  // double minDist=1e3;
807  // for (int i = 0, n = limitPlot_->GetN(); i < n; ++i) {
808  // double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i], ey = limitPlot_->GetErrorY(i);
809  // if (verbose > 0) std::cout << " r " << x << (CLs_ ? ", CLs = " : ", CLsplusb = ") << y << " +/- " << ey << std::endl;
810  // if (y-3*ey >= clsTarget) { rMin = x; clsMin = CLs_t(y,ey); }
811  // if (y+3*ey <= clsTarget) { rMax = x; clsMax = CLs_t(y,ey); }
812  // if (fabs(y-clsTarget) < minDist) { limit = x; minDist = fabs(y-clsTarget); }
813  // }
814  // if (verbose > 0) std::cout << " after scan x ~ " << limit << ", bounds [ " << rMin << ", " << rMax << "]" << std::endl;
815  // limitErr = std::max(limit-rMin, rMax-limit);
816  // expoFit.SetRange(rMin,rMax);
817 
818  // if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
819  // if (verbose > 1) std::cout << " reached accuracy " << limitErr << " below " << std::max(rAbsAccuracy_, rRelAccuracy_ * limit) << std::endl;
820  // done = true;
821  // }
822  // } else {
823 
824  fLimitPlot.reset(new TGraphErrors());
825 
826  if (fVerbose > 0) std::cout << "Search for upper limit to the limit" << std::endl;
827  for (int tries = 0; tries < 6; ++tries) {
828  //clsMax = eval(w, mc_s, mc_b, data, rMax);
829  if (! RunOnePoint(rMax) ) {
830  oocoutE((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Hypotest failed" << std::endl;
831  return false;
832  }
833  clsMax = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
834  if (clsMax.first == 0 || clsMax.first + 3 * fabs(clsMax.second) < clsTarget ) break;
835  rMax += rMax;
836  if (tries == 5) {
837  oocoutE((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Cannot set higher limit: at " << r->GetName()
838  << " = " << rMax << " still get "
839  << (fUseCLs ? "CLs" : "CLsplusb") << " = " << clsMax.first << std::endl;
840  return false;
841  }
842  }
843  if (fVerbose > 0) {
844  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Search for lower limit to the limit" << std::endl;
845  }
846  //clsMin = (fUseCLs && rMin == 0 ? CLs_t(1,0) : eval(w, mc_s, mc_b, data, rMin));
847  if ( fUseCLs && rMin == 0 ) {
848  clsMin = CLs_t(1,0);
849  }
850  else {
851  if (! RunOnePoint(rMin) ) return false;
852  clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
853  }
854  if (clsMin.first != 1 && clsMin.first - 3 * fabs(clsMin.second) < clsTarget) {
855  if (fUseCLs) {
856  rMin = 0;
857  clsMin = CLs_t(1,0); // this is always true for CLs
858  } else {
859  rMin = -rMax / 4;
860  for (int tries = 0; tries < 6; ++tries) {
861  if (! RunOnePoint(rMax) ) return false;
862  clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
863  if (clsMin.first == 1 || clsMin.first - 3 * fabs(clsMin.second) > clsTarget) break;
864  rMin += rMin;
865  if (tries == 5) {
866  oocoutE((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Cannot set lower limit: at " << r->GetName()
867  << " = " << rMin << " still get " << (fUseCLs ? "CLs" : "CLsplusb")
868  << " = " << clsMin.first << std::endl;
869  return false;
870  }
871  }
872  }
873  }
874 
875  if (fVerbose > 0)
876  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Now doing proper bracketing & bisection" << std::endl;
877  do {
878 
879  // break loop in case max toys is reached
880  if (fMaxToys > 0 && fTotalToysRun > fMaxToys ) {
881  oocoutW((TObject*)0,Eval) << "HypoTestInverter::RunLimit - maximum number of toys reached " << std::endl;
882  done = false; break;
883  }
884 
885 
886  // determine point by bisection or interpolation
887  limit = 0.5*(rMin+rMax); limitErr = 0.5*(rMax-rMin);
888  if (fgAlgo == "logSecant" && clsMax.first != 0) {
889  double logMin = log(clsMin.first), logMax = log(clsMax.first), logTarget = log(clsTarget);
890  limit = rMin + (rMax-rMin) * (logTarget - logMin)/(logMax - logMin);
891  if (clsMax.second != 0 && clsMin.second != 0) {
892  limitErr = hypot((logTarget-logMax) * (clsMin.second/clsMin.first), (logTarget-logMin) * (clsMax.second/clsMax.first));
893  limitErr *= (rMax-rMin)/((logMax-logMin)*(logMax-logMin));
894  }
895  }
896  r->setError(limitErr);
897 
898  // exit if reached accuracy on r
899  if (limitErr < std::max(absAccuracy, relAccuracy * limit)) {
900  if (fVerbose > 1)
901  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - reached accuracy " << limitErr
902  << " below " << std::max(absAccuracy, relAccuracy * limit) << std::endl;
903  done = true; break;
904  }
905 
906  // evaluate point
907 
908  //clsMid = eval(w, mc_s, mc_b, data, limit, true, clsTarget);
909  if (! RunOnePoint(limit, true, clsTarget) ) return false;
910  clsMid = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
911 
912  if (clsMid.second == -1) {
913  std::cerr << "Hypotest failed" << std::endl;
914  return false;
915  }
916 
917  // if sufficiently far away, drop one of the points
918  if (fabs(clsMid.first-clsTarget) >= 2*clsMid.second) {
919  if ((clsMid.first>clsTarget) == (clsMax.first>clsTarget)) {
920  rMax = limit; clsMax = clsMid;
921  } else {
922  rMin = limit; clsMin = clsMid;
923  }
924  } else {
925  if (fVerbose > 0) std::cout << "Trying to move the interval edges closer" << std::endl;
926  double rMinBound = rMin, rMaxBound = rMax;
927  // try to reduce the size of the interval
928  while (clsMin.second == 0 || fabs(rMin-limit) > std::max(absAccuracy, relAccuracy * limit)) {
929  rMin = 0.5*(rMin+limit);
930  if (!RunOnePoint(rMin,true, clsTarget) ) return false;
931  clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
932  //clsMin = eval(w, mc_s, mc_b, data, rMin, true, clsTarget);
933  if (fabs(clsMin.first-clsTarget) <= 2*clsMin.second) break;
934  rMinBound = rMin;
935  }
936  while (clsMax.second == 0 || fabs(rMax-limit) > std::max(absAccuracy, relAccuracy * limit)) {
937  rMax = 0.5*(rMax+limit);
938 // clsMax = eval(w, mc_s, mc_b, data, rMax, true, clsTarget);
939  if (!RunOnePoint(rMax,true,clsTarget) ) return false;
940  clsMax = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
941  if (fabs(clsMax.first-clsTarget) <= 2*clsMax.second) break;
942  rMaxBound = rMax;
943  }
944  expoFit.SetRange(rMinBound,rMaxBound);
945  break;
946  }
947  } while (true);
948 
949  if (!done) { // didn't reach accuracy with scan, now do fit
950  if (fVerbose) {
951  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Before fit --- \n";
952  std::cout << "Limit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " [" << rMin << ", " << rMax << "]\n";
953  }
954 
955  expoFit.FixParameter(0,clsTarget);
956  expoFit.SetParameter(1,log(clsMax.first/clsMin.first)/(rMax-rMin));
957  expoFit.SetParameter(2,limit);
958  double rMinBound, rMaxBound; expoFit.GetRange(rMinBound, rMaxBound);
959  limitErr = std::max(fabs(rMinBound-limit), fabs(rMaxBound-limit));
960  int npoints = 0;
961 
962  HypoTestInverterPlot plot("plot","plot",fResults);
963  fLimitPlot.reset(plot.MakePlot() );
964 
965 
966  for (int j = 0; j < fLimitPlot->GetN(); ++j) {
967  if (fLimitPlot->GetX()[j] >= rMinBound && fLimitPlot->GetX()[j] <= rMaxBound) npoints++;
968  }
969  for (int i = 0, imax = /*(readHybridResults_ ? 0 : */ 8; i <= imax; ++i, ++npoints) {
970  fLimitPlot->Sort();
971  fLimitPlot->Fit(&expoFit,(fVerbose <= 1 ? "QNR EX0" : "NR EXO"));
972  if (fVerbose) {
973  oocoutI((TObject*)0,Eval) << "Fit to " << npoints << " points: " << expoFit.GetParameter(2) << " +/- " << expoFit.GetParError(2) << std::endl;
974  }
975  if ((rMin < expoFit.GetParameter(2)) && (expoFit.GetParameter(2) < rMax) && (expoFit.GetParError(2) < 0.5*(rMaxBound-rMinBound))) {
976  // sanity check fit result
977  limit = expoFit.GetParameter(2);
978  limitErr = expoFit.GetParError(2);
979  if (limitErr < std::max(absAccuracy, relAccuracy * limit)) break;
980  }
981  // add one point in the interval.
982  double rTry = RooRandom::uniform()*(rMaxBound-rMinBound)+rMinBound;
983  if (i != imax) {
984  if (!RunOnePoint(rTry,true,clsTarget) ) return false;
985  //eval(w, mc_s, mc_b, data, rTry, true, clsTarget);
986  }
987 
988  }
989  }
990 
991 //if (!plot_.empty() && fLimitPlot.get()) {
992  if (fLimitPlot.get() && fLimitPlot->GetN() > 0) {
993  //new TCanvas("c1","c1");
994  fLimitPlot->Sort();
995  fLimitPlot->SetLineWidth(2);
996  double xmin = r->getMin(), xmax = r->getMax();
997  for (int j = 0; j < fLimitPlot->GetN(); ++j) {
998  if (fLimitPlot->GetY()[j] > 1.4*clsTarget || fLimitPlot->GetY()[j] < 0.6*clsTarget) continue;
999  xmin = std::min(fLimitPlot->GetX()[j], xmin);
1000  xmax = std::max(fLimitPlot->GetX()[j], xmax);
1001  }
1002  fLimitPlot->GetXaxis()->SetRangeUser(xmin,xmax);
1003  fLimitPlot->GetYaxis()->SetRangeUser(0.5*clsTarget, 1.5*clsTarget);
1004  fLimitPlot->Draw("AP");
1005  expoFit.Draw("SAME");
1006  TLine line(fLimitPlot->GetX()[0], clsTarget, fLimitPlot->GetX()[fLimitPlot->GetN()-1], clsTarget);
1008  line.DrawLine(limit, 0, limit, fLimitPlot->GetY()[0]);
1010  line.DrawLine(limit-limitErr, 0, limit-limitErr, fLimitPlot->GetY()[0]);
1011  line.DrawLine(limit+limitErr, 0, limit+limitErr, fLimitPlot->GetY()[0]);
1012  //c1->Print(plot_.c_str());
1013  }
1014 
1015  oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Result: \n"
1016  << "\tLimit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " @ " << (1-fSize) * 100 << "% CL\n";
1017  if (fVerbose > 1) oocoutI((TObject*)0,Eval) << "Total toys: " << fTotalToysRun << std::endl;
1018 
1019  // set value in results
1020  fResults->fUpperLimit = limit;
1021  fResults->fUpperLimitError = limitErr;
1022  fResults->fFittedUpperLimit = true;
1023  // lower limit are always min of p value
1026  fResults->fFittedLowerLimit = true;
1027 
1028  return true;
1029 }
1030 
1031 ////////////////////////////////////////////////////////////////////////////////
1032 /// get the distribution of lower limit
1033 /// if rebuild = false (default) it will re-use the results of the scan done
1034 /// for obtained the observed limit and no extra toys will be generated
1035 /// if rebuild a new set of B toys will be done and the procedure will be repeated
1036 /// for each toy
1037 
1039  if (!rebuild) {
1040  if (!fResults) {
1041  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::GetLowerLimitDistribution(false) - result not existing\n";
1042  return 0;
1043  }
1045  }
1046 
1047  TList * clsDist = 0;
1048  TList * clsbDist = 0;
1049  if (fUseCLs) clsDist = &fResults->fExpPValues;
1050  else clsbDist = &fResults->fExpPValues;
1051 
1052  return RebuildDistributions(false, nToys,clsDist, clsbDist);
1053 
1054 }
1055 
1056 ////////////////////////////////////////////////////////////////////////////////
1057 /// get the distribution of lower limit
1058 /// if rebuild = false (default) it will re-use the results of the scan done
1059 /// for obtained the observed limit and no extra toys will be generated
1060 /// if rebuild a new set of B toys will be done and the procedure will be repeated
1061 /// for each toy
1062 /// The nuisance parameter value used for rebuild is the current one in the model
1063 /// so it is user responsibility to set to the desired value (nomi
1064 
1066  if (!rebuild) {
1067  if (!fResults) {
1068  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::GetUpperLimitDistribution(false) - result not existing\n";
1069  return 0;
1070  }
1072  }
1073 
1074  TList * clsDist = 0;
1075  TList * clsbDist = 0;
1076  if (fUseCLs) clsDist = &fResults->fExpPValues;
1077  else clsbDist = &fResults->fExpPValues;
1078 
1079  return RebuildDistributions(true, nToys,clsDist, clsbDist);
1080 }
1081 
1082 ////////////////////////////////////////////////////////////////////////////////
1083 
1085  if (fCalculator0) fCalculator0->SetData(data);
1086 }
1087 
1088 ////////////////////////////////////////////////////////////////////////////////
1089 /// rebuild the sampling distributions by
1090 /// generating some toys and find for each of them a new upper limit
1091 /// Return the upper limit distribution and optionally also the pValue distributions for Cls, Clsb and Clbxs
1092 /// as a TList for each scanned point
1093 /// The method uses the present parameter value. It is user responsibility to give the current parameters to rebuild the distributions
1094 /// It returns a upper or lower limit distribution depending on the isUpper flag, however it computes also the lower limit distribution and it is saved in the
1095 /// output file as an histogram
1096 
1097 SamplingDistribution * HypoTestInverter::RebuildDistributions(bool isUpper, int nToys, TList * clsDist, TList * clsbDist, TList * clbDist, const char *outputfile) {
1098 
1099  if (!fScannedVariable || !fCalculator0) return 0;
1100  // get first background snapshot
1101  const ModelConfig * bModel = fCalculator0->GetAlternateModel();
1102  const ModelConfig * sbModel = fCalculator0->GetNullModel();
1103  if (!bModel || ! sbModel) return 0;
1104  RooArgSet paramPoint;
1105  if (!sbModel->GetParametersOfInterest()) return 0;
1106  paramPoint.add(*sbModel->GetParametersOfInterest());
1107 
1108  const RooArgSet * poibkg = bModel->GetSnapshot();
1109  if (!poibkg) {
1110  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RebuildDistribution - background snapshot not existing"
1111  << " assume is for POI = 0" << std::endl;
1113  paramPoint = RooArgSet(*fScannedVariable);
1114  }
1115  else
1116  paramPoint = *poibkg;
1117  // generate data at bkg parameter point
1118 
1119  ToyMCSampler * toymcSampler = dynamic_cast<ToyMCSampler *>(fCalculator0->GetTestStatSampler() );
1120  if (!toymcSampler) {
1121  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RebuildDistribution - no toy MC sampler existing" << std::endl;
1122  return 0;
1123  }
1124  // set up test stat sampler in case of asymptotic calculator
1125  if (dynamic_cast<RooStats::AsymptoticCalculator*>(fCalculator0) ) {
1126  toymcSampler->SetObservables(*sbModel->GetObservables() );
1127  toymcSampler->SetParametersForTestStat(*sbModel->GetParametersOfInterest());
1128  toymcSampler->SetPdf(*sbModel->GetPdf());
1129  toymcSampler->SetNuisanceParameters(*sbModel->GetNuisanceParameters());
1130  if (sbModel->GetGlobalObservables() ) toymcSampler->SetGlobalObservables(*sbModel->GetGlobalObservables() );
1131  // set number of events
1132  if (!sbModel->GetPdf()->canBeExtended())
1133  toymcSampler->SetNEventsPerToy(1);
1134  }
1135 
1136  // loop on data to generate
1137  int nPoints = fNBins;
1138 
1139  bool storePValues = clsDist || clsbDist || clbDist;
1140  if (fNBins <=0 && storePValues) {
1141  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RebuildDistribution - cannot return p values distribution with the auto scan" << std::endl;
1142  storePValues = false;
1143  nPoints = 0;
1144  }
1145 
1146  if (storePValues) {
1147  if (fResults) nPoints = fResults->ArraySize();
1148  if (nPoints <=0) {
1149  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - result is not existing and number of point to scan is not set"
1150  << std::endl;
1151  return 0;
1152  }
1153  }
1154 
1155  if (nToys <= 0) nToys = 100; // default value
1156 
1157  std::vector<std::vector<double> > CLs_values(nPoints);
1158  std::vector<std::vector<double> > CLsb_values(nPoints);
1159  std::vector<std::vector<double> > CLb_values(nPoints);
1160 
1161  if (storePValues) {
1162  for (int i = 0; i < nPoints; ++i) {
1163  CLs_values[i].reserve(nToys);
1164  CLb_values[i].reserve(nToys);
1165  CLsb_values[i].reserve(nToys);
1166  }
1167  }
1168 
1169  std::vector<double> limit_values; limit_values.reserve(nToys);
1170 
1171  oocoutI((TObject*)0,InputArguments) << "HypoTestInverter - rebuilding the p value distributions by generating ntoys = "
1172  << nToys << std::endl;
1173 
1174 
1175  oocoutI((TObject*)0,InputArguments) << "Rebuilding using parameter of interest point: ";
1177  if (sbModel->GetNuisanceParameters() ) {
1178  oocoutI((TObject*)0,InputArguments) << "And using nuisance parameters: ";
1180  }
1181  // save all parameters to restore them later
1182  assert(bModel->GetPdf() );
1183  assert(bModel->GetObservables() );
1184  RooArgSet * allParams = bModel->GetPdf()->getParameters( *bModel->GetObservables() );
1185  RooArgSet saveParams;
1186  allParams->snapshot(saveParams);
1187 
1188  TFile * fileOut = TFile::Open(outputfile,"RECREATE");
1189  if (!fileOut) {
1190  oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - RebuildDistributions - Error opening file " << outputfile
1191  << " - the resulting limits will not be stored" << std::endl;
1192  }
1193  // create temporary histograms to store the limit result
1194  TH1D * hL = new TH1D("lowerLimitDist","Rebuilt lower limit distribution",100,1.,0.);
1195  TH1D * hU = new TH1D("upperLimitDist","Rebuilt upper limit distribution",100,1.,0.);
1196  TH1D * hN = new TH1D("nObs","Observed events",100,1.,0.);
1197  hL->SetBuffer(2*nToys);
1198  hU->SetBuffer(2*nToys);
1199  std::vector<TH1*> hCLb;
1200  std::vector<TH1*> hCLsb;
1201  std::vector<TH1*> hCLs;
1202  if (storePValues) {
1203  for (int i = 0; i < nPoints; ++i) {
1204  hCLb.push_back(new TH1D(TString::Format("CLbDist_bin%d",i),"CLb distribution",100,1.,0.));
1205  hCLs.push_back(new TH1D(TString::Format("ClsDist_bin%d",i),"CLs distribution",100,1.,0.));
1206  hCLsb.push_back(new TH1D(TString::Format("CLsbDist_bin%d",i),"CLs+b distribution",100,1.,0.));
1207  }
1208  }
1209 
1210 
1211  // loop now on the toys
1212  for (int itoy = 0; itoy < nToys; ++itoy) {
1213 
1214  oocoutP((TObject*)0,Eval) << "\nHypoTestInverter - RebuildDistributions - running toy # " << itoy << " / "
1215  << nToys << std::endl;
1216 
1217 
1218  printf("\n\nshnapshot of s+b model \n");
1219  sbModel->GetSnapshot()->Print("v");
1220 
1221  // reset parameters to initial values to be sure in case they are not reset
1222  if (itoy> 0) *allParams = saveParams;
1223 
1224  // need to set th epdf to clear the cache in ToyMCSampler
1225  // pdf we must use is background pdf
1226  toymcSampler->SetPdf(*bModel->GetPdf() );
1227 
1228 
1229  RooAbsData * bkgdata = toymcSampler->GenerateToyData(paramPoint);
1230 
1231  double nObs = bkgdata->sumEntries();
1232  // for debugging in case of number counting models
1233  if (bkgdata->numEntries() ==1 && !bModel->GetPdf()->canBeExtended()) {
1234  oocoutP((TObject*)0,Generation) << "Generate observables are : ";
1235  RooArgList genObs(*bkgdata->get(0));
1237  nObs = 0;
1238  for (int i = 0; i < genObs.getSize(); ++i) {
1239  RooRealVar * x = dynamic_cast<RooRealVar*>(&genObs[i]);
1240  if (x) nObs += x->getVal();
1241  }
1242  }
1243  hN->Fill(nObs);
1244 
1245  // by copying I will have the same min/max as previous ones
1246  HypoTestInverter inverter = *this;
1247  inverter.SetData(*bkgdata);
1248 
1249  // print global observables
1250  auto gobs = bModel->GetPdf()->getVariables()->selectCommon(* sbModel->GetGlobalObservables() );
1251  gobs->Print("v");
1252 
1253  HypoTestInverterResult * r = inverter.GetInterval();
1254 
1255  if (r == 0) continue;
1256 
1257  double value = (isUpper) ? r->UpperLimit() : r->LowerLimit();
1258  limit_values.push_back( value );
1259  hU->Fill(r->UpperLimit() );
1260  hL->Fill(r->LowerLimit() );
1261 
1262 
1263  std::cout << "The computed upper limit for toy #" << itoy << " is " << value << std::endl;
1264 
1265  // write every 10 toys
1266  if (itoy%10 == 0 || itoy == nToys-1) {
1267  hU->Write("",TObject::kOverwrite);
1268  hL->Write("",TObject::kOverwrite);
1269  hN->Write("",TObject::kOverwrite);
1270  }
1271 
1272  if (!storePValues) continue;
1273 
1274  if (nPoints < r->ArraySize()) {
1275  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter: skip extra points" << std::endl;
1276  }
1277  else if (nPoints > r->ArraySize()) {
1278  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter: missing some points" << std::endl;
1279  }
1280 
1281 
1282  for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
1283  HypoTestResult * hr = r->GetResult(ipoint);
1284  if (hr) {
1285  CLs_values[ipoint].push_back( hr->CLs() );
1286  CLsb_values[ipoint].push_back( hr->CLsplusb() );
1287  CLb_values[ipoint].push_back( hr->CLb() );
1288  hCLs[ipoint]->Fill( hr->CLs() );
1289  hCLb[ipoint]->Fill( hr->CLb() );
1290  hCLsb[ipoint]->Fill( hr->CLsplusb() );
1291  }
1292  else {
1293  oocoutW((TObject*)0,InputArguments) << "HypoTestInverter: missing result for point: x = "
1294  << fResults->GetXValue(ipoint) << std::endl;
1295  }
1296  }
1297  // write every 10 toys
1298  if (itoy%10 == 0 || itoy == nToys-1) {
1299  for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
1300  hCLs[ipoint]->Write("",TObject::kOverwrite);
1301  hCLb[ipoint]->Write("",TObject::kOverwrite);
1302  hCLsb[ipoint]->Write("",TObject::kOverwrite);
1303  }
1304  }
1305 
1306 
1307  delete r;
1308  delete bkgdata;
1309  }
1310 
1311 
1312  if (storePValues) {
1313  if (clsDist) clsDist->SetOwner(true);
1314  if (clbDist) clbDist->SetOwner(true);
1315  if (clsbDist) clsbDist->SetOwner(true);
1316 
1317  oocoutI((TObject*)0,InputArguments) << "HypoTestInverter: storing rebuilt p values " << std::endl;
1318 
1319  for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
1320  if (clsDist) {
1321  TString name = TString::Format("CLs_distrib_%d",ipoint);
1322  clsDist->Add( new SamplingDistribution(name,name,CLs_values[ipoint] ) );
1323  }
1324  if (clbDist) {
1325  TString name = TString::Format("CLb_distrib_%d",ipoint);
1326  clbDist->Add( new SamplingDistribution(name,name,CLb_values[ipoint] ) );
1327  }
1328  if (clsbDist) {
1329  TString name = TString::Format("CLsb_distrib_%d",ipoint);
1330  clsbDist->Add( new SamplingDistribution(name,name,CLsb_values[ipoint] ) );
1331  }
1332  }
1333  }
1334 
1335  if (fileOut) {
1336  fileOut->Close();
1337  }
1338  else {
1339  // delete all the histograms
1340  delete hL;
1341  delete hU;
1342  for (int i = 0; i < nPoints && storePValues; ++i) {
1343  delete hCLs[i];
1344  delete hCLb[i];
1345  delete hCLsb[i];
1346  }
1347  }
1348 
1349  const char * disName = (isUpper) ? "upperLimit_dist" : "lowerLimit_dist";
1350  return new SamplingDistribution(disName, disName, limit_values);
1351 }
virtual Double_t sumEntries() const =0
virtual ~HypoTestInverter()
destructor (delete the HypoTestInverterResult)
void CreateResults() const
create a new HypoTestInverterResult to hold all computed results
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:1910
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:794
bool RunOnePoint(double thisX, bool adaptive=false, double clTarget=-1) const
run only one point at the given POI value
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
void SetBackgroundAsAlt(Bool_t l=kTRUE)
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3274
float xmin
Definition: THbookFile.cxx:93
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
const ModelConfig * GetAlternateModel(void) const
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:249
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3944
A class for performing a hypothesis test inversion by scanning the hypothesis test results of a HypoT...
IntervalCalculator is an interface class for a tools which produce RooStats ConfIntervals.
virtual void SetParametersForTestStat(const RooArgSet &nullpoi)
Definition: ToyMCSampler.h:155
TLine * line
HypoTestResult * Eval(HypoTestCalculatorGeneric &hc, bool adaptive, double clsTarget) const
Run the Hypothesis test at a previous configured point (internal function called by RunOnePoint) ...
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
virtual Double_t CLb() const
Convert NullPValue into a "confidence level".
virtual TestStatistic * GetTestStatistic() const =0
virtual void SetGlobalObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:175
RooAbsCollection * selectCommon(const RooAbsCollection &refColl) const
Create a subset of the current collection, consisting only of those elements that are contained as we...
Definition: Rtypes.h:64
bool RunFixedScan(int nBins, double xMin, double xMax, bool scanLog=false) const
Run a fixed scan.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
#define oocoutF(o, a)
Definition: RooMsgService.h:49
std::unique_ptr< TGraphErrors > fLimitPlot
TestStatSampler is an interface class for a tools which produce RooStats SamplingDistributions.
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
#define oocoutI(o, a)
Definition: RooMsgService.h:45
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:53
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results ...
HypoTestResult is a base class for results from hypothesis tests.
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition: TF1.cxx:3538
virtual TLine * DrawLine(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Draw this line with new coordinates.
Definition: TLine.cxx:89
SamplingDistribution * GetUpperLimitDistribution(bool rebuild=false, int nToys=100)
get the distribution of lower limit if rebuild = false (default) it will re-use the results of the sc...
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
Definition: ToyMCSampler.h:149
HypoTestResult * GetResult(int index) const
return a pointer to the i^th result object
overwrite existing object with same name
Definition: TObject.h:88
Basic string class.
Definition: TString.h:131
static void CheckInputModels(const HypoTestCalculatorGeneric &hc, const RooRealVar &scanVar)
check the model given the given hypotestcalculator
static unsigned int fgNToys
SamplingDistribution * GetAltDistribution(void) const
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1336
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:196
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition: TF1.cxx:1929
Bool_t IsNaN(Double_t x)
Definition: TMath.h:882
Py_ssize_t ArraySize(const std::string &name)
Definition: Utility.cxx:806
STL namespace.
Double_t CLbError() const
The error on the "confidence level" of the null hypothesis.
HypoTestInverterResult * fResults
virtual void SetBuffer(Int_t buffersize, Option_t *option="")
Set the maximum number of entries to be kept in the buffer.
Definition: TH1.cxx:7921
HypoTestInverter & operator=(const HypoTestInverter &rhs)
assignment operator NOTE: this class does not copy the contained result and the HypoTestCalculator, but only the pointers It requires the original HTI to be alive
virtual void SetTestStatistic(TestStatistic *testStatistic)=0
Common base class for the Hypothesis Test Calculators.
#define oocoutP(o, a)
Definition: RooMsgService.h:46
Double_t x[n]
Definition: legend1.C:17
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
Definition: TString.cxx:2311
Double_t CLsError() const
The error on the ratio .
#define oocoutE(o, a)
Definition: RooMsgService.h:48
Bool_t GetPValueIsRightTail(void) const
SamplingDistribution * GetLowerLimitDistribution() const
get expected lower limit distributions implemented using interpolation The size for the sampling dist...
void UseCLs(bool on=true)
flag to switch between using CLsb (default) or CLs as confidence level
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
static void SetCloseProof(Bool_t flag)
set flag to close proof for every new run
virtual void SetData(RooAbsData &data)
Set the DataSet.
HypoTestInverter()
default constructor (doesn&#39;t do anything)
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
A doubly linked list.
Definition: TList.h:44
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Definition: TMath.h:412
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Definition: TMath.h:418
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:225
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
virtual HypoTestResult * GetHypoTest() const
inherited methods from HypoTestCalculator interface
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
TestStatistic * GetTestStatistic() const
return the test statistic which is or will be used by the class
ROOT::R::TRInterface & r
Definition: Object.C:4
TGraphErrors * MakePlot(Option_t *opt="")
return a TGraphErrors with the obtained observed p-values resultinf from the scan By default (Option ...
RooAbsArg * first() const
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:821
Double_t GetTestStatisticData(void) const
int fTotalToysRun
plot of limits
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
void PrintListContent(const RooArgList &l, std::ostream &os=std::cout)
virtual void FixParameter(Int_t ipar, Double_t value)
Fix the value of a parameter The specified value will be used in a fit operation. ...
Definition: TF1.cxx:1561
A simple line.
Definition: TLine.h:22
TList fExpPValues
list of HypoTestResult for each point
Int_t GetSize() const
size of samples
RooRealVar * fScannedVariable
pointer to the generic hypotest calculator used
SamplingDistribution * GetNullDistribution(void) const
HypoTestCalculatorGeneric * fCalculator0
SamplingDistribution * GetUpperLimitDistribution() const
get expected upper limit distributions implemented using interpolation
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
float xmax
Definition: THbookFile.cxx:93
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:71
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
void SetTestStatisticData(const Double_t tsd)
void Close(Option_t *option="") override
Close a file.
Definition: TFile.cxx:874
This class simply holds a sampling distribution of some test statistic.
virtual void Append(const HypoTestResult *other)
add values from another HypoTestResult
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:236
Same purpose as HybridCalculatorOriginal, but different implementation.
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:97
#define h(i)
Definition: RSha256.hxx:106
virtual Double_t CLs() const
is simply (not a method, but a quantity)
virtual Double_t AlternatePValue() const
Return p-value for alternate hypothesis.
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:134
virtual Double_t ConfidenceLevel() const
Get the Confidence level for the test.
virtual Double_t CLsplusb() const
Convert AlternatePValue into a "confidence level".
Namespace for the RooStats classes.
Definition: Asimov.h:19
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
#define ClassImp(name)
Definition: Rtypes.h:361
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:83
RooAbsArg * find(const char *name) const
Find object with given name in list.
HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence inte...
virtual HypoTestInverterResult * GetInterval() const
Run a fixed scan or the automatic scan depending on the configuration.
virtual void SetObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:173
Class to plot an HypoTestInverterResult, result of the HypoTestInverter calculator.
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don&#39;t match any of...
Definition: RooAbsArg.cxx:545
int type
Definition: TGX11.cxx:120
virtual void SetPdf(RooAbsPdf &pdf)
Definition: ToyMCSampler.h:160
Double_t CLsplusbError() const
The error on the "confidence level" of the alternative hypothesis.
virtual RooAbsData * GenerateToyData(RooArgSet &paramPoint, RooAbsPdf &pdf) const
Definition: ToyMCSampler.h:110
const ModelConfig * GetNullModel(void) const
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition: TAttLine.h:42
bool SetTestStatistic(TestStatistic &stat)
set the test statistic to use
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TNamed.cxx:74
#define oocoutW(o, a)
Definition: RooMsgService.h:47
Mother of all ROOT objects.
Definition: TObject.h:37
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return NULL if not existing)
Definition: ModelConfig.h:255
virtual void Add(TObject *obj)
Definition: TList.h:87
double fLowerLimitError
interpolation option (linear or spline)
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:506
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
virtual void GetRange(Double_t *xmin, Double_t *xmax) const
Return range of a generic N-D function.
Definition: TF1.cxx:2282
1-Dim function class
Definition: TF1.h:210
std::unique_ptr< HypoTestCalculatorGeneric > fHC
Does a frequentist hypothesis test.
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
void Clear()
delete contained result and graph
static RooRealVar * GetVariableToScan(const HypoTestCalculatorGeneric &hc)
get the variable to scan try first with null model if not go to alternate model
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:628
bool RunLimit(double &limit, double &limitErr, double absTol=0, double relTol=0, const double *hint=0) const
Run an automatic scan until the desired accuracy is reached.
const RooArgSet * GetSnapshot() const
get RooArgSet for parameters for a particular hypothesis (return NULL if not existing) ...
SamplingDistribution * RebuildDistributions(bool isUpper=true, int nToys=100, TList *clsDist=0, TList *clsbDist=0, TList *clbDist=0, const char *outputfile="HypoTestInverterRebuiltDist.root")
rebuild the sampling distributions by generating some toys and find for each of them a new upper limi...
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:358
double exp(double)
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
virtual void SetNuisanceParameters(const RooArgSet &np)
Definition: ToyMCSampler.h:171
virtual void SetData(RooAbsData &)
Set the DataSet ( add to the the workspace if not already there ?)
static void CloseProof(Option_t *option="s")
close all proof connections
Definition: ProofConfig.h:92
TestStatistic is an interface class to provide a facility for construction test statistics distributi...
Definition: TestStatistic.h:31
int ArraySize() const
number of entries in the results array
SamplingDistribution * GetLowerLimitDistribution(bool rebuild=false, int nToys=100)
get the distribution of lower limit if rebuild = false (default) it will re-use the results of the sc...
Hypothesis Test Calculator based on the asymptotic formulae for the profile likelihood ratio...
virtual void SetConfidenceLevel(Double_t cl)
set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval) ...
char name[80]
Definition: TGX11.cxx:109
double log(double)
TestStatSampler * GetTestStatSampler(void) const
Returns instance of TestStatSampler.
size_t fSize
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
std::vector< double > fXValues
number of points used to build expected p-values