#include "RooStats/ToyMCSampler.h"
#ifndef ROO_MSG_SERVICE
#include "RooMsgService.h"
#endif
#ifndef ROO_DATA_HIST
#include "RooDataHist.h"
#endif
#ifndef ROO_REAL_VAR
#include "RooRealVar.h"
#endif
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooRandom.h"
#include "RooStudyManager.h"
#include "RooStats/ToyMCStudy.h"
#include "RooStats/DetailedOutputAggregator.h"
#include "RooStats/RooStatsUtils.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "TMath.h"
using namespace RooFit;
using namespace std;
ClassImp(RooStats::ToyMCSampler)
namespace RooStats {
void NuisanceParametersSampler::NextPoint(RooArgSet& nuisPoint, Double_t& weight) {
if (fIndex >= fNToys) {
Refresh();
fIndex = 0;
}
nuisPoint = *fPoints->get(fIndex++);
weight = fPoints->weight();
if(fPoints->weight() == 0.0) {
oocoutI((TObject*)NULL,Generation) << "Weight 0 encountered. Skipping." << endl;
NextPoint(nuisPoint, weight);
}
}
void NuisanceParametersSampler::Refresh() {
if (!fPrior || !fParams) return;
if (fPoints) delete fPoints;
if (fExpected) {
oocoutI((TObject*)NULL,InputArguments) << "Using expected nuisance parameters." << endl;
int nBins = fNToys;
TIter it2 = fParams->createIterator();
RooRealVar *myarg2;
while ((myarg2 = dynamic_cast<RooRealVar*>(it2.Next()))) {
myarg2->setBins(nBins);
}
fPoints = fPrior->generate(
*fParams,
AllBinned(),
ExpectedData(),
NumEvents(1)
);
if(fPoints->numEntries() != fNToys) {
fNToys = fPoints->numEntries();
oocoutI((TObject*)NULL,InputArguments) <<
"Adjusted number of toys to number of bins of nuisance parameters: " << fNToys << endl;
}
}else{
oocoutI((TObject*)NULL,InputArguments) << "Using randomized nuisance parameters." << endl;
fPoints = fPrior->generate(*fParams, fNToys);
}
}
Bool_t ToyMCSampler::fgAlwaysUseMultiGen = kFALSE ;
void ToyMCSampler::SetAlwaysUseMultiGen(Bool_t flag) { fgAlwaysUseMultiGen = flag ; }
ToyMCSampler::ToyMCSampler() : fSamplingDistName("SD"), fNToys(1)
{
fPdf = NULL;
fParametersForTestStat = NULL;
fPriorNuisance = NULL;
fNuisancePars = NULL;
fObservables = NULL;
fGlobalObservables = NULL;
fSize = 0.05;
fNEvents = 0;
fGenerateBinned = kFALSE;
fGenerateBinnedTag = "";
fGenerateAutoBinned = kTRUE;
fExpectedNuisancePar = kFALSE;
fToysInTails = 0.0;
fMaxToys = RooNumber::infinity();
fAdaptiveLowLimit = -RooNumber::infinity();
fAdaptiveHighLimit = RooNumber::infinity();
fProtoData = NULL;
fProofConfig = NULL;
fNuisanceParametersSampler = NULL;
_allVars = NULL ;
_gs1 = NULL ;
_gs2 = NULL ;
_gs3 = NULL ;
_gs4 = NULL ;
RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);
fUseMultiGen = kFALSE ;
}
ToyMCSampler::ToyMCSampler(TestStatistic &ts, Int_t ntoys) :
fSamplingDistName(ts.GetVarName()), fNToys(ntoys)
{
fPdf = NULL;
fParametersForTestStat = NULL;
fPriorNuisance = NULL;
fNuisancePars = NULL;
fObservables = NULL;
fGlobalObservables = NULL;
fSize = 0.05;
fNEvents = 0;
fGenerateBinned = kFALSE;
fGenerateBinnedTag = "";
fGenerateAutoBinned = kTRUE;
fExpectedNuisancePar = kFALSE;
fToysInTails = 0.0;
fMaxToys = RooNumber::infinity();
fAdaptiveLowLimit = -RooNumber::infinity();
fAdaptiveHighLimit = RooNumber::infinity();
fProtoData = NULL;
fProofConfig = NULL;
fNuisanceParametersSampler = NULL;
_allVars = NULL ;
_gs1 = NULL ;
_gs2 = NULL ;
_gs3 = NULL ;
_gs4 = NULL ;
RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);
fUseMultiGen = kFALSE ;
AddTestStatistic(&ts);
}
ToyMCSampler::~ToyMCSampler() {
if(fNuisanceParametersSampler) delete fNuisanceParametersSampler;
ClearCache();
}
Bool_t ToyMCSampler::CheckConfig(void) {
bool goodConfig = true;
if(fTestStatistics.size() == 0 || fTestStatistics[0] == NULL) { ooccoutE((TObject*)NULL,InputArguments) << "Test statistic not set." << endl; goodConfig = false; }
if(!fObservables) { ooccoutE((TObject*)NULL,InputArguments) << "Observables not set." << endl; goodConfig = false; }
if(!fParametersForTestStat) { ooccoutE((TObject*)NULL,InputArguments) << "Parameter values used to evaluate the test statistic are not set." << endl; goodConfig = false; }
if(!fPdf) { ooccoutE((TObject*)NULL,InputArguments) << "Pdf not set." << endl; goodConfig = false; }
return goodConfig;
}
RooArgList* ToyMCSampler::EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi) {
DetailedOutputAggregator detOutAgg;
const RooArgList* allTS = EvaluateAllTestStatistics(data, poi, detOutAgg);
if (!allTS) return 0;
return dynamic_cast<RooArgList*>(allTS->snapshot());
}
const RooArgList* ToyMCSampler::EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi, DetailedOutputAggregator& detOutAgg) {
RooArgSet *allVars = fPdf ? fPdf->getVariables() : 0;
RooArgSet *saveAll = allVars ? dynamic_cast<RooArgSet*>(allVars->snapshot()) : 0;
for( unsigned int i = 0; i < fTestStatistics.size(); i++ ) {
if( fTestStatistics[i] == NULL ) continue;
TString name( TString::Format("%s_TS%u", fSamplingDistName.c_str(), i) );
RooArgSet* parForTS = dynamic_cast<RooArgSet*>(poi.snapshot());
RooRealVar ts( name, fTestStatistics[i]->GetVarName(), fTestStatistics[i]->Evaluate( data, *parForTS ) );
RooArgList tset(ts);
detOutAgg.AppendArgSet(&tset);
delete parForTS;
if (const RooArgSet* detOut = fTestStatistics[i]->GetDetailedOutput()) {
name.Append("_");
detOutAgg.AppendArgSet(detOut, name);
}
if (saveAll) *allVars = *saveAll;
}
delete saveAll;
delete allVars;
return detOutAgg.GetAsArgList();
}
SamplingDistribution* ToyMCSampler::GetSamplingDistribution(RooArgSet& paramPointIn) {
if(fTestStatistics.size() > 1) {
oocoutW((TObject*)NULL, InputArguments) << "Multiple test statistics defined, but only one distribution will be returned." << endl;
for( unsigned int i=0; i < fTestStatistics.size(); i++ ) {
oocoutW((TObject*)NULL, InputArguments) << " \t test statistic: " << fTestStatistics[i] << endl;
}
}
RooDataSet* r = GetSamplingDistributions(paramPointIn);
if(r == NULL || r->numEntries() == 0) {
oocoutW((TObject*)NULL, Generation) << "no sampling distribution generated" << endl;
return NULL;
}
SamplingDistribution* samp = new SamplingDistribution( r->GetName(), r->GetTitle(), *r );
delete r;
return samp;
}
RooDataSet* ToyMCSampler::GetSamplingDistributions(RooArgSet& paramPointIn)
{
if(!fProofConfig)
return GetSamplingDistributionsSingleWorker(paramPointIn);
CheckConfig();
if(fToysInTails) {
fToysInTails = 0;
oocoutW((TObject*)NULL, InputArguments)
<< "Adaptive sampling in ToyMCSampler is not supported for parallel runs."
<< endl;
}
Int_t totToys = fNToys;
fNToys = (int)ceil((double)fNToys / (double)fProofConfig->GetNExperiments());
ToyMCStudy* toymcstudy = new ToyMCStudy ;
toymcstudy->SetToyMCSampler(*this);
toymcstudy->SetParamPoint(paramPointIn);
toymcstudy->SetRandomSeed(RooRandom::randomGenerator()->Integer(TMath::Limits<unsigned int>::Max() ) );
RooWorkspace w(fProofConfig->GetWorkspace());
RooStudyManager studymanager(w, *toymcstudy);
studymanager.runProof(fProofConfig->GetNExperiments(), fProofConfig->GetHost(), fProofConfig->GetShowGui());
RooDataSet* output = toymcstudy->merge();
fNToys = totToys;
delete toymcstudy;
return output;
}
RooDataSet* ToyMCSampler::GetSamplingDistributionsSingleWorker(RooArgSet& paramPointIn)
{
ClearCache();
CheckConfig();
RooArgSet *paramPoint = (RooArgSet*) paramPointIn.snapshot();
RooArgSet *allVars = fPdf->getVariables();
RooArgSet *saveAll = (RooArgSet*) allVars->snapshot();
DetailedOutputAggregator detOutAgg;
Double_t toysInTails = 0.0;
for (Int_t i = 0; i < fMaxToys; ++i) {
if (toysInTails >= fToysInTails && i+1 > fNToys) break;
if ( i% 500 == 0 && i>0 ) {
oocoutP((TObject*)0,Generation) << "generated toys: " << i << " / " << fNToys;
if (fToysInTails) ooccoutP((TObject*)0,Generation) << " (tails: " << toysInTails << " / " << fToysInTails << ")" << std::endl;
else ooccoutP((TObject*)0,Generation) << endl;
}
Double_t valueFirst = -999.0, weight = 1.0;
*allVars = *saveAll;
RooAbsData* toydata = GenerateToyData(*paramPoint, weight);
*allVars = *fParametersForTestStat;
const RooArgList* allTS = EvaluateAllTestStatistics(*toydata, *fParametersForTestStat, detOutAgg);
if (allTS->getSize() > Int_t(fTestStatistics.size()))
detOutAgg.AppendArgSet( fGlobalObservables, "globObs_" );
if (RooRealVar* firstTS = dynamic_cast<RooRealVar*>(allTS->first()))
valueFirst = firstTS->getVal();
delete toydata;
if(valueFirst != valueFirst) {
oocoutW((TObject*)NULL, Generation) << "skip: " << valueFirst << ", " << weight << endl;
continue;
}
detOutAgg.CommitSet(weight);
if (valueFirst <= fAdaptiveLowLimit || valueFirst >= fAdaptiveHighLimit) {
if(weight >= 0.) toysInTails += weight;
else toysInTails += 1.;
}
}
*allVars = *saveAll;
delete saveAll;
delete allVars;
delete paramPoint;
return detOutAgg.GetAsDataSet(fSamplingDistName, fSamplingDistName);
}
void ToyMCSampler::GenerateGlobalObservables(RooAbsPdf& pdf) const {
if(!fGlobalObservables || fGlobalObservables->getSize()==0) {
ooccoutE((TObject*)NULL,InputArguments) << "Global Observables not set." << endl;
return;
}
if (fUseMultiGen || fgAlwaysUseMultiGen) {
RooSimultaneous* simPdf = dynamic_cast<RooSimultaneous*>( &pdf );
if (!simPdf) {
RooDataSet *one = pdf.generate(*fGlobalObservables, 1);
const RooArgSet *values = one->get(0);
if (!_allVars) {
_allVars = pdf.getVariables();
}
*_allVars = *values;
delete one;
} else {
if (_pdfList.size() == 0) {
RooCategory& channelCat = (RooCategory&)simPdf->indexCat();
int nCat = channelCat.numTypes();
for (int i=0; i < nCat; ++i){
channelCat.setIndex(i);
RooAbsPdf* pdftmp = simPdf->getPdf(channelCat.getLabel());
assert(pdftmp);
RooArgSet* globtmp = pdftmp->getObservables(*fGlobalObservables);
RooAbsPdf::GenSpec* gs = pdftmp->prepareMultiGen(*globtmp, NumEvents(1));
_pdfList.push_back(pdftmp);
_obsList.push_back(globtmp);
_gsList.push_back(gs);
}
}
list<RooArgSet*>::iterator oiter = _obsList.begin();
list<RooAbsPdf::GenSpec*>::iterator giter = _gsList.begin();
for (list<RooAbsPdf*>::iterator iter = _pdfList.begin(); iter != _pdfList.end(); ++iter, ++giter, ++oiter) {
RooDataSet* tmp = (*iter)->generate(**giter);
**oiter = *tmp->get(0);
delete tmp;
}
}
} else {
RooDataSet* one = pdf.generateSimGlobal( *fGlobalObservables, 1 );
const RooArgSet *values = one->get(0);
RooArgSet* allVars = pdf.getVariables();
*allVars = *values;
delete allVars;
delete one;
}
}
RooAbsData* ToyMCSampler::GenerateToyData(RooArgSet& paramPoint, double& weight, RooAbsPdf& pdf) const {
if(!fObservables) {
ooccoutE((TObject*)NULL,InputArguments) << "Observables not set." << endl;
return NULL;
}
RooArgSet* allVars = fPdf->getVariables();
*allVars = paramPoint;
if(!fNuisanceParametersSampler && fPriorNuisance && fNuisancePars) {
fNuisanceParametersSampler = new NuisanceParametersSampler(fPriorNuisance, fNuisancePars, fNToys, fExpectedNuisancePar);
if ((fUseMultiGen || fgAlwaysUseMultiGen) && fNuisanceParametersSampler )
oocoutI((TObject*)NULL,InputArguments) << "Cannot use multigen when nuisance parameters vary for every toy" << endl;
}
RooArgSet observables(*fObservables);
if(fGlobalObservables && fGlobalObservables->getSize()) {
observables.remove(*fGlobalObservables);
GenerateGlobalObservables(pdf);
}
const RooArgSet* saveVars = (const RooArgSet*)allVars->snapshot();
if(fNuisanceParametersSampler) {
RooArgSet allVarsMinusParamPoint(*allVars);
allVarsMinusParamPoint.remove(paramPoint, kFALSE, kTRUE);
fNuisanceParametersSampler->NextPoint(allVarsMinusParamPoint, weight);
}else{
weight = 1.0;
}
RooAbsData *data = Generate(pdf, observables);
*allVars = *saveVars;
delete allVars;
delete saveVars;
return data;
}
RooAbsData* ToyMCSampler::Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooDataSet* protoData, int forceEvents) const {
if(fProtoData) {
protoData = fProtoData;
forceEvents = protoData->numEntries();
}
RooAbsData *data = NULL;
int events = forceEvents;
if(events == 0) events = fNEvents;
bool useMultiGen = (fUseMultiGen || fgAlwaysUseMultiGen) && !fNuisanceParametersSampler;
if(events == 0) {
if( pdf.canBeExtended() && pdf.expectedEvents(observables) > 0) {
if(fGenerateBinned) {
if(protoData) data = pdf.generate(observables, AllBinned(), Extended(), ProtoData(*protoData, true, true));
else data = pdf.generate(observables, AllBinned(), Extended());
}else{
if(protoData) {
if (useMultiGen) {
if (!_gs2) { _gs2 = pdf.prepareMultiGen(observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true)) ; }
data = pdf.generate(*_gs2) ;
} else {
data = pdf.generate (observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true));
}
}
else {
if (useMultiGen) {
if (!_gs1) { _gs1 = pdf.prepareMultiGen(observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag) ) ; }
data = pdf.generate(*_gs1) ;
} else {
data = pdf.generate (observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag) );
}
}
}
}else{
oocoutE((TObject*)0,InputArguments)
<< "ToyMCSampler: Error : pdf is not extended and number of events per toy is zero"
<< endl;
}
}else{
if(fGenerateBinned) {
if(protoData) data = pdf.generate(observables, events, AllBinned(), ProtoData(*protoData, true, true));
else data = pdf.generate(observables, events, AllBinned());
}else{
if(protoData) {
if (useMultiGen) {
if (!_gs3) { _gs3 = pdf.prepareMultiGen(observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true)); }
data = pdf.generate(*_gs3) ;
} else {
data = pdf.generate (observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true));
}
} else {
if (useMultiGen) {
if (!_gs4) { _gs4 = pdf.prepareMultiGen(observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag)); }
data = pdf.generate(*_gs4) ;
} else {
data = pdf.generate (observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag));
}
}
}
}
return data;
}
SamplingDistribution* ToyMCSampler::AppendSamplingDistribution(
RooArgSet& allParameters,
SamplingDistribution* last,
Int_t additionalMC)
{
Int_t tmp = fNToys;
fNToys = additionalMC;
SamplingDistribution* newSamples = GetSamplingDistribution(allParameters);
fNToys = tmp;
if(last){
last->Add(newSamples);
delete newSamples;
return last;
}
return newSamples;
}
void ToyMCSampler::ClearCache() {
if (_gs1) delete _gs1;
_gs1 = 0;
if (_gs2) delete _gs2;
_gs2 = 0;
if (_gs3) delete _gs3;
_gs3 = 0;
if (_gs4) delete _gs4;
_gs4 = 0;
if (_pdfList.size() > 0) {
std::list<RooArgSet*>::iterator oiter = _obsList.begin();
for (std::list<RooAbsPdf::GenSpec*>::iterator giter = _gsList.begin(); giter != _gsList.end(); ++giter, ++oiter) {
delete *oiter;
delete *giter;
}
_pdfList.clear();
_obsList.clear();
_gsList.clear();
}
if (_allVars) delete _allVars;
_allVars = 0;
}
}