// @(#)root/roostats:$Id$
// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
/*************************************************************************
* Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
* For the list of contributors see $ROOTSYS/README/CREDITS. *
*************************************************************************/
////////////////////////////////////////////////////////////////////////////////
// include other header files
#include "RooAbsData.h"
#
#include "TMath.h"
#include "RooStats/HybridResult.h"
#include "RooStats/HypoTestInverter.h"
#include <cassert>
#include <cmath>
#include "TF1.h"
#include "TFile.h"
#include "TH1.h"
#include "TLine.h"
#include "TCanvas.h"
#include "TGraphErrors.h"
#include "RooRealVar.h"
#include "RooArgSet.h"
#include "RooAbsPdf.h"
#include "RooRandom.h"
#include "RooConstVar.h"
#include "RooMsgService.h"
#include "RooStats/ModelConfig.h"
#include "RooStats/HybridCalculator.h"
#include "RooStats/FrequentistCalculator.h"
#include "RooStats/AsymptoticCalculator.h"
#include "RooStats/SimpleLikelihoodRatioTestStat.h"
#include "RooStats/RatioOfProfiledLikelihoodsTestStat.h"
#include "RooStats/ProfileLikelihoodTestStat.h"
#include "RooStats/ToyMCSampler.h"
#include "RooStats/HypoTestPlot.h"
#include "RooStats/HypoTestInverterPlot.h"
#include "RooStats/ProofConfig.h"
ClassImp(RooStats::HypoTestInverter)
using namespace RooStats;
using namespace std;
// static variable definitions
double HypoTestInverter::fgCLAccuracy = 0.005;
unsigned int HypoTestInverter::fgNToys = 500;
double HypoTestInverter::fgAbsAccuracy = 0.05;
double HypoTestInverter::fgRelAccuracy = 0.05;
std::string HypoTestInverter::fgAlgo = "logSecant";
bool HypoTestInverter::fgCloseProof = false;
// helper class to wrap the functionality of the various HypoTestCalculators
template<class HypoTestType>
struct HypoTestWrapper {
static void SetToys(HypoTestType * h, int toyNull, int toyAlt) { h->SetToys(toyNull,toyAlt); }
};
void HypoTestInverter::SetCloseProof(Bool_t flag) {
// set flag to close proof for every new run
fgCloseProof = flag;
}
RooRealVar * HypoTestInverter::GetVariableToScan(const HypoTestCalculatorGeneric &hc) {
// get the variable to scan
// try first with null model if not go to alternate model
RooRealVar * varToScan = 0;
const ModelConfig * mc = hc.GetNullModel();
if (mc) {
const RooArgSet * poi = mc->GetParametersOfInterest();
if (poi) varToScan = dynamic_cast<RooRealVar*> (poi->first() );
}
if (!varToScan) {
mc = hc.GetAlternateModel();
if (mc) {
const RooArgSet * poi = mc->GetParametersOfInterest();
if (poi) varToScan = dynamic_cast<RooRealVar*> (poi->first() );
}
}
return varToScan;
}
void HypoTestInverter::CheckInputModels(const HypoTestCalculatorGeneric &hc,const RooRealVar & scanVariable) {
// check the model given the given hypotestcalculator
const ModelConfig * modelSB = hc.GetNullModel();
const ModelConfig * modelB = hc.GetAlternateModel();
if (!modelSB || ! modelB)
oocoutF((TObject*)0,InputArguments) << "HypoTestInverter - model are not existing" << std::endl;
assert(modelSB && modelB);
oocoutI((TObject*)0,InputArguments) << "HypoTestInverter ---- Input models: \n"
<< "\t\t using as S+B (null) model : "
<< modelSB->GetName() << "\n"
<< "\t\t using as B (alternate) model : "
<< modelB->GetName() << "\n" << std::endl;
// check if scanVariable is included in B model pdf
RooAbsPdf * bPdf = modelB->GetPdf();
const RooArgSet * bObs = modelB->GetObservables();
if (!bPdf || !bObs) {
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - B model has no pdf or observables defined" << std::endl;
return;
}
RooArgSet * bParams = bPdf->getParameters(*bObs);
if (!bParams) {
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - pdf of B model has no parameters" << std::endl;
return;
}
if (bParams->find(scanVariable.GetName() ) ) {
const RooArgSet * poiB = modelB->GetSnapshot();
if (!poiB || !poiB->find(scanVariable.GetName()) ||
( (RooRealVar*) poiB->find(scanVariable.GetName()) )->getVal() != 0 )
oocoutW((TObject*)0,InputArguments) << "HypoTestInverter - using a B model with POI "
<< scanVariable.GetName() << " not equal to zero "
<< " user must check input model configurations " << endl;
if (poiB) delete poiB;
}
delete bParams;
}
HypoTestInverter::HypoTestInverter( ) :
fTotalToysRun(0),
fMaxToys(0),
fCalculator0(0),
fScannedVariable(0),
fResults(0),
fUseCLs(false),
fScanLog(false),
fSize(0),
fVerbose(0),
fCalcType(kUndefined),
fNBins(0), fXmin(1), fXmax(1),
fNumErr(0)
{
// default constructor (doesn't do anything)
}
HypoTestInverter::HypoTestInverter( HypoTestCalculatorGeneric& hc,
RooRealVar* scannedVariable, double size ) :
fTotalToysRun(0),
fMaxToys(0),
fCalculator0(0),
fScannedVariable(scannedVariable),
fResults(0),
fUseCLs(false),
fScanLog(false),
fSize(size),
fVerbose(0),
fCalcType(kUndefined),
fNBins(0), fXmin(1), fXmax(1),
fNumErr(0)
{
// Constructor from a HypoTestCalculatorGeneric
// The HypoTest calculator must be a FrequentistCalculator or HybridCalculator type
// Other type of calculators are not supported.
// The calculator must be created before by using the S+B model for the null and
// the B model for the alt
// If no variable to scan are given they are assumed to be the first variable
// from the parameter of interests of the null model
if (!fScannedVariable) {
fScannedVariable = HypoTestInverter::GetVariableToScan(hc);
}
if (!fScannedVariable)
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
else
CheckInputModels(hc,*fScannedVariable);
HybridCalculator * hybCalc = dynamic_cast<HybridCalculator*>(&hc);
if (hybCalc) {
fCalcType = kHybrid;
fCalculator0 = hybCalc;
return;
}
FrequentistCalculator * freqCalc = dynamic_cast<FrequentistCalculator*>(&hc);
if (freqCalc) {
fCalcType = kFrequentist;
fCalculator0 = freqCalc;
return;
}
AsymptoticCalculator * asymCalc = dynamic_cast<AsymptoticCalculator*>(&hc);
if (asymCalc) {
fCalcType = kAsymptotic;
fCalculator0 = asymCalc;
return;
}
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Type of hypotest calculator is not supported " <<std::endl;
fCalculator0 = &hc;
}
HypoTestInverter::HypoTestInverter( HybridCalculator& hc,
RooRealVar* scannedVariable, double size ) :
fTotalToysRun(0),
fMaxToys(0),
fCalculator0(&hc),
fScannedVariable(scannedVariable),
fResults(0),
fUseCLs(false),
fScanLog(false),
fSize(size),
fVerbose(0),
fCalcType(kHybrid),
fNBins(0), fXmin(1), fXmax(1),
fNumErr(0)
{
// Constructor from a reference to a HybridCalculator
// The calculator must be created before by using the S+B model for the null and
// the B model for the alt
// If no variable to scan are given they are assumed to be the first variable
// from the parameter of interests of the null model
if (!fScannedVariable) {
fScannedVariable = GetVariableToScan(hc);
}
if (!fScannedVariable)
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
else
CheckInputModels(hc,*fScannedVariable);
}
HypoTestInverter::HypoTestInverter( FrequentistCalculator& hc,
RooRealVar* scannedVariable, double size ) :
fTotalToysRun(0),
fMaxToys(0),
fCalculator0(&hc),
fScannedVariable(scannedVariable),
fResults(0),
fUseCLs(false),
fScanLog(false),
fSize(size),
fVerbose(0),
fCalcType(kFrequentist),
fNBins(0), fXmin(1), fXmax(1),
fNumErr(0)
{
// Constructor from a reference to a FrequentistCalculator
// The calculator must be created before by using the S+B model for the null and
// the B model for the alt
// If no variable to scan are given they are assumed to be the first variable
// from the parameter of interests of the null model
if (!fScannedVariable) {
fScannedVariable = GetVariableToScan(hc);
}
if (!fScannedVariable)
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
else
CheckInputModels(hc,*fScannedVariable);
}
HypoTestInverter::HypoTestInverter( AsymptoticCalculator& hc,
RooRealVar* scannedVariable, double size ) :
fTotalToysRun(0),
fMaxToys(0),
fCalculator0(&hc),
fScannedVariable(scannedVariable),
fResults(0),
fUseCLs(false),
fScanLog(false),
fSize(size),
fVerbose(0),
fCalcType(kAsymptotic),
fNBins(0), fXmin(1), fXmax(1),
fNumErr(0)
{
// Constructor from a reference to a AsymptoticCalculator
// The calculator must be created before by using the S+B model for the null and
// the B model for the alt
// If no variable to scan are given they are assumed to be the first variable
// from the parameter of interests of the null model
if (!fScannedVariable) {
fScannedVariable = GetVariableToScan(hc);
}
if (!fScannedVariable)
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
else
CheckInputModels(hc,*fScannedVariable);
}
////////////////////////////////////////////////////////////////////////////////
/// Constructor from a model for B model and a model for S+B.
/// An HypoTestCalculator (Hybrid of Frequentis) will be created using the
/// S+B model as the null and the B model as the alternate
/// If no variable to scan are given they are assumed to be the first variable
/// from the parameter of interests of the null model
HypoTestInverter::HypoTestInverter( RooAbsData& data, ModelConfig &sbModel, ModelConfig &bModel,
RooRealVar * scannedVariable, ECalculatorType type, double size) :
fTotalToysRun(0),
fMaxToys(0),
fCalculator0(0),
fScannedVariable(scannedVariable),
fResults(0),
fUseCLs(false),
fScanLog(false),
fSize(size),
fVerbose(0),
fCalcType(type),
fNBins(0), fXmin(1), fXmax(1),
fNumErr(0)
{
if(fCalcType==kFrequentist) fHC.reset(new FrequentistCalculator(data, bModel, sbModel));
if(fCalcType==kHybrid) fHC.reset( new HybridCalculator(data, bModel, sbModel)) ;
if(fCalcType==kAsymptotic) fHC.reset( new AsymptoticCalculator(data, bModel, sbModel));
fCalculator0 = fHC.get();
// get scanned variabke
if (!fScannedVariable) {
fScannedVariable = GetVariableToScan(*fCalculator0);
}
if (!fScannedVariable)
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - Cannot guess the variable to scan " << std::endl;
else
CheckInputModels(*fCalculator0,*fScannedVariable);
}
HypoTestInverter::HypoTestInverter(const HypoTestInverter & rhs) :
IntervalCalculator(),
fTotalToysRun(0),
fCalculator0(0), fScannedVariable(0), // add these for Coverity
fResults(0)
{
// copy-constructor
// NOTE: this class does not copy the contained result and
// the HypoTestCalculator, but only the pointers
// It requires the original HTI to be alive
(*this) = rhs;
}
HypoTestInverter & 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
if (this == &rhs) return *this;
fTotalToysRun = 0;
fMaxToys = rhs.fMaxToys;
fCalculator0 = rhs.fCalculator0;
fScannedVariable = rhs.fScannedVariable;
fUseCLs = rhs.fUseCLs;
fScanLog = rhs.fScanLog;
fSize = rhs.fSize;
fVerbose = rhs.fVerbose;
fCalcType = rhs.fCalcType;
fNBins = rhs.fNBins;
fXmin = rhs.fXmin;
fXmax = rhs.fXmax;
fNumErr = rhs.fNumErr;
return *this;
}
HypoTestInverter::~HypoTestInverter()
{
// destructor (delete the HypoTestInverterResult)
if (fResults) delete fResults;
fCalculator0 = 0;
}
TestStatistic * HypoTestInverter::GetTestStatistic( ) const
{
// return the test statistic which is or will be used by the class
if(fCalculator0 && fCalculator0->GetTestStatSampler()){
return fCalculator0->GetTestStatSampler()->GetTestStatistic();
}
else
return 0;
}
bool HypoTestInverter::SetTestStatistic(TestStatistic& stat)
{
// set the test statistic to use
if(fCalculator0 && fCalculator0->GetTestStatSampler()){
fCalculator0->GetTestStatSampler()->SetTestStatistic(&stat);
return true;
}
else return false;
}
void HypoTestInverter::Clear() {
// delete contained result and graph
if (fResults) delete fResults;
fResults = 0;
fLimitPlot.reset(nullptr);
}
void HypoTestInverter::CreateResults() const {
// create a new HypoTestInverterResult to hold all computed results
if (fResults == 0) {
TString results_name = "result_";
results_name += fScannedVariable->GetName();
fResults = new HypoTestInverterResult(results_name,*fScannedVariable,ConfidenceLevel());
TString title = "HypoTestInverter Result For ";
title += fScannedVariable->GetName();
fResults->SetTitle(title);
}
fResults->UseCLs(fUseCLs);
fResults->SetConfidenceLevel(1.-fSize);
// check if one or two sided scan
if (fCalculator0) {
// if asymptotic calculator
AsymptoticCalculator * ac = dynamic_cast<AsymptoticCalculator*>(fCalculator0);
if (ac)
fResults->fIsTwoSided = ac->IsTwoSided();
else {
// in case of the other calculators
TestStatSampler * sampler = fCalculator0->GetTestStatSampler();
if (sampler) {
ProfileLikelihoodTestStat * pl = dynamic_cast<ProfileLikelihoodTestStat*>(sampler->GetTestStatistic());
if (pl && pl->IsTwoSided() ) fResults->fIsTwoSided = true;
}
}
}
}
HypoTestInverterResult* HypoTestInverter::GetInterval() const {
// Run a fixed scan or the automatic scan depending on the configuration
// Return if needed a copy of the result object which will be managed by the user
// if having a result with at least one point return it
if (fResults && fResults->ArraySize() >= 1) {
oocoutI((TObject*)0,Eval) << "HypoTestInverter::GetInterval - return an already existing interval " << std::endl;
return (HypoTestInverterResult*)(fResults->Clone());
}
if (fNBins > 0) {
oocoutI((TObject*)0,Eval) << "HypoTestInverter::GetInterval - run a fixed scan" << std::endl;
bool ret = RunFixedScan(fNBins, fXmin, fXmax, fScanLog);
if (!ret)
oocoutE((TObject*)0,Eval) << "HypoTestInverter::GetInterval - error running a fixed scan " << std::endl;
}
else {
oocoutI((TObject*)0,Eval) << "HypoTestInverter::GetInterval - run an automatic scan" << std::endl;
double limit(0),err(0);
bool ret = RunLimit(limit,err);
if (!ret)
oocoutE((TObject*)0,Eval) << "HypoTestInverter::GetInterval - error running an auto scan " << std::endl;
}
if (fgCloseProof) ProofConfig::CloseProof();
return (HypoTestInverterResult*) (fResults->Clone());
}
HypoTestResult * HypoTestInverter::Eval(HypoTestCalculatorGeneric &hc, bool adaptive, double clsTarget) const {
// Run the Hypothesis test at a previous configured point
//(internal function called by RunOnePoint)
//for debug
// std::cout << ">>>>>>>>>>> " << std::endl;
// std::cout << "alternate model " << std::endl;
// hc.GetAlternateModel()->GetNuisanceParameters()->Print("V");
// hc.GetAlternateModel()->GetParametersOfInterest()->Print("V");
// std::cout << "Null model " << std::endl;
// hc.GetNullModel()->GetNuisanceParameters()->Print("V");
// hc.GetNullModel()->GetParametersOfInterest()->Print("V");
// std::cout << "<<<<<<<<<<<<<<< " << std::endl;
// run the hypothesis test
HypoTestResult * hcResult = hc.GetHypoTest();
if (hcResult == 0) {
oocoutE((TObject*)0,Eval) << "HypoTestInverter::Eval - HypoTest failed" << std::endl;
return hcResult;
}
// since the b model is the alt need to set the flag
hcResult->SetBackgroundAsAlt(true);
// bool flipPvalue = false;
// if (flipPValues)
// hcResult->SetPValueIsRightTail(!hcResult->GetPValueIsRightTail());
// adjust for some numerical error in discrete models and == is not anymore
if (hcResult->GetPValueIsRightTail() )
hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()-fNumErr); // issue with < vs <= in discrete models
else
hcResult->SetTestStatisticData(hcResult->GetTestStatisticData()+fNumErr); // issue with < vs <= in discrete models
double clsMid = (fUseCLs ? hcResult->CLs() : hcResult->CLsplusb());
double clsMidErr = (fUseCLs ? hcResult->CLsError() : hcResult->CLsplusbError());
//if (fVerbose) std::cout << (fUseCLs ? "\tCLs = " : "\tCLsplusb = ") << clsMid << " +/- " << clsMidErr << std::endl;
if (adaptive) {
if (fCalcType == kHybrid) HypoTestWrapper<HybridCalculator>::SetToys((HybridCalculator*)&hc, fUseCLs ? fgNToys : 1, 4*fgNToys);
if (fCalcType == kFrequentist) HypoTestWrapper<FrequentistCalculator>::SetToys((FrequentistCalculator*)&hc, fUseCLs ? fgNToys : 1, 4*fgNToys);
while (clsMidErr >= fgCLAccuracy && (clsTarget == -1 || fabs(clsMid-clsTarget) < 3*clsMidErr) ) {
std::unique_ptr<HypoTestResult> more(hc.GetHypoTest());
// if (flipPValues)
// more->SetPValueIsRightTail(!more->GetPValueIsRightTail());
hcResult->Append(more.get());
clsMid = (fUseCLs ? hcResult->CLs() : hcResult->CLsplusb());
clsMidErr = (fUseCLs ? hcResult->CLsError() : hcResult->CLsplusbError());
if (fVerbose) std::cout << (fUseCLs ? "\tCLs = " : "\tCLsplusb = ") << clsMid << " +/- " << clsMidErr << std::endl;
}
}
if (fVerbose ) {
oocoutP((TObject*)0,Eval) << "P values for " << fScannedVariable->GetName() << " = " <<
fScannedVariable->getVal() << "\n" <<
"\tCLs = " << hcResult->CLs() << " +/- " << hcResult->CLsError() << "\n" <<
"\tCLb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << "\n" <<
"\tCLsplusb = " << hcResult->CLsplusb() << " +/- " << hcResult->CLsplusbError() << "\n" <<
std::endl;
}
if (fCalcType == kFrequentist || fCalcType == kHybrid) {
fTotalToysRun += (hcResult->GetAltDistribution()->GetSize() + hcResult->GetNullDistribution()->GetSize());
// set sampling distribution name
TString nullDistName = TString::Format("%s_%s_%4.2f",hcResult->GetNullDistribution()->GetName(),
fScannedVariable->GetName(), fScannedVariable->getVal() );
TString altDistName = TString::Format("%s_%s_%4.2f",hcResult->GetAltDistribution()->GetName(),
fScannedVariable->GetName(), fScannedVariable->getVal() );
hcResult->GetNullDistribution()->SetName(nullDistName);
hcResult->GetAltDistribution()->SetName(altDistName);
}
return hcResult;
}
bool HypoTestInverter::RunFixedScan( int nBins, double xMin, double xMax, bool scanLog ) const
{
// Run a Fixed scan in npoints between min and max
CreateResults();
// interpolate the limits
fResults->fFittedLowerLimit = false;
fResults->fFittedUpperLimit = false;
// safety checks
if ( nBins<=0 ) {
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - Please provide nBins>0\n";
return false;
}
if ( nBins==1 && xMin!=xMax ) {
oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - nBins==1 -> I will run for xMin (" << xMin << ")\n";
}
if ( xMin==xMax && nBins>1 ) {
oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - xMin==xMax -> I will enforce nBins==1\n";
nBins = 1;
}
if ( xMin>xMax ) {
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - Please provide xMin ("
<< xMin << ") smaller that xMax (" << xMax << ")\n";
return false;
}
if (xMin < fScannedVariable->getMin()) {
xMin = fScannedVariable->getMin();
oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - xMin < lower bound, use xmin = "
<< xMin << std::endl;
}
if (xMax > fScannedVariable->getMax()) {
xMax = fScannedVariable->getMax();
oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RunFixedScan - xMax > upper bound, use xmax = "
<< xMax << std::endl;
}
double thisX = xMin;
for (int i=0; i<nBins; i++) {
if (i > 0) { // avoids case of nBins = 1
if (scanLog)
thisX = exp( log(xMin) + i*(log(xMax)-log(xMin))/(nBins-1) ); // scan in log x
else
thisX = xMin + i*(xMax-xMin)/(nBins-1); // linear scan in x
}
bool status = RunOnePoint(thisX);
// check if failed status
if ( status==false ) {
std::cout << "\t\tLoop interrupted because of failed status\n";
return false;
}
}
return true;
}
bool HypoTestInverter::RunOnePoint( double rVal, bool adaptive, double clTarget) const
{
// run only one point at the given POI value
CreateResults();
// check if rVal is in the range specified for fScannedVariable
if ( rVal < fScannedVariable->getMin() ) {
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunOnePoint - Out of range: using the lower bound "
<< fScannedVariable->getMin()
<< " on the scanned variable rather than " << rVal<< "\n";
rVal = fScannedVariable->getMin();
}
if ( rVal > fScannedVariable->getMax() ) {
// print a message when you have a significative difference since rval is computed
if ( rVal > fScannedVariable->getMax()*(1.+1.E-12) )
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RunOnePoint - Out of range: using the upper bound "
<< fScannedVariable->getMax()
<< " on the scanned variable rather than " << rVal<< "\n";
rVal = fScannedVariable->getMax();
}
// save old value
double oldValue = fScannedVariable->getVal();
// evaluate hybrid calculator at a single point
fScannedVariable->setVal(rVal);
// need to set value of rval in hybridcalculator
// assume null model is S+B and alternate is B only
const ModelConfig * sbModel = fCalculator0->GetNullModel();
RooArgSet poi; poi.add(*sbModel->GetParametersOfInterest());
// set poi to right values
poi = RooArgSet(*fScannedVariable);
const_cast<ModelConfig*>(sbModel)->SetSnapshot(poi);
if (fVerbose > 0)
oocoutP((TObject*)0,Eval) << "Running for " << fScannedVariable->GetName() << " = " << fScannedVariable->getVal() << endl;
// compute the results
HypoTestResult* result = Eval(*fCalculator0,adaptive,clTarget);
if (!result) {
oocoutE((TObject*)0,Eval) << "HypoTestInverter - Error running point " << fScannedVariable->GetName() << " = " <<
fScannedVariable->getVal() << endl;
return false;
}
// in case of a dummy result
if (TMath::IsNaN(result->NullPValue() ) && TMath::IsNaN(result->AlternatePValue() ) ) {
oocoutW((TObject*)0,Eval) << "HypoTestInverter - Skip invalid result for point " << fScannedVariable->GetName() << " = " <<
fScannedVariable->getVal() << endl;
return true; // need to return true to avoid breaking the scan loop
}
double lastXtested;
if ( fResults->ArraySize()!=0 ) lastXtested = fResults->GetXValue(fResults->ArraySize()-1);
else lastXtested = -999;
if ( (std::abs(rVal) < 1 && TMath::AreEqualAbs(rVal, lastXtested,1.E-12) ) ||
(std::abs(rVal) >= 1 && TMath::AreEqualRel(rVal, lastXtested,1.E-12) ) ) {
oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunOnePoint - Merge with previous result for "
<< fScannedVariable->GetName() << " = " << rVal << std::endl;
HypoTestResult* prevResult = fResults->GetResult(fResults->ArraySize()-1);
if (prevResult && prevResult->GetNullDistribution() && prevResult->GetAltDistribution()) {
prevResult->Append(result);
delete result; // we can delete the result
}
else {
// if it was empty we re-use it
oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunOnePoint - replace previous empty result\n";
fResults->fYObjects.Remove( prevResult);
fResults->fYObjects.Add(result);
}
} else {
// fill the results in the HypoTestInverterResult array
fResults->fXValues.push_back(rVal);
fResults->fYObjects.Add(result);
}
// std::cout << "computed value for poi " << rVal << " : " << fResults->GetYValue(fResults->ArraySize()-1)
// << " +/- " << fResults->GetYError(fResults->ArraySize()-1) << endl;
fScannedVariable->setVal(oldValue);
return true;
}
bool HypoTestInverter::RunLimit(double &limit, double &limitErr, double absAccuracy, double relAccuracy, const double*hint) const {
// run an automatic scan until the desired accurancy is reached
// Start by default from the full interval (min,max) of the POI and then via bisection find the line crossing
// the target line
// Optionally an hint can be provided and the scan will be done closer to that value
// If by bisection the desired accuracy will not be reached a fit to the points is performed
// routine from G. Petrucciani (from HiggsCombination CMS package)
// bool HybridNew::runLimit(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) {
RooRealVar *r = fScannedVariable;
//w->loadSnapshot("clean");
if ((hint != 0) && (*hint > r->getMin())) {
r->setMax(std::min<double>(3.0 * (*hint), r->getMax()));
r->setMin(std::max<double>(0.3 * (*hint), r->getMin()));
oocoutI((TObject*)0,InputArguments) << "HypoTestInverter::RunLimit - Use hint value " << *hint
<< " search in interval " << r->getMin() << " , " << r->getMax() << std::endl;
}
// if not specified use the default values for rel and absolute accuracy
if (absAccuracy <= 0) absAccuracy = fgAbsAccuracy;
if (relAccuracy <= 0) relAccuracy = fgRelAccuracy;
typedef std::pair<double,double> CLs_t;
double clsTarget = fSize;
CLs_t clsMin(1,0), clsMax(0,0), clsMid(0,0);
double rMin = r->getMin(), rMax = r->getMax();
limit = 0.5*(rMax + rMin);
limitErr = 0.5*(rMax - rMin);
bool done = false;
TF1 expoFit("expoFit","[0]*exp([1]*(x-[2]))", rMin, rMax);
// if (readHybridResults_) {
// if (verbose > 0) std::cout << "Search for upper limit using pre-computed grid of p-values" << std::endl;
// readAllToysFromFile();
// double minDist=1e3;
// for (int i = 0, n = limitPlot_->GetN(); i < n; ++i) {
// double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i], ey = limitPlot_->GetErrorY(i);
// if (verbose > 0) std::cout << " r " << x << (CLs_ ? ", CLs = " : ", CLsplusb = ") << y << " +/- " << ey << std::endl;
// if (y-3*ey >= clsTarget) { rMin = x; clsMin = CLs_t(y,ey); }
// if (y+3*ey <= clsTarget) { rMax = x; clsMax = CLs_t(y,ey); }
// if (fabs(y-clsTarget) < minDist) { limit = x; minDist = fabs(y-clsTarget); }
// }
// if (verbose > 0) std::cout << " after scan x ~ " << limit << ", bounds [ " << rMin << ", " << rMax << "]" << std::endl;
// limitErr = std::max(limit-rMin, rMax-limit);
// expoFit.SetRange(rMin,rMax);
// if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
// if (verbose > 1) std::cout << " reached accuracy " << limitErr << " below " << std::max(rAbsAccuracy_, rRelAccuracy_ * limit) << std::endl;
// done = true;
// }
// } else {
fLimitPlot.reset(new TGraphErrors());
if (fVerbose > 0) std::cout << "Search for upper limit to the limit" << std::endl;
for (int tries = 0; tries < 6; ++tries) {
//clsMax = eval(w, mc_s, mc_b, data, rMax);
if (! RunOnePoint(rMax) ) {
oocoutE((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Hypotest failed" << std::endl;
return false;
}
clsMax = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
if (clsMax.first == 0 || clsMax.first + 3 * fabs(clsMax.second) < clsTarget ) break;
rMax += rMax;
if (tries == 5) {
oocoutE((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Cannot set higher limit: at " << r->GetName()
<< " = " << rMax << " still get "
<< (fUseCLs ? "CLs" : "CLsplusb") << " = " << clsMax.first << std::endl;
return false;
}
}
if (fVerbose > 0) {
oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Search for lower limit to the limit" << std::endl;
}
//clsMin = (fUseCLs && rMin == 0 ? CLs_t(1,0) : eval(w, mc_s, mc_b, data, rMin));
if ( fUseCLs && rMin == 0 ) {
clsMin = CLs_t(1,0);
}
else {
if (! RunOnePoint(rMin) ) return false;
clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
}
if (clsMin.first != 1 && clsMin.first - 3 * fabs(clsMin.second) < clsTarget) {
if (fUseCLs) {
rMin = 0;
clsMin = CLs_t(1,0); // this is always true for CLs
} else {
rMin = -rMax / 4;
for (int tries = 0; tries < 6; ++tries) {
if (! RunOnePoint(rMax) ) return false;
clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
if (clsMin.first == 1 || clsMin.first - 3 * fabs(clsMin.second) > clsTarget) break;
rMin += rMin;
if (tries == 5) {
oocoutE((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Cannot set lower limit: at " << r->GetName()
<< " = " << rMin << " still get " << (fUseCLs ? "CLs" : "CLsplusb")
<< " = " << clsMin.first << std::endl;
return false;
}
}
}
}
if (fVerbose > 0)
oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Now doing proper bracketing & bisection" << std::endl;
do {
// break loop in case max toys is reached
if (fMaxToys > 0 && fTotalToysRun > fMaxToys ) {
oocoutW((TObject*)0,Eval) << "HypoTestInverter::RunLimit - maximum number of toys reached " << std::endl;
done = false; break;
}
// determine point by bisection or interpolation
limit = 0.5*(rMin+rMax); limitErr = 0.5*(rMax-rMin);
if (fgAlgo == "logSecant" && clsMax.first != 0) {
double logMin = log(clsMin.first), logMax = log(clsMax.first), logTarget = log(clsTarget);
limit = rMin + (rMax-rMin) * (logTarget - logMin)/(logMax - logMin);
if (clsMax.second != 0 && clsMin.second != 0) {
limitErr = hypot((logTarget-logMax) * (clsMin.second/clsMin.first), (logTarget-logMin) * (clsMax.second/clsMax.first));
limitErr *= (rMax-rMin)/((logMax-logMin)*(logMax-logMin));
}
}
r->setError(limitErr);
// exit if reached accuracy on r
if (limitErr < std::max(absAccuracy, relAccuracy * limit)) {
if (fVerbose > 1)
oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - reached accuracy " << limitErr
<< " below " << std::max(absAccuracy, relAccuracy * limit) << std::endl;
done = true; break;
}
// evaluate point
//clsMid = eval(w, mc_s, mc_b, data, limit, true, clsTarget);
if (! RunOnePoint(limit, true, clsTarget) ) return false;
clsMid = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
if (clsMid.second == -1) {
std::cerr << "Hypotest failed" << std::endl;
return false;
}
// if sufficiently far away, drop one of the points
if (fabs(clsMid.first-clsTarget) >= 2*clsMid.second) {
if ((clsMid.first>clsTarget) == (clsMax.first>clsTarget)) {
rMax = limit; clsMax = clsMid;
} else {
rMin = limit; clsMin = clsMid;
}
} else {
if (fVerbose > 0) std::cout << "Trying to move the interval edges closer" << std::endl;
double rMinBound = rMin, rMaxBound = rMax;
// try to reduce the size of the interval
while (clsMin.second == 0 || fabs(rMin-limit) > std::max(absAccuracy, relAccuracy * limit)) {
rMin = 0.5*(rMin+limit);
if (!RunOnePoint(rMin,true, clsTarget) ) return false;
clsMin = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
//clsMin = eval(w, mc_s, mc_b, data, rMin, true, clsTarget);
if (fabs(clsMin.first-clsTarget) <= 2*clsMin.second) break;
rMinBound = rMin;
}
while (clsMax.second == 0 || fabs(rMax-limit) > std::max(absAccuracy, relAccuracy * limit)) {
rMax = 0.5*(rMax+limit);
// clsMax = eval(w, mc_s, mc_b, data, rMax, true, clsTarget);
if (!RunOnePoint(rMax,true,clsTarget) ) return false;
clsMax = std::make_pair( fResults->GetLastYValue(), fResults->GetLastYError() );
if (fabs(clsMax.first-clsTarget) <= 2*clsMax.second) break;
rMaxBound = rMax;
}
expoFit.SetRange(rMinBound,rMaxBound);
break;
}
} while (true);
if (!done) { // didn't reach accuracy with scan, now do fit
if (fVerbose) {
oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Before fit --- \n";
std::cout << "Limit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " [" << rMin << ", " << rMax << "]\n";
}
expoFit.FixParameter(0,clsTarget);
expoFit.SetParameter(1,log(clsMax.first/clsMin.first)/(rMax-rMin));
expoFit.SetParameter(2,limit);
double rMinBound, rMaxBound; expoFit.GetRange(rMinBound, rMaxBound);
limitErr = std::max(fabs(rMinBound-limit), fabs(rMaxBound-limit));
int npoints = 0;
HypoTestInverterPlot plot("plot","plot",fResults);
fLimitPlot.reset(plot.MakePlot() );
for (int j = 0; j < fLimitPlot->GetN(); ++j) {
if (fLimitPlot->GetX()[j] >= rMinBound && fLimitPlot->GetX()[j] <= rMaxBound) npoints++;
}
for (int i = 0, imax = /*(readHybridResults_ ? 0 : */ 8; i <= imax; ++i, ++npoints) {
fLimitPlot->Sort();
fLimitPlot->Fit(&expoFit,(fVerbose <= 1 ? "QNR EX0" : "NR EXO"));
if (fVerbose) {
oocoutI((TObject*)0,Eval) << "Fit to " << npoints << " points: " << expoFit.GetParameter(2) << " +/- " << expoFit.GetParError(2) << std::endl;
}
if ((rMin < expoFit.GetParameter(2)) && (expoFit.GetParameter(2) < rMax) && (expoFit.GetParError(2) < 0.5*(rMaxBound-rMinBound))) {
// sanity check fit result
limit = expoFit.GetParameter(2);
limitErr = expoFit.GetParError(2);
if (limitErr < std::max(absAccuracy, relAccuracy * limit)) break;
}
// add one point in the interval.
double rTry = RooRandom::uniform()*(rMaxBound-rMinBound)+rMinBound;
if (i != imax) {
if (!RunOnePoint(rTry,true,clsTarget) ) return false;
//eval(w, mc_s, mc_b, data, rTry, true, clsTarget);
}
}
}
//if (!plot_.empty() && fLimitPlot.get()) {
if (fLimitPlot.get() && fLimitPlot->GetN() > 0) {
//new TCanvas("c1","c1");
fLimitPlot->Sort();
fLimitPlot->SetLineWidth(2);
double xmin = r->getMin(), xmax = r->getMax();
for (int j = 0; j < fLimitPlot->GetN(); ++j) {
if (fLimitPlot->GetY()[j] > 1.4*clsTarget || fLimitPlot->GetY()[j] < 0.6*clsTarget) continue;
xmin = std::min(fLimitPlot->GetX()[j], xmin);
xmax = std::max(fLimitPlot->GetX()[j], xmax);
}
fLimitPlot->GetXaxis()->SetRangeUser(xmin,xmax);
fLimitPlot->GetYaxis()->SetRangeUser(0.5*clsTarget, 1.5*clsTarget);
fLimitPlot->Draw("AP");
expoFit.Draw("SAME");
TLine line(fLimitPlot->GetX()[0], clsTarget, fLimitPlot->GetX()[fLimitPlot->GetN()-1], clsTarget);
line.SetLineColor(kRed); line.SetLineWidth(2); line.Draw();
line.DrawLine(limit, 0, limit, fLimitPlot->GetY()[0]);
line.SetLineWidth(1); line.SetLineStyle(2);
line.DrawLine(limit-limitErr, 0, limit-limitErr, fLimitPlot->GetY()[0]);
line.DrawLine(limit+limitErr, 0, limit+limitErr, fLimitPlot->GetY()[0]);
//c1->Print(plot_.c_str());
}
oocoutI((TObject*)0,Eval) << "HypoTestInverter::RunLimit - Result: \n"
<< "\tLimit: " << r->GetName() << " < " << limit << " +/- " << limitErr << " @ " << (1-fSize) * 100 << "% CL\n";
if (fVerbose > 1) oocoutI((TObject*)0,Eval) << "Total toys: " << fTotalToysRun << std::endl;
// set value in results
fResults->fUpperLimit = limit;
fResults->fUpperLimitError = limitErr;
fResults->fFittedUpperLimit = true;
// lower limit are always min of p value
fResults->fLowerLimit = fScannedVariable->getMin();
fResults->fLowerLimitError = 0;
fResults->fFittedLowerLimit = true;
return true;
}
SamplingDistribution * HypoTestInverter::GetLowerLimitDistribution(bool rebuild, int nToys) {
// get the distribution of lower limit
// if rebuild = false (default) it will re-use the results of the scan done
// for obtained the observed limit and no extra toys will be generated
// if rebuild a new set of B toys will be done and the procedure will be repeted
// for each toy
if (!rebuild) {
if (!fResults) {
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::GetLowerLimitDistribution(false) - result not existing\n";
return 0;
}
return fResults->GetLowerLimitDistribution();
}
TList * clsDist = 0;
TList * clsbDist = 0;
if (fUseCLs) clsDist = &fResults->fExpPValues;
else clsbDist = &fResults->fExpPValues;
return RebuildDistributions(false, nToys,clsDist, clsbDist);
}
SamplingDistribution * HypoTestInverter::GetUpperLimitDistribution(bool rebuild, int nToys) {
// get the distribution of lower limit
// if rebuild = false (default) it will re-use the results of the scan done
// for obtained the observed limit and no extra toys will be generated
// if rebuild a new set of B toys will be done and the procedure will be repeted
// for each toy
// The nuisance parameter value used for rebuild is the current one in the model
// so it is user responsability to set to the desired value (nomi
if (!rebuild) {
if (!fResults) {
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::GetUpperLimitDistribution(false) - result not existing\n";
return 0;
}
return fResults->GetUpperLimitDistribution();
}
TList * clsDist = 0;
TList * clsbDist = 0;
if (fUseCLs) clsDist = &fResults->fExpPValues;
else clsbDist = &fResults->fExpPValues;
return RebuildDistributions(true, nToys,clsDist, clsbDist);
}
void HypoTestInverter::SetData(RooAbsData & data) {
if (fCalculator0) fCalculator0->SetData(data);
}
SamplingDistribution * HypoTestInverter::RebuildDistributions(bool isUpper, int nToys, TList * clsDist, TList * clsbDist, TList * clbDist, const char *outputfile) {
// rebuild the sampling distributions by
// generating some toys and find for each of theam a new upper limit
// Return the upper limit distribution and optionally also the pValue distributions for Cls, Clsb and Clbxs
// as a TList for each scanned point
// The method uses the present parameter value. It is user responsability to give the current parameters to rebuild the distributions
// 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
// output file as an histogram
if (!fScannedVariable || !fCalculator0) return 0;
// get first background snapshot
const ModelConfig * bModel = fCalculator0->GetAlternateModel();
const ModelConfig * sbModel = fCalculator0->GetNullModel();
if (!bModel || ! sbModel) return 0;
RooArgSet paramPoint;
if (!sbModel->GetParametersOfInterest()) return 0;
paramPoint.add(*sbModel->GetParametersOfInterest());
const RooArgSet * poibkg = bModel->GetSnapshot();
if (!poibkg) {
oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RebuildDistribution - bakground snapshot not existing"
<< " assume is for POI = 0" << std::endl;
fScannedVariable->setVal(0);
paramPoint = RooArgSet(*fScannedVariable);
}
else
paramPoint = *poibkg;
// generate data at bkg parameter point
ToyMCSampler * toymcSampler = dynamic_cast<ToyMCSampler *>(fCalculator0->GetTestStatSampler() );
if (!toymcSampler) {
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter::RebuildDistribution - no toy MC sampler existing" << std::endl;
return 0;
}
// set up test stat sampler in case of asymptotic calculator
if (dynamic_cast<RooStats::AsymptoticCalculator*>(fCalculator0) ) {
toymcSampler->SetObservables(*sbModel->GetObservables() );
toymcSampler->SetParametersForTestStat(*sbModel->GetParametersOfInterest());
toymcSampler->SetPdf(*sbModel->GetPdf());
toymcSampler->SetNuisanceParameters(*sbModel->GetNuisanceParameters());
if (sbModel->GetGlobalObservables() ) toymcSampler->SetGlobalObservables(*sbModel->GetGlobalObservables() );
// set number of events
if (!sbModel->GetPdf()->canBeExtended())
toymcSampler->SetNEventsPerToy(1);
}
// loop on data to generate
int nPoints = fNBins;
bool storePValues = clsDist || clsbDist || clbDist;
if (fNBins <=0 && storePValues) {
oocoutW((TObject*)0,InputArguments) << "HypoTestInverter::RebuildDistribution - cannot return p values distribution with the auto scan" << std::endl;
storePValues = false;
nPoints = 0;
}
if (storePValues) {
if (fResults) nPoints = fResults->ArraySize();
if (nPoints <=0) {
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - result is not existing and number of point to scan is not set"
<< std::endl;
return 0;
}
}
if (nToys <= 0) nToys = 100; // default value
std::vector<std::vector<double> > CLs_values(nPoints);
std::vector<std::vector<double> > CLsb_values(nPoints);
std::vector<std::vector<double> > CLb_values(nPoints);
if (storePValues) {
for (int i = 0; i < nPoints; ++i) {
CLs_values[i].reserve(nToys);
CLb_values[i].reserve(nToys);
CLsb_values[i].reserve(nToys);
}
}
std::vector<double> limit_values; limit_values.reserve(nToys);
oocoutI((TObject*)0,InputArguments) << "HypoTestInverter - rebuilding the p value distributions by generating ntoys = "
<< nToys << std::endl;
oocoutI((TObject*)0,InputArguments) << "Rebuilding using parameter of interest point: ";
RooStats::PrintListContent(paramPoint, oocoutI((TObject*)0,InputArguments) );
if (sbModel->GetNuisanceParameters() ) {
oocoutI((TObject*)0,InputArguments) << "And using nuisance parameters: ";
RooStats::PrintListContent(*sbModel->GetNuisanceParameters(), oocoutI((TObject*)0,InputArguments) );
}
// save all parameters to restore them later
assert(bModel->GetPdf() );
assert(bModel->GetObservables() );
RooArgSet * allParams = bModel->GetPdf()->getParameters( *bModel->GetObservables() );
RooArgSet saveParams;
allParams->snapshot(saveParams);
TFile * fileOut = TFile::Open(outputfile,"RECREATE");
if (!fileOut) {
oocoutE((TObject*)0,InputArguments) << "HypoTestInverter - RebuildDistributions - Error opening file " << outputfile
<< " - the resulting limits will not be stored" << std::endl;
}
// create temporary histograms to store the limit result
TH1D * hL = new TH1D("lowerLimitDist","Rebuilt lower limit distribution",100,1.,0.);
TH1D * hU = new TH1D("upperLimitDist","Rebuilt upper limit distribution",100,1.,0.);
TH1D * hN = new TH1D("nObs","Observed events",100,1.,0.);
hL->SetBuffer(2*nToys);
hU->SetBuffer(2*nToys);
std::vector<TH1*> hCLb;
std::vector<TH1*> hCLsb;
std::vector<TH1*> hCLs;
if (storePValues) {
for (int i = 0; i < nPoints; ++i) {
hCLb.push_back(new TH1D(TString::Format("CLbDist_bin%d",i),"CLb distribution",100,1.,0.));
hCLs.push_back(new TH1D(TString::Format("ClsDist_bin%d",i),"CLs distribution",100,1.,0.));
hCLsb.push_back(new TH1D(TString::Format("CLsbDist_bin%d",i),"CLs+b distribution",100,1.,0.));
}
}
// loop now on the toys
for (int itoy = 0; itoy < nToys; ++itoy) {
oocoutP((TObject*)0,Eval) << "\nHypoTestInverter - RebuildDistributions - running toy # " << itoy << " / "
<< nToys << std::endl;
printf("\n\nshnapshot of s+b model \n");
sbModel->GetSnapshot()->Print("v");
// reset parameters to initial values to be sure in case they are not reset
if (itoy> 0) *allParams = saveParams;
// need to set th epdf to clear the cache in ToyMCSampler
// pdf we must use is background pdf
toymcSampler->SetPdf(*bModel->GetPdf() );
RooAbsData * bkgdata = toymcSampler->GenerateToyData(paramPoint);
double nObs = bkgdata->sumEntries();
// for debugging in case of number counting models
if (bkgdata->numEntries() ==1 && !bModel->GetPdf()->canBeExtended()) {
oocoutP((TObject*)0,Generation) << "Generate observables are : ";
RooArgList genObs(*bkgdata->get(0));
RooStats::PrintListContent(genObs, oocoutP((TObject*)0,Generation) );
nObs = 0;
for (int i = 0; i < genObs.getSize(); ++i) {
RooRealVar * x = dynamic_cast<RooRealVar*>(&genObs[i]);
if (x) nObs += x->getVal();
}
}
hN->Fill(nObs);
// by copying I will have the same min/max as previous ones
HypoTestInverter inverter = *this;
inverter.SetData(*bkgdata);
// print global observables
auto gobs = bModel->GetPdf()->getVariables()->selectCommon(* sbModel->GetGlobalObservables() );
gobs->Print("v");
HypoTestInverterResult * r = inverter.GetInterval();
if (r == 0) continue;
double value = (isUpper) ? r->UpperLimit() : r->LowerLimit();
limit_values.push_back( value );
hU->Fill(r->UpperLimit() );
hL->Fill(r->LowerLimit() );
std::cout << "The computed upper limit for toy #" << itoy << " is " << value << std::endl;
// write every 10 toys
if (itoy%10 == 0 || itoy == nToys-1) {
hU->Write("",TObject::kOverwrite);
hL->Write("",TObject::kOverwrite);
hN->Write("",TObject::kOverwrite);
}
if (!storePValues) continue;
if (nPoints < r->ArraySize()) {
oocoutW((TObject*)0,InputArguments) << "HypoTestInverter: skip extra points" << std::endl;
}
else if (nPoints > r->ArraySize()) {
oocoutW((TObject*)0,InputArguments) << "HypoTestInverter: missing some points" << std::endl;
}
for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
HypoTestResult * hr = r->GetResult(ipoint);
if (hr) {
CLs_values[ipoint].push_back( hr->CLs() );
CLsb_values[ipoint].push_back( hr->CLsplusb() );
CLb_values[ipoint].push_back( hr->CLb() );
hCLs[ipoint]->Fill( hr->CLs() );
hCLb[ipoint]->Fill( hr->CLb() );
hCLsb[ipoint]->Fill( hr->CLsplusb() );
}
else {
oocoutW((TObject*)0,InputArguments) << "HypoTestInverter: missing result for point: x = "
<< fResults->GetXValue(ipoint) << std::endl;
}
}
// write every 10 toys
if (itoy%10 == 0 || itoy == nToys-1) {
for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
hCLs[ipoint]->Write("",TObject::kOverwrite);
hCLb[ipoint]->Write("",TObject::kOverwrite);
hCLsb[ipoint]->Write("",TObject::kOverwrite);
}
}
delete r;
delete bkgdata;
}
if (storePValues) {
if (clsDist) clsDist->SetOwner(true);
if (clbDist) clbDist->SetOwner(true);
if (clsbDist) clsbDist->SetOwner(true);
oocoutI((TObject*)0,InputArguments) << "HypoTestInverter: storing rebuilt p values " << std::endl;
for (int ipoint = 0; ipoint < nPoints; ++ipoint) {
if (clsDist) {
TString name = TString::Format("CLs_distrib_%d",ipoint);
clsDist->Add( new SamplingDistribution(name,name,CLs_values[ipoint] ) );
}
if (clbDist) {
TString name = TString::Format("CLb_distrib_%d",ipoint);
clbDist->Add( new SamplingDistribution(name,name,CLb_values[ipoint] ) );
}
if (clsbDist) {
TString name = TString::Format("CLsb_distrib_%d",ipoint);
clsbDist->Add( new SamplingDistribution(name,name,CLsb_values[ipoint] ) );
}
}
}
if (fileOut) {
fileOut->Close();
}
else {
// delete all the histograms
delete hL;
delete hU;
for (int i = 0; i < nPoints && storePValues; ++i) {
delete hCLs[i];
delete hCLb[i];
delete hCLsb[i];
}
}
const char * disName = (isUpper) ? "upperLimit_dist" : "lowerLimit_dist";
return new SamplingDistribution(disName, disName, limit_values);
}