#ifndef ROOSTATS_ToyMCSampler
#define ROOSTATS_ToyMCSampler
#ifndef ROOT_Rtypes
#include "Rtypes.h"
#endif
#include <vector>
#include <sstream>
#include "RooStats/TestStatSampler.h"
#include "RooStats/SamplingDistribution.h"
#include "RooStats/TestStatistic.h"
#include "RooStats/ModelConfig.h"
#include "RooStats/ProofConfig.h"
#include "RooWorkspace.h"
#include "RooMsgService.h"
#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
namespace RooStats {
class DetailedOutputAggregator;
class NuisanceParametersSampler {
public:
NuisanceParametersSampler(RooAbsPdf *prior=NULL, const RooArgSet *parameters=NULL, Int_t nToys=1000, Bool_t asimov=kFALSE) :
fPrior(prior),
fParams(parameters),
fNToys(nToys),
fExpected(asimov),
fPoints(NULL),
fIndex(0)
{
if(prior) Refresh();
}
virtual ~NuisanceParametersSampler() {
if(fPoints) { delete fPoints; fPoints = NULL; }
}
void NextPoint(RooArgSet& nuisPoint, Double_t& weight);
protected:
void Refresh();
private:
RooAbsPdf *fPrior;
const RooArgSet *fParams;
Int_t fNToys;
Bool_t fExpected;
RooAbsData *fPoints;
Int_t fIndex;
};
class ToyMCSampler: public TestStatSampler {
public:
ToyMCSampler();
ToyMCSampler(TestStatistic &ts, Int_t ntoys);
virtual ~ToyMCSampler();
static void SetAlwaysUseMultiGen(Bool_t flag);
void SetUseMultiGen(Bool_t flag) { fUseMultiGen = flag ; }
virtual SamplingDistribution* GetSamplingDistribution(RooArgSet& paramPoint);
virtual RooDataSet* GetSamplingDistributions(RooArgSet& paramPoint);
virtual RooDataSet* GetSamplingDistributionsSingleWorker(RooArgSet& paramPoint);
virtual SamplingDistribution* AppendSamplingDistribution(
RooArgSet& allParameters,
SamplingDistribution* last,
Int_t additionalMC
);
virtual void AddTestStatistic(TestStatistic* t = NULL) {
if( t == NULL ) {
oocoutI((TObject*)0,InputArguments) << "No test statistic given. Doing nothing." << std::endl;
return;
}
fTestStatistics.push_back( t );
}
virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, RooAbsPdf& pdf) const {
if(fExpectedNuisancePar) oocoutE((TObject*)NULL,InputArguments) << "ToyMCSampler: using expected nuisance parameters but ignoring weight. Use GetSamplingDistribution(paramPoint, weight) instead." << std::endl;
double weight;
return GenerateToyData(paramPoint, weight, pdf);
}
virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint) const { return GenerateToyData(paramPoint,*fPdf); }
virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight, RooAbsPdf& pdf) const;
virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight) const { return GenerateToyData(paramPoint,weight,*fPdf); }
virtual void GenerateGlobalObservables(RooAbsPdf& pdf) const;
virtual Double_t EvaluateTestStatistic(RooAbsData& data, RooArgSet& nullPOI, int i ) {
return fTestStatistics[i]->Evaluate(data, nullPOI);
}
virtual Double_t EvaluateTestStatistic(RooAbsData& data, RooArgSet& nullPOI) { return EvaluateTestStatistic( data,nullPOI, 0 ); }
virtual RooArgList* EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi);
virtual TestStatistic* GetTestStatistic(unsigned int i) const {
if( fTestStatistics.size() <= i ) return NULL;
return fTestStatistics[i];
}
virtual TestStatistic* GetTestStatistic(void) const { return GetTestStatistic(0); }
virtual Double_t ConfidenceLevel() const { return 1. - fSize; }
virtual void Initialize(
RooAbsArg& ,
RooArgSet& ,
RooArgSet&
) {}
virtual Int_t GetNToys(void) { return fNToys; }
virtual void SetNToys(const Int_t ntoy) { fNToys = ntoy; }
virtual void SetNEventsPerToy(const Int_t nevents) {
fNEvents = nevents;
}
virtual void SetParametersForTestStat(const RooArgSet& nullpoi) {
if( fParametersForTestStat ) delete fParametersForTestStat;
fParametersForTestStat = (const RooArgSet*)nullpoi.snapshot();
}
virtual void SetPdf(RooAbsPdf& pdf) { fPdf = &pdf; ClearCache(); }
virtual void SetPriorNuisance(RooAbsPdf* pdf) {
fPriorNuisance = pdf;
if (fNuisanceParametersSampler) {
delete fNuisanceParametersSampler;
fNuisanceParametersSampler = NULL;
}
}
virtual void SetNuisanceParameters(const RooArgSet& np) { fNuisancePars = &np; }
virtual void SetObservables(const RooArgSet& o) { fObservables = &o; }
virtual void SetGlobalObservables(const RooArgSet& o) { fGlobalObservables = &o; }
virtual void SetTestSize(Double_t size) { fSize = size; }
virtual void SetConfidenceLevel(Double_t cl) { fSize = 1. - cl; }
virtual void SetTestStatistic(TestStatistic *testStatistic, unsigned int i) {
if( fTestStatistics.size() < i ) {
oocoutE((TObject*)NULL,InputArguments) << "Cannot set test statistic for this index." << std::endl;
return;
}
if( fTestStatistics.size() == i)
fTestStatistics.push_back(testStatistic);
else
fTestStatistics[i] = testStatistic;
}
virtual void SetTestStatistic(TestStatistic *t) { return SetTestStatistic(t,0); }
virtual void SetExpectedNuisancePar(Bool_t i = kTRUE) { fExpectedNuisancePar = i; }
virtual void SetAsimovNuisancePar(Bool_t i = kTRUE) { fExpectedNuisancePar = i; }
Bool_t CheckConfig(void);
void SetGenerateBinned(bool binned = true) { fGenerateBinned = binned; }
void SetGenerateBinnedTag( const char* binnedTag = "" ) { fGenerateBinnedTag = binnedTag; }
void SetGenerateAutoBinned( Bool_t autoBinned = kTRUE ) { fGenerateAutoBinned = autoBinned; }
void SetSamplingDistName(const char* name) { if(name) fSamplingDistName = name; }
std::string GetSamplingDistName(void) { return fSamplingDistName; }
void SetMaxToys(Double_t t) { fMaxToys = t; }
void SetToysLeftTail(Double_t toys, Double_t threshold) {
fToysInTails = toys;
fAdaptiveLowLimit = threshold;
fAdaptiveHighLimit = RooNumber::infinity();
}
void SetToysRightTail(Double_t toys, Double_t threshold) {
fToysInTails = toys;
fAdaptiveHighLimit = threshold;
fAdaptiveLowLimit = -RooNumber::infinity();
}
void SetToysBothTails(Double_t toys, Double_t low_threshold, Double_t high_threshold) {
fToysInTails = toys;
fAdaptiveHighLimit = high_threshold;
fAdaptiveLowLimit = low_threshold;
}
void SetProofConfig(ProofConfig *pc = NULL) { fProofConfig = pc; }
void SetProtoData(const RooDataSet* d) { fProtoData = d; }
protected:
const RooArgList* EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi, DetailedOutputAggregator& detOutAgg);
RooAbsData* Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooDataSet *protoData=NULL, int forceEvents=0) const;
virtual void ClearCache();
RooAbsPdf *fPdf;
const RooArgSet* fParametersForTestStat;
std::vector<TestStatistic*> fTestStatistics;
std::string fSamplingDistName;
RooAbsPdf *fPriorNuisance;
const RooArgSet *fNuisancePars;
const RooArgSet *fObservables;
const RooArgSet *fGlobalObservables;
Int_t fNToys;
Int_t fNEvents;
Double_t fSize;
Bool_t fExpectedNuisancePar;
Bool_t fGenerateBinned;
TString fGenerateBinnedTag;
Bool_t fGenerateAutoBinned;
Double_t fToysInTails;
Double_t fMaxToys;
Double_t fAdaptiveLowLimit;
Double_t fAdaptiveHighLimit;
const RooDataSet *fProtoData;
ProofConfig *fProofConfig;
mutable NuisanceParametersSampler *fNuisanceParametersSampler;
mutable RooArgSet* _allVars ;
mutable std::list<RooAbsPdf*> _pdfList ;
mutable std::list<RooArgSet*> _obsList ;
mutable std::list<RooAbsPdf::GenSpec*> _gsList ;
mutable RooAbsPdf::GenSpec* _gs1 ;
mutable RooAbsPdf::GenSpec* _gs2 ;
mutable RooAbsPdf::GenSpec* _gs3 ;
mutable RooAbsPdf::GenSpec* _gs4 ;
static Bool_t fgAlwaysUseMultiGen ;
Bool_t fUseMultiGen ;
protected:
ClassDef(ToyMCSampler,3)
};
}
#endif