#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 NuisanceParametersSampler;
class ToyMCSampler: public TestStatSampler {
public:
ToyMCSampler() :
fTestStat(NULL), fSamplingDistName("temp"), fNToys(1)
{
fPdf = NULL;
fPriorNuisance = NULL;
fNullPOI = NULL;
fNuisancePars = NULL;
fObservables = NULL;
fGlobalObservables = NULL;
fSize = 0.05;
fNEvents = 0;
fGenerateBinned = kFALSE;
fExpectedNuisancePar = kFALSE;
fToysInTails = 0.0;
fMaxToys = RooNumber::infinity();
fAdaptiveLowLimit = -RooNumber::infinity();
fAdaptiveHighLimit = RooNumber::infinity();
fImportanceDensity = NULL;
fImportanceSnapshot = NULL;
fProtoData = NULL;
fProofConfig = NULL;
fNuisanceParametersSampler = NULL;
_allVars = NULL ;
_gs1 = NULL ;
_gs2 = NULL ;
_gs3 = NULL ;
_gs4 = NULL ;
fUseMultiGen = kFALSE ;
}
ToyMCSampler(TestStatistic &ts, Int_t ntoys) :
fTestStat(&ts), fSamplingDistName(ts.GetVarName()), fNToys(ntoys)
{
fPdf = NULL;
fPriorNuisance = NULL;
fNullPOI = NULL;
fNuisancePars = NULL;
fObservables = NULL;
fGlobalObservables = NULL;
fSize = 0.05;
fNEvents = 0;
fGenerateBinned = kFALSE;
fExpectedNuisancePar = kFALSE;
fToysInTails = 0.0;
fMaxToys = RooNumber::infinity();
fAdaptiveLowLimit = -RooNumber::infinity();
fAdaptiveHighLimit = RooNumber::infinity();
fImportanceDensity = NULL;
fImportanceSnapshot = NULL;
fProtoData = NULL;
fProofConfig = NULL;
fNuisanceParametersSampler = NULL;
_allVars = NULL ;
_gs1 = NULL ;
_gs2 = NULL ;
_gs3 = NULL ;
_gs4 = NULL ;
fUseMultiGen = kFALSE ;
}
virtual ~ToyMCSampler();
static void SetAlwaysUseMultiGen(Bool_t flag) { fgAlwaysUseMultiGen = flag ; }
void SetUseMultiGen(Bool_t flag) { fUseMultiGen = flag ; }
virtual SamplingDistribution* GetSamplingDistribution(RooArgSet& paramPoint);
virtual SamplingDistribution* GetSamplingDistributionSingleWorker(RooArgSet& paramPoint);
virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint) const {
if(fExpectedNuisancePar) oocoutE((TObject*)NULL,InputArguments) << "ToyMCSampler: using expected nuisance parameters but ignoring weight. Use GetSamplingDistribution(paramPoint, weight) instead." << endl;
double weight;
return GenerateToyData(paramPoint, weight);
}
virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight) const;
virtual RooAbsData* GenerateToyDataImportanceSampling(RooArgSet& paramPoint, double& weight) const;
virtual void GenerateGlobalObservables(void) const;
virtual SamplingDistribution* 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;
}
virtual Double_t EvaluateTestStatistic(RooAbsData& data, RooArgSet& nullPOI) {
return fTestStat->Evaluate(data, nullPOI);
}
virtual TestStatistic* GetTestStatistic() const { return fTestStat; }
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 (fNullPOI) delete fNullPOI; fNullPOI = (RooArgSet*)nullpoi.snapshot();
}
virtual void SetPdf(RooAbsPdf& pdf) { fPdf = &pdf; }
virtual void SetPriorNuisance(RooAbsPdf* pdf) { fPriorNuisance = pdf; }
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) { fTestStat = testStatistic; }
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 SetSamplingDistName(const char* name) { if(name) fSamplingDistName = name; }
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 SetImportanceDensity(RooAbsPdf *p) {
if(p) oocoutW((TObject*)NULL,InputArguments) << "ToyMCSampler Importance Sampling: This is in beta." << endl;
fImportanceDensity = p;
}
void SetImportanceSnapshot(const RooArgSet &s) { fImportanceSnapshot = &s; }
void SetProofConfig(ProofConfig *pc = NULL) { fProofConfig = pc; }
void SetProtoData(const RooDataSet* d) { fProtoData = d; }
protected:
RooAbsData* Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooDataSet *protoData=NULL, int forceEvents=0) const;
TestStatistic *fTestStat;
RooAbsPdf *fPdf;
string fSamplingDistName;
RooAbsPdf *fPriorNuisance;
RooArgSet *fNullPOI;
const RooArgSet *fNuisancePars;
const RooArgSet *fObservables;
const RooArgSet *fGlobalObservables;
Int_t fNToys;
Int_t fNEvents;
Double_t fSize;
Bool_t fExpectedNuisancePar;
Bool_t fGenerateBinned;
Double_t fToysInTails;
Double_t fMaxToys;
Double_t fAdaptiveLowLimit;
Double_t fAdaptiveHighLimit;
RooAbsPdf *fImportanceDensity;
const RooArgSet *fImportanceSnapshot;
const RooDataSet *fProtoData;
ProofConfig *fProofConfig;
mutable NuisanceParametersSampler *fNuisanceParametersSampler;
mutable RooArgSet* _allVars ;
mutable list<RooAbsPdf*> _pdfList ;
mutable list<RooArgSet*> _obsList ;
mutable 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,2)
};
}
#endif