Logo ROOT   6.08/07
Reference Guide
ToyMCSampler.h
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Sven Kreiss and Kyle Cranmer June 2010
3 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
4 // Additions and modifications by Mario Pelliccioni
5 /*************************************************************************
6  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
7  * All rights reserved. *
8  * *
9  * For the licensing terms see $ROOTSYS/LICENSE. *
10  * For the list of contributors see $ROOTSYS/README/CREDITS. *
11  *************************************************************************/
12 
13 #ifndef ROOSTATS_ToyMCSampler
14 #define ROOSTATS_ToyMCSampler
15 
16 
17 #ifndef ROOT_Rtypes
18 #include "Rtypes.h"
19 #endif
20 
21 #include <vector>
22 #include <sstream>
23 
26 #include "RooStats/TestStatistic.h"
27 #include "RooStats/ModelConfig.h"
28 #include "RooStats/ProofConfig.h"
29 
30 #include "RooWorkspace.h"
31 #include "RooMsgService.h"
32 #include "RooAbsPdf.h"
33 #include "RooRealVar.h"
34 
35 #include "RooDataSet.h"
36 
37 
38 
39 namespace RooStats {
40 
41  class DetailedOutputAggregator;
42 
43 
44 /**
45  Helper class for ToyMCSampler. Handles all of the nuisance parameter related
46  functions. Once instantiated, it gives a new nuisance parameter point
47  at each call to nextPoint(...).
48  Only used inside ToyMCSampler, ie "private" in the cxx file
49 */
50 
52 
53  public:
54  NuisanceParametersSampler(RooAbsPdf *prior=NULL, const RooArgSet *parameters=NULL, Int_t nToys=1000, Bool_t asimov=kFALSE) :
55  fPrior(prior),
56  fParams(parameters),
57  fNToys(nToys),
58  fExpected(asimov),
59  fPoints(NULL),
60  fIndex(0)
61  {
62  if(prior) Refresh();
63  }
65  if(fPoints) { delete fPoints; fPoints = NULL; }
66  }
67 
68  void NextPoint(RooArgSet& nuisPoint, Double_t& weight);
69 
70  protected:
71  void Refresh();
72 
73  private:
74  RooAbsPdf *fPrior; // prior for nuisance parameters
75  const RooArgSet *fParams; // nuisance parameters
78 
79  RooAbsData *fPoints; // generated nuisance parameter points
80  Int_t fIndex; // current index in fPoints array
81 };
82 
83 
84 /**
85 
86 ToyMCSampler is an implementation of the TestStatSampler interface.
87 It generates Toy Monte Carlo for a given parameter point and evaluates a
88 TestStatistic.
89 
90 For parallel runs, ToyMCSampler can be given an instance of ProofConfig
91 and then run in parallel using proof or proof-lite. Internally, it uses
92 ToyMCStudy with the RooStudyManager.
93 
94 \ingroup Roostats
95 
96 */
97 
98 
100 
101  public:
102 
103  ToyMCSampler();
104  ToyMCSampler(TestStatistic &ts, Int_t ntoys);
105  virtual ~ToyMCSampler();
106 
107  static void SetAlwaysUseMultiGen(Bool_t flag);
108 
109  void SetUseMultiGen(Bool_t flag) { fUseMultiGen = flag ; }
110 
111  // main interface
112  virtual SamplingDistribution* GetSamplingDistribution(RooArgSet& paramPoint);
113  virtual RooDataSet* GetSamplingDistributions(RooArgSet& paramPoint);
114  virtual RooDataSet* GetSamplingDistributionsSingleWorker(RooArgSet& paramPoint);
115 
116  virtual SamplingDistribution* AppendSamplingDistribution(
117  RooArgSet& allParameters,
118  SamplingDistribution* last,
119  Int_t additionalMC
120  );
121 
122 
123  // The pdf can be NULL in which case the density from SetPdf()
124  // is used. The snapshot and TestStatistic is also optional.
125  virtual void AddTestStatistic(TestStatistic* t = NULL) {
126  if( t == NULL ) {
127  oocoutI((TObject*)0,InputArguments) << "No test statistic given. Doing nothing." << std::endl;
128  return;
129  }
130 
131  //if( t == NULL && fTestStatistics.size() >= 1 ) t = fTestStatistics[0];
132 
133  fTestStatistics.push_back( t );
134  }
135 
136 
137 
138  // generates toy data
139  // without weight
140  virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, RooAbsPdf& pdf) const {
141  if(fExpectedNuisancePar) oocoutE((TObject*)NULL,InputArguments) << "ToyMCSampler: using expected nuisance parameters but ignoring weight. Use GetSamplingDistribution(paramPoint, weight) instead." << std::endl;
142  double weight;
143  return GenerateToyData(paramPoint, weight, pdf);
144  }
145  virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint) const { return GenerateToyData(paramPoint,*fPdf); }
146  // with weight
147  virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight, RooAbsPdf& pdf) const;
148  virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight) const { return GenerateToyData(paramPoint,weight,*fPdf); }
149 
150  // generate global observables
151  virtual void GenerateGlobalObservables(RooAbsPdf& pdf) const;
152 
153 
154  // Main interface to evaluate the test statistic on a dataset
156  return fTestStatistics[i]->Evaluate(data, nullPOI);
157  }
158  virtual Double_t EvaluateTestStatistic(RooAbsData& data, RooArgSet& nullPOI) { return EvaluateTestStatistic( data,nullPOI, 0 ); }
159  virtual RooArgList* EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi);
160 
161 
162  virtual TestStatistic* GetTestStatistic(unsigned int i) const {
163  if( fTestStatistics.size() <= i ) return NULL;
164  return fTestStatistics[i];
165  }
166  virtual TestStatistic* GetTestStatistic(void) const { return GetTestStatistic(0); }
167 
168  virtual Double_t ConfidenceLevel() const { return 1. - fSize; }
169  virtual void Initialize(
170  RooAbsArg& /*testStatistic*/,
171  RooArgSet& /*paramsOfInterest*/,
172  RooArgSet& /*nuisanceParameters*/
173  ) {}
174 
175  virtual Int_t GetNToys(void) { return fNToys; }
176  virtual void SetNToys(const Int_t ntoy) { fNToys = ntoy; }
177  virtual void SetNEventsPerToy(const Int_t nevents) {
178  // Forces n events even for extended PDFs. Set NEvents=0 to
179  // use the Poisson distributed events from the extended PDF.
180  fNEvents = nevents;
181  }
182 
183 
184  // Set the Pdf, add to the the workspace if not already there
185  virtual void SetParametersForTestStat(const RooArgSet& nullpoi) {
186  if( fParametersForTestStat ) delete fParametersForTestStat;
187  fParametersForTestStat = (const RooArgSet*)nullpoi.snapshot();
188  }
189 
190  virtual void SetPdf(RooAbsPdf& pdf) { fPdf = &pdf; ClearCache(); }
191 
192  // How to randomize the prior. Set to NULL to deactivate randomization.
193  virtual void SetPriorNuisance(RooAbsPdf* pdf) {
194  fPriorNuisance = pdf;
195  if (fNuisanceParametersSampler) {
196  delete fNuisanceParametersSampler;
197  fNuisanceParametersSampler = NULL;
198  }
199  }
200  // specify the nuisance parameters (eg. the rest of the parameters)
201  virtual void SetNuisanceParameters(const RooArgSet& np) { fNuisancePars = &np; }
202  // specify the observables in the dataset (needed to evaluate the test statistic)
203  virtual void SetObservables(const RooArgSet& o) { fObservables = &o; }
204  // specify the conditional observables
205  virtual void SetGlobalObservables(const RooArgSet& o) { fGlobalObservables = &o; }
206 
207 
208  // set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
209  virtual void SetTestSize(Double_t size) { fSize = size; }
210  // set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
211  virtual void SetConfidenceLevel(Double_t cl) { fSize = 1. - cl; }
212 
213  // Set the TestStatistic (want the argument to be a function of the data & parameter points
214  virtual void SetTestStatistic(TestStatistic *testStatistic, unsigned int i) {
215  if( fTestStatistics.size() < i ) {
216  oocoutE((TObject*)NULL,InputArguments) << "Cannot set test statistic for this index." << std::endl;
217  return;
218  }
219  if( fTestStatistics.size() == i)
220  fTestStatistics.push_back(testStatistic);
221  else
222  fTestStatistics[i] = testStatistic;
223  }
224  virtual void SetTestStatistic(TestStatistic *t) { return SetTestStatistic(t,0); }
225 
226  virtual void SetExpectedNuisancePar(Bool_t i = kTRUE) { fExpectedNuisancePar = i; }
227  virtual void SetAsimovNuisancePar(Bool_t i = kTRUE) { fExpectedNuisancePar = i; }
228 
229  // Checks for sufficient information to do a GetSamplingDistribution(...).
230  Bool_t CheckConfig(void);
231 
232  // control to use bin data generation (=> see RooFit::AllBinned() option)
233  void SetGenerateBinned(bool binned = true) { fGenerateBinned = binned; }
234  // name of the tag for individual components to be generated binned (=> see RooFit::GenBinned() option)
235  void SetGenerateBinnedTag( const char* binnedTag = "" ) { fGenerateBinnedTag = binnedTag; }
236  // set auto binned generation (=> see RooFit::AutoBinned() option)
237  void SetGenerateAutoBinned( Bool_t autoBinned = kTRUE ) { fGenerateAutoBinned = autoBinned; }
238 
239  // Set the name of the sampling distribution used for plotting
240  void SetSamplingDistName(const char* name) { if(name) fSamplingDistName = name; }
241  std::string GetSamplingDistName(void) { return fSamplingDistName; }
242 
243  // This option forces a maximum number of total toys.
244  void SetMaxToys(Double_t t) { fMaxToys = t; }
245 
246  void SetToysLeftTail(Double_t toys, Double_t threshold) {
247  fToysInTails = toys;
248  fAdaptiveLowLimit = threshold;
249  fAdaptiveHighLimit = RooNumber::infinity();
250  }
251  void SetToysRightTail(Double_t toys, Double_t threshold) {
252  fToysInTails = toys;
253  fAdaptiveHighLimit = threshold;
254  fAdaptiveLowLimit = -RooNumber::infinity();
255  }
256  void SetToysBothTails(Double_t toys, Double_t low_threshold, Double_t high_threshold) {
257  fToysInTails = toys;
258  fAdaptiveHighLimit = high_threshold;
259  fAdaptiveLowLimit = low_threshold;
260  }
261 
262  // calling with argument or NULL deactivates proof
263  void SetProofConfig(ProofConfig *pc = NULL) { fProofConfig = pc; }
264 
265  void SetProtoData(const RooDataSet* d) { fProtoData = d; }
266 
267  protected:
268 
269  const RooArgList* EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi, DetailedOutputAggregator& detOutAgg);
270 
271  // helper for GenerateToyData
272  RooAbsData* Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooDataSet *protoData=NULL, int forceEvents=0) const;
273 
274  // helper method for clearing the cache
275  virtual void ClearCache();
276 
277 
278  // densities, snapshots, and test statistics to reweight to
279  RooAbsPdf *fPdf; // model (can be alt or null)
281  std::vector<TestStatistic*> fTestStatistics;
282 
283  std::string fSamplingDistName; // name of the model
284  RooAbsPdf *fPriorNuisance; // prior pdf for nuisance parameters
288  Int_t fNToys; // number of toys to generate
289  Int_t fNEvents; // number of events per toy (may be ignored depending on settings)
291  Bool_t fExpectedNuisancePar; // whether to use expectation values for nuisance parameters (ie Asimov data set)
295 
296  // minimum no of toys in tails for adaptive sampling
297  // (taking weights into account, therefore double)
298  // Default: 0.0 which means no adaptive sampling
300  // maximum no of toys
301  // (taking weights into account, therefore double)
303  // tails
306 
307  const RooDataSet *fProtoData; // in dev
308 
310 
312 
313  // objects below cache information and are mutable and non-persistent
314  mutable RooArgSet* _allVars ; //!
315  mutable std::list<RooAbsPdf*> _pdfList ; //!
316  mutable std::list<RooArgSet*> _obsList ; //!
317  mutable std::list<RooAbsPdf::GenSpec*> _gsList ; //!
318  mutable RooAbsPdf::GenSpec* _gs1 ; //! GenSpec #1
319  mutable RooAbsPdf::GenSpec* _gs2 ; //! GenSpec #2
320  mutable RooAbsPdf::GenSpec* _gs3 ; //! GenSpec #3
321  mutable RooAbsPdf::GenSpec* _gs4 ; //! GenSpec #4
322 
323  static Bool_t fgAlwaysUseMultiGen ; // Use PrepareMultiGen always
324  Bool_t fUseMultiGen ; // Use PrepareMultiGen?
325 
326  protected:
327  ClassDef(ToyMCSampler,3) // A simple implementation of the TestStatSampler interface
328 };
329 }
330 
331 
332 #endif
Holds configuration options for proof and proof-lite.
Definition: ProofConfig.h:49
virtual void SetParametersForTestStat(const RooArgSet &nullpoi)
Definition: ToyMCSampler.h:185
virtual void SetTestStatistic(TestStatistic *t)
Definition: ToyMCSampler.h:224
void SetToysLeftTail(Double_t toys, Double_t threshold)
Definition: ToyMCSampler.h:246
virtual void SetGlobalObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:205
virtual Double_t EvaluateTestStatistic(RooAbsData &data, RooArgSet &nullPOI, int i)
Definition: ToyMCSampler.h:155
TestStatSampler is an interface class for a tools which produce RooStats SamplingDistributions.
virtual Double_t ConfidenceLevel() const
Definition: ToyMCSampler.h:168
#define oocoutI(o, a)
Definition: RooMsgService.h:45
RooAbsPdf::GenSpec * _gs1
Definition: ToyMCSampler.h:318
virtual void SetExpectedNuisancePar(Bool_t i=kTRUE)
Definition: ToyMCSampler.h:226
virtual void SetNEventsPerToy(const Int_t nevents)
Definition: ToyMCSampler.h:177
Basic string class.
Definition: TString.h:137
std::string fSamplingDistName
Definition: ToyMCSampler.h:283
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
NuisanceParametersSampler * fNuisanceParametersSampler
Definition: ToyMCSampler.h:311
virtual void SetConfidenceLevel(Double_t cl)
Definition: ToyMCSampler.h:211
virtual RooAbsData * GenerateToyData(RooArgSet &paramPoint) const
Definition: ToyMCSampler.h:145
virtual TestStatistic * GetTestStatistic(unsigned int i) const
Definition: ToyMCSampler.h:162
virtual void SetTestStatistic(TestStatistic *testStatistic, unsigned int i)
Definition: ToyMCSampler.h:214
virtual RooAbsData * GenerateToyData(RooArgSet &paramPoint, double &weight) const
Definition: ToyMCSampler.h:148
void NextPoint(RooArgSet &nuisPoint, Double_t &weight)
#define ClassDef(name, id)
Definition: Rtypes.h:254
void SetToysRightTail(Double_t toys, Double_t threshold)
Definition: ToyMCSampler.h:251
void SetToysBothTails(Double_t toys, Double_t low_threshold, Double_t high_threshold)
Definition: ToyMCSampler.h:256
#define oocoutE(o, a)
Definition: RooMsgService.h:48
void SetGenerateAutoBinned(Bool_t autoBinned=kTRUE)
Definition: ToyMCSampler.h:237
NuisanceParametersSampler(RooAbsPdf *prior=NULL, const RooArgSet *parameters=NULL, Int_t nToys=1000, Bool_t asimov=kFALSE)
Definition: ToyMCSampler.h:54
std::list< RooArgSet * > _obsList
Definition: ToyMCSampler.h:316
virtual void SetNToys(const Int_t ntoy)
Definition: ToyMCSampler.h:176
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
void SetGenerateBinnedTag(const char *binnedTag="")
Definition: ToyMCSampler.h:235
virtual void Initialize(RooAbsArg &, RooArgSet &, RooArgSet &)
Definition: ToyMCSampler.h:169
virtual Int_t GetNToys(void)
Definition: ToyMCSampler.h:175
const RooDataSet * fProtoData
Definition: ToyMCSampler.h:307
void SetProofConfig(ProofConfig *pc=NULL)
Definition: ToyMCSampler.h:263
ProofConfig * fProofConfig
Definition: ToyMCSampler.h:309
static Double_t infinity()
Return internal infinity representation.
Definition: RooNumber.cxx:49
std::vector< TestStatistic * > fTestStatistics
Definition: ToyMCSampler.h:281
const RooArgSet * fGlobalObservables
Definition: ToyMCSampler.h:287
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:99
void SetSamplingDistName(const char *name)
Definition: ToyMCSampler.h:240
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
This class simply holds a sampling distribution of some test statistic.
virtual void SetAsimovNuisancePar(Bool_t i=kTRUE)
Definition: ToyMCSampler.h:227
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooAbsPdf::GenSpec * _gs4
GenSpec #3.
Definition: ToyMCSampler.h:321
double Double_t
Definition: RtypesCore.h:55
virtual void SetObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:203
void SetMaxToys(Double_t t)
Definition: ToyMCSampler.h:244
void SetProtoData(const RooDataSet *d)
Definition: ToyMCSampler.h:265
std::list< RooAbsPdf::GenSpec * > _gsList
Definition: ToyMCSampler.h:317
virtual void SetPdf(RooAbsPdf &pdf)
Definition: ToyMCSampler.h:190
static Bool_t fgAlwaysUseMultiGen
GenSpec #4.
Definition: ToyMCSampler.h:323
virtual RooAbsData * GenerateToyData(RooArgSet &paramPoint, RooAbsPdf &pdf) const
Definition: ToyMCSampler.h:140
virtual void SetPriorNuisance(RooAbsPdf *pdf)
Definition: ToyMCSampler.h:193
Mother of all ROOT objects.
Definition: TObject.h:37
const RooArgSet * fObservables
Definition: ToyMCSampler.h:286
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
const RooArgSet * fParametersForTestStat
Definition: ToyMCSampler.h:280
#define NULL
Definition: Rtypes.h:82
virtual void SetTestSize(Double_t size)
Definition: ToyMCSampler.h:209
void SetUseMultiGen(Bool_t flag)
Definition: ToyMCSampler.h:109
RooAbsPdf::GenSpec * _gs3
GenSpec #2.
Definition: ToyMCSampler.h:320
void SetGenerateBinned(bool binned=true)
Definition: ToyMCSampler.h:233
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
RooAbsPdf::GenSpec * _gs2
GenSpec #1.
Definition: ToyMCSampler.h:319
std::string GetSamplingDistName(void)
Definition: ToyMCSampler.h:241
const Bool_t kTRUE
Definition: Rtypes.h:91
const RooArgSet * fNuisancePars
Definition: ToyMCSampler.h:285
virtual void SetNuisanceParameters(const RooArgSet &np)
Definition: ToyMCSampler.h:201
TestStatistic is an interface class to provide a facility for construction test statistics distributi...
Definition: TestStatistic.h:33
virtual void AddTestStatistic(TestStatistic *t=NULL)
Definition: ToyMCSampler.h:125
Helper class for ToyMCSampler.
Definition: ToyMCSampler.h:51
RooAbsPdf * fPriorNuisance
Definition: ToyMCSampler.h:284
std::list< RooAbsPdf * > _pdfList
Definition: ToyMCSampler.h:315
char name[80]
Definition: TGX11.cxx:109
virtual TestStatistic * GetTestStatistic(void) const
Definition: ToyMCSampler.h:166
virtual Double_t EvaluateTestStatistic(RooAbsData &data, RooArgSet &nullPOI)
Definition: ToyMCSampler.h:158