ROOT  6.06/09
Reference Guide
ToyMCSampler.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Sven Kreiss June 2010
3 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
4 /*************************************************************************
5  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #include "RooStats/ToyMCSampler.h"
13 
14 #ifndef ROO_MSG_SERVICE
15 #include "RooMsgService.h"
16 #endif
17 
18 #ifndef ROO_DATA_HIST
19 #include "RooDataHist.h"
20 #endif
21 
22 #ifndef ROO_REAL_VAR
23 #include "RooRealVar.h"
24 #endif
25 
26 #include "TCanvas.h"
27 #include "RooPlot.h"
28 #include "RooRandom.h"
29 
30 #include "RooStudyManager.h"
31 #include "RooStats/ToyMCStudy.h"
33 #include "RooStats/RooStatsUtils.h"
34 #include "RooSimultaneous.h"
35 #include "RooCategory.h"
36 
37 #include "TMath.h"
38 
39 
40 using namespace RooFit;
41 using namespace std;
42 
43 
45 
46 namespace RooStats {
47 
48 void NuisanceParametersSampler::NextPoint(RooArgSet& nuisPoint, Double_t& weight) {
49  // Assigns new nuisance parameter point to members of nuisPoint.
50  // nuisPoint can be more objects than just the nuisance
51  // parameters.
52 
53  // check whether to get new set of nuisanceParPoints
54  if (fIndex >= fNToys) {
55  Refresh();
56  fIndex = 0;
57  }
58 
59  // get value
60  nuisPoint = *fPoints->get(fIndex++);
61  weight = fPoints->weight();
62 
63  // check whether result will have any influence
64  if(fPoints->weight() == 0.0) {
65  oocoutI((TObject*)NULL,Generation) << "Weight 0 encountered. Skipping." << endl;
66  NextPoint(nuisPoint, weight);
67  }
68 }
69 void NuisanceParametersSampler::Refresh() {
70  // Creates the initial set of nuisance parameter points. It also refills the
71  // set with new parameter points if called repeatedly. This helps with
72  // adaptive sampling as the required number of nuisance parameter points
73  // might increase during the run.
74 
75  if (!fPrior || !fParams) return;
76 
77  if (fPoints) delete fPoints;
78 
79  if (fExpected) {
80  // UNDER CONSTRUCTION
81  oocoutI((TObject*)NULL,InputArguments) << "Using expected nuisance parameters." << endl;
82 
83  int nBins = fNToys;
84 
85  // From FeldmanCousins.cxx:
86  // set nbins for the POI
87  TIter it2 = fParams->createIterator();
88  RooRealVar *myarg2;
89  while ((myarg2 = dynamic_cast<RooRealVar*>(it2.Next()))) {
90  myarg2->setBins(nBins);
91  }
92 
93 
94  fPoints = fPrior->generate(
95  *fParams,
96  AllBinned(),
97  ExpectedData(),
98  NumEvents(1) // for Asimov set, this is only a scale factor
99  );
100  if(fPoints->numEntries() != fNToys) {
101  fNToys = fPoints->numEntries();
103  "Adjusted number of toys to number of bins of nuisance parameters: " << fNToys << endl;
104  }
105 
106 /*
107  // check
108  TCanvas *c1 = new TCanvas;
109  RooPlot *p = dynamic_cast<RooRealVar*>(fParams->first())->frame();
110  fPoints->plotOn(p);
111  p->Draw();
112  for(int x=0; x < fPoints->numEntries(); x++) {
113  fPoints->get(x)->Print("v");
114  cout << fPoints->weight() << endl;
115  }
116 */
117 
118  }else{
119  oocoutI((TObject*)NULL,InputArguments) << "Using randomized nuisance parameters." << endl;
120 
121  fPoints = fPrior->generate(*fParams, fNToys);
122  }
123 }
124 
125 
126 
127 Bool_t ToyMCSampler::fgAlwaysUseMultiGen = kFALSE ;
128 
129 
130 void ToyMCSampler::SetAlwaysUseMultiGen(Bool_t flag) { fgAlwaysUseMultiGen = flag ; }
131 
132 
133 
134 ToyMCSampler::ToyMCSampler() : fSamplingDistName("SD"), fNToys(1)
135 {
136  // Proof constructor. Do not use.
137 
138  fPdf = NULL;
139  fParametersForTestStat = NULL;
140  fPriorNuisance = NULL;
141  fNuisancePars = NULL;
142  fObservables = NULL;
143  fGlobalObservables = NULL;
144 
145  fSize = 0.05;
146  fNEvents = 0;
147  fGenerateBinned = kFALSE;
148  fGenerateBinnedTag = "";
149  fGenerateAutoBinned = kTRUE;
150  fExpectedNuisancePar = kFALSE;
151 
152  fToysInTails = 0.0;
153  fMaxToys = RooNumber::infinity();
154  fAdaptiveLowLimit = -RooNumber::infinity();
155  fAdaptiveHighLimit = RooNumber::infinity();
156 
157  fProtoData = NULL;
158 
159  fProofConfig = NULL;
160  fNuisanceParametersSampler = NULL;
161 
162  _allVars = NULL ;
163  _gs1 = NULL ;
164  _gs2 = NULL ;
165  _gs3 = NULL ;
166  _gs4 = NULL ;
167 
168  //suppress messages for num integration of Roofit
170 
171  fUseMultiGen = kFALSE ;
172 }
173 
174 ToyMCSampler::ToyMCSampler(TestStatistic &ts, Int_t ntoys) :
175  fSamplingDistName(ts.GetVarName()), fNToys(ntoys)
176 {
177  fPdf = NULL;
178  fParametersForTestStat = NULL;
179  fPriorNuisance = NULL;
180  fNuisancePars = NULL;
181  fObservables = NULL;
182  fGlobalObservables = NULL;
183 
184  fSize = 0.05;
185  fNEvents = 0;
186  fGenerateBinned = kFALSE;
187  fGenerateBinnedTag = "";
188  fGenerateAutoBinned = kTRUE;
189  fExpectedNuisancePar = kFALSE;
190 
191  fToysInTails = 0.0;
192  fMaxToys = RooNumber::infinity();
193  fAdaptiveLowLimit = -RooNumber::infinity();
194  fAdaptiveHighLimit = RooNumber::infinity();
195 
196  fProtoData = NULL;
197 
198  fProofConfig = NULL;
199  fNuisanceParametersSampler = NULL;
200 
201  _allVars = NULL ;
202  _gs1 = NULL ;
203  _gs2 = NULL ;
204  _gs3 = NULL ;
205  _gs4 = NULL ;
206 
207  //suppress messages for num integration of Roofit
209 
210  fUseMultiGen = kFALSE ;
211 
212  AddTestStatistic(&ts);
213 }
214 
215 
216 ToyMCSampler::~ToyMCSampler() {
217  if(fNuisanceParametersSampler) delete fNuisanceParametersSampler;
218 
219  ClearCache();
220 }
221 
222 
223 Bool_t ToyMCSampler::CheckConfig(void) {
224  // only checks, no guessing/determination (do this in calculators,
225  // e.g. using ModelConfig::GuessObsAndNuisance(...))
226  bool goodConfig = true;
227 
228  if(fTestStatistics.size() == 0 || fTestStatistics[0] == NULL) { ooccoutE((TObject*)NULL,InputArguments) << "Test statistic not set." << endl; goodConfig = false; }
229  if(!fObservables) { ooccoutE((TObject*)NULL,InputArguments) << "Observables not set." << endl; goodConfig = false; }
230  if(!fParametersForTestStat) { ooccoutE((TObject*)NULL,InputArguments) << "Parameter values used to evaluate the test statistic are not set." << endl; goodConfig = false; }
231  if(!fPdf) { ooccoutE((TObject*)NULL,InputArguments) << "Pdf not set." << endl; goodConfig = false; }
232 
233 
234  //ooccoutI((TObject*)NULL,InputArguments) << "ToyMCSampler configuration:" << endl;
235  //ooccoutI((TObject*)NULL,InputArguments) << "Pdf from SetPdf: " << fPdf << endl;
236  // for( unsigned int i=0; i < fTestStatistics.size(); i++ ) {
237  // ooccoutI((TObject*)NULL,InputArguments) << "test statistics["<<i<<"]: " << fTestStatistics[i] << endl;
238  // }
239  //ooccoutI((TObject*)NULL,InputArguments) << endl;
240 
241  return goodConfig;
242 }
243 
244 
245 RooArgList* ToyMCSampler::EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi) {
246  // Evaluate all test statistics, returning result and any detailed output.
247  // PDF parameter values are saved in case they are modified by
248  // TestStatistic::Evaluate (eg. SimpleLikelihoodRatioTestStat).
249  DetailedOutputAggregator detOutAgg;
250  const RooArgList* allTS = EvaluateAllTestStatistics(data, poi, detOutAgg);
251  if (!allTS) return 0;
252  // no need to delete allTS, it is deleted in destructor of detOutAgg
253  return dynamic_cast<RooArgList*>(allTS->snapshot());
254 }
255 
256 
257 const RooArgList* ToyMCSampler::EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi, DetailedOutputAggregator& detOutAgg) {
258  RooArgSet *allVars = fPdf ? fPdf->getVariables() : 0;
259  RooArgSet *saveAll = allVars ? dynamic_cast<RooArgSet*>(allVars->snapshot()) : 0;
260  for( unsigned int i = 0; i < fTestStatistics.size(); i++ ) {
261  if( fTestStatistics[i] == NULL ) continue;
262  TString name( TString::Format("%s_TS%u", fSamplingDistName.c_str(), i) );
263  RooArgSet* parForTS = dynamic_cast<RooArgSet*>(poi.snapshot());
264  RooRealVar ts( name, fTestStatistics[i]->GetVarName(), fTestStatistics[i]->Evaluate( data, *parForTS ) );
265  RooArgList tset(ts);
266  detOutAgg.AppendArgSet(&tset);
267  delete parForTS;
268  if (const RooArgSet* detOut = fTestStatistics[i]->GetDetailedOutput()) {
269  name.Append("_");
270  detOutAgg.AppendArgSet(detOut, name);
271  }
272  if (saveAll) *allVars = *saveAll; // restore values, perhaps modified by fTestStatistics[i]->Evaluate()
273  }
274  delete saveAll;
275  delete allVars;
276  return detOutAgg.GetAsArgList();
277 }
278 
279 
280 SamplingDistribution* ToyMCSampler::GetSamplingDistribution(RooArgSet& paramPointIn) {
281  if(fTestStatistics.size() > 1) {
282  oocoutW((TObject*)NULL, InputArguments) << "Multiple test statistics defined, but only one distribution will be returned." << endl;
283  for( unsigned int i=0; i < fTestStatistics.size(); i++ ) {
284  oocoutW((TObject*)NULL, InputArguments) << " \t test statistic: " << fTestStatistics[i] << endl;
285  }
286  }
287 
288  RooDataSet* r = GetSamplingDistributions(paramPointIn);
289  if(r == NULL || r->numEntries() == 0) {
290  oocoutW((TObject*)NULL, Generation) << "no sampling distribution generated" << endl;
291  return NULL;
292  }
293 
294  SamplingDistribution* samp = new SamplingDistribution( r->GetName(), r->GetTitle(), *r );
295  delete r;
296  return samp;
297 }
298 
299 
300 RooDataSet* ToyMCSampler::GetSamplingDistributions(RooArgSet& paramPointIn)
301 {
302  // Use for serial and parallel runs.
303 
304  // ======= S I N G L E R U N ? =======
305  if(!fProofConfig)
306  return GetSamplingDistributionsSingleWorker(paramPointIn);
307 
308 
309  // ======= P A R A L L E L R U N =======
310  if (!CheckConfig()){
312  << "Bad COnfiguration in ToyMCSampler "
313  << endl;
314  return nullptr;
315  }
316 
317  // turn adaptive sampling off if given
318  if(fToysInTails) {
319  fToysInTails = 0;
321  << "Adaptive sampling in ToyMCSampler is not supported for parallel runs."
322  << endl;
323  }
324 
325  // adjust number of toys on the slaves to keep the total number of toys constant
326  Int_t totToys = fNToys;
327  fNToys = (int)ceil((double)fNToys / (double)fProofConfig->GetNExperiments()); // round up
328 
329  // create the study instance for parallel processing
330  ToyMCStudy* toymcstudy = new ToyMCStudy ;
331  toymcstudy->SetToyMCSampler(*this);
332  toymcstudy->SetParamPoint(paramPointIn);
333  toymcstudy->SetRandomSeed(RooRandom::randomGenerator()->Integer(TMath::Limits<unsigned int>::Max() ) );
334 
335  // temporary workspace for proof to avoid messing with TRef
336  RooWorkspace w(fProofConfig->GetWorkspace());
337  RooStudyManager studymanager(w, *toymcstudy);
338  studymanager.runProof(fProofConfig->GetNExperiments(), fProofConfig->GetHost(), fProofConfig->GetShowGui());
339 
340  RooDataSet* output = toymcstudy->merge();
341 
342  // reset the number of toys
343  fNToys = totToys;
344 
345  delete toymcstudy;
346  return output;
347 }
348 
349 RooDataSet* ToyMCSampler::GetSamplingDistributionsSingleWorker(RooArgSet& paramPointIn)
350 {
351  // This is the main function for serial runs. It is called automatically
352  // from inside GetSamplingDistribution when no ProofConfig is given.
353  // You should not call this function yourself. This function should
354  // be used by ToyMCStudy on the workers (ie. when you explicitly want
355  // a serial run although ProofConfig is present).
356 
357  // Make sure the cache is clear. It is important to clear it hear, because
358  // the cache might be invalid even when just the firstPOI was changed, for which
359  // no accessor has to be called. (Fixes a bug when ToyMCSampler is
360  // used with the Neyman Construction)
361  ClearCache();
362 
363  if (!CheckConfig()){
365  << "Bad COnfiguration in ToyMCSampler "
366  << endl;
367  return nullptr;
368  }
369 
370 
371  // important to cache the paramPoint b/c test statistic might
372  // modify it from event to event
373  RooArgSet *paramPoint = (RooArgSet*) paramPointIn.snapshot();
374  RooArgSet *allVars = fPdf->getVariables();
375  RooArgSet *saveAll = (RooArgSet*) allVars->snapshot();
376 
377 
378  DetailedOutputAggregator detOutAgg;
379 
380  // counts the number of toys in the limits set for adaptive sampling
381  // (taking weights into account; always on first test statistic)
382  Double_t toysInTails = 0.0;
383 
384  for (Int_t i = 0; i < fMaxToys; ++i) {
385  // need to check at the beginning for case that zero toys are requested
386  if (toysInTails >= fToysInTails && i+1 > fNToys) break;
387 
388  // status update
389  if ( i% 500 == 0 && i>0 ) {
390  oocoutP((TObject*)0,Generation) << "generated toys: " << i << " / " << fNToys;
391  if (fToysInTails) ooccoutP((TObject*)0,Generation) << " (tails: " << toysInTails << " / " << fToysInTails << ")" << std::endl;
392  else ooccoutP((TObject*)0,Generation) << endl;
393  }
394 
395 
396  // TODO: change this treatment to keep track of all values so that the threshold
397  // for adaptive sampling is counted for all distributions and not just the
398  // first one.
399  Double_t valueFirst = -999.0, weight = 1.0;
400 
401  // set variables to requested parameter point
402  *allVars = *saveAll; // important for example for SimpleLikelihoodRatioTestStat
403 
404  RooAbsData* toydata = GenerateToyData(*paramPoint, weight);
405 
406  *allVars = *fParametersForTestStat;
407 
408  const RooArgList* allTS = EvaluateAllTestStatistics(*toydata, *fParametersForTestStat, detOutAgg);
409  if (allTS->getSize() > Int_t(fTestStatistics.size()))
410  detOutAgg.AppendArgSet( fGlobalObservables, "globObs_" );
411  if (RooRealVar* firstTS = dynamic_cast<RooRealVar*>(allTS->first()))
412  valueFirst = firstTS->getVal();
413 
414  delete toydata;
415 
416  // check for nan
417  if(valueFirst != valueFirst) {
418  oocoutW((TObject*)NULL, Generation) << "skip: " << valueFirst << ", " << weight << endl;
419  continue;
420  }
421 
422  detOutAgg.CommitSet(weight);
423 
424  // adaptive sampling checks
425  if (valueFirst <= fAdaptiveLowLimit || valueFirst >= fAdaptiveHighLimit) {
426  if(weight >= 0.) toysInTails += weight;
427  else toysInTails += 1.;
428  }
429  }
430 
431  // clean up
432  *allVars = *saveAll;
433  delete saveAll;
434  delete allVars;
435  delete paramPoint;
436 
437  return detOutAgg.GetAsDataSet(fSamplingDistName, fSamplingDistName);
438 }
439 
440 void ToyMCSampler::GenerateGlobalObservables(RooAbsPdf& pdf) const {
441 
442 
443  if(!fGlobalObservables || fGlobalObservables->getSize()==0) {
444  ooccoutE((TObject*)NULL,InputArguments) << "Global Observables not set." << endl;
445  return;
446  }
447 
448 
449  if (fUseMultiGen || fgAlwaysUseMultiGen) {
450 
451  // generate one set of global observables and assign it
452  // has problem for sim pdfs
453  RooSimultaneous* simPdf = dynamic_cast<RooSimultaneous*>( &pdf );
454  if (!simPdf) {
455  RooDataSet *one = pdf.generate(*fGlobalObservables, 1);
456 
457  const RooArgSet *values = one->get(0);
458  if (!_allVars) {
459  _allVars = pdf.getVariables();
460  }
461  *_allVars = *values;
462  delete one;
463 
464  } else {
465 
466  if (_pdfList.size() == 0) {
467  RooCategory& channelCat = (RooCategory&)simPdf->indexCat();
468  int nCat = channelCat.numTypes();
469  for (int i=0; i < nCat; ++i){
470  channelCat.setIndex(i);
471  RooAbsPdf* pdftmp = simPdf->getPdf(channelCat.getLabel());
472  assert(pdftmp);
473  RooArgSet* globtmp = pdftmp->getObservables(*fGlobalObservables);
474  RooAbsPdf::GenSpec* gs = pdftmp->prepareMultiGen(*globtmp, NumEvents(1));
475  _pdfList.push_back(pdftmp);
476  _obsList.push_back(globtmp);
477  _gsList.push_back(gs);
478  }
479  }
480 
481  list<RooArgSet*>::iterator oiter = _obsList.begin();
482  list<RooAbsPdf::GenSpec*>::iterator giter = _gsList.begin();
483  for (list<RooAbsPdf*>::iterator iter = _pdfList.begin(); iter != _pdfList.end(); ++iter, ++giter, ++oiter) {
484  //RooDataSet* tmp = (*iter)->generate(**oiter,1) ;
485  RooDataSet* tmp = (*iter)->generate(**giter);
486  **oiter = *tmp->get(0);
487  delete tmp;
488  }
489  }
490 
491 
492  } else {
493 
494  // not using multigen for global observables
495  RooDataSet* one = pdf.generateSimGlobal( *fGlobalObservables, 1 );
496  const RooArgSet *values = one->get(0);
497  RooArgSet* allVars = pdf.getVariables();
498  *allVars = *values;
499  delete allVars;
500  delete one;
501 
502  }
503 }
504 
505 RooAbsData* ToyMCSampler::GenerateToyData(RooArgSet& paramPoint, double& weight, RooAbsPdf& pdf) const {
506  // This method generates a toy data set for the given parameter point taking
507  // global observables into account.
508  // The values of the generated global observables remain in the pdf's variables.
509  // They have to have those values for the subsequent evaluation of the
510  // test statistics.
511 
512  if(!fObservables) {
513  ooccoutE((TObject*)NULL,InputArguments) << "Observables not set." << endl;
514  return NULL;
515  }
516 
517  // assign input paramPoint
518  RooArgSet* allVars = fPdf->getVariables();
519  *allVars = paramPoint;
520 
521 
522  // create nuisance parameter points
523  if(!fNuisanceParametersSampler && fPriorNuisance && fNuisancePars) {
524  fNuisanceParametersSampler = new NuisanceParametersSampler(fPriorNuisance, fNuisancePars, fNToys, fExpectedNuisancePar);
525  if ((fUseMultiGen || fgAlwaysUseMultiGen) && fNuisanceParametersSampler )
526  oocoutI((TObject*)NULL,InputArguments) << "Cannot use multigen when nuisance parameters vary for every toy" << endl;
527  }
528 
529  // generate global observables
530  RooArgSet observables(*fObservables);
531  if(fGlobalObservables && fGlobalObservables->getSize()) {
532  observables.remove(*fGlobalObservables);
533  GenerateGlobalObservables(pdf);
534  }
535 
536  // save values to restore later.
537  // but this must remain after(!) generating global observables
538  const RooArgSet* saveVars = (const RooArgSet*)allVars->snapshot();
539 
540  if(fNuisanceParametersSampler) { // use nuisance parameters?
541  // Construct a set of nuisance parameters that has the parameters
542  // in the input paramPoint removed. Therefore, no parameter in
543  // paramPoint is randomized.
544  // Therefore when a parameter is given (should be held fixed),
545  // but is also in the list of nuisance parameters, the parameter
546  // will be held fixed. This is useful for debugging to hold single
547  // parameters fixed although under "normal" circumstances it is
548  // randomized.
549  RooArgSet allVarsMinusParamPoint(*allVars);
550  allVarsMinusParamPoint.remove(paramPoint, kFALSE, kTRUE); // match by name
551 
552  // get nuisance parameter point and weight
553  fNuisanceParametersSampler->NextPoint(allVarsMinusParamPoint, weight);
554 
555 
556  }else{
557  weight = 1.0;
558  }
559 
560  RooAbsData *data = Generate(pdf, observables);
561 
562  // We generated the data with the randomized nuisance parameter (if hybrid)
563  // but now we are setting the nuisance parameters back to where they were.
564  *allVars = *saveVars;
565  delete allVars;
566  delete saveVars;
567 
568  return data;
569 }
570 
571 
572 
573 RooAbsData* ToyMCSampler::Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooDataSet* protoData, int forceEvents) const {
574  // This is the generate function to use in the context of the ToyMCSampler
575  // instead of the standard RooAbsPdf::generate(...).
576  // It takes into account whether the number of events is given explicitly
577  // or whether it should use the expected number of events. It also takes
578  // into account the option to generate a binned data set (ie RooDataHist).
579 
580  if(fProtoData) {
581  protoData = fProtoData;
582  forceEvents = protoData->numEntries();
583  }
584 
585  RooAbsData *data = NULL;
586  int events = forceEvents;
587  if(events == 0) events = fNEvents;
588 
589  // cannot use multigen when the nuisance parameters change for every toy
590  bool useMultiGen = (fUseMultiGen || fgAlwaysUseMultiGen) && !fNuisanceParametersSampler;
591 
592  if(events == 0) {
593  if( pdf.canBeExtended() && pdf.expectedEvents(observables) > 0) {
594  if(fGenerateBinned) {
595  if(protoData) data = pdf.generate(observables, AllBinned(), Extended(), ProtoData(*protoData, true, true));
596  else data = pdf.generate(observables, AllBinned(), Extended());
597  }else{
598  if(protoData) {
599  if (useMultiGen) {
600  if (!_gs2) { _gs2 = pdf.prepareMultiGen(observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true)) ; }
601  data = pdf.generate(*_gs2) ;
602  } else {
603  data = pdf.generate (observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true));
604  }
605  }
606  else {
607  if (useMultiGen) {
608  if (!_gs1) { _gs1 = pdf.prepareMultiGen(observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag) ) ; }
609  data = pdf.generate(*_gs1) ;
610  } else {
611  data = pdf.generate (observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag) );
612  }
613 
614  }
615  }
616  }else{
618  << "ToyMCSampler: Error : pdf is not extended and number of events per toy is zero"
619  << endl;
620  }
621  }else{
622  if(fGenerateBinned) {
623  if(protoData) data = pdf.generate(observables, events, AllBinned(), ProtoData(*protoData, true, true));
624  else data = pdf.generate(observables, events, AllBinned());
625  }else{
626  if(protoData) {
627  if (useMultiGen) {
628  if (!_gs3) { _gs3 = pdf.prepareMultiGen(observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true)); }
629  data = pdf.generate(*_gs3) ;
630  } else {
631  data = pdf.generate (observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true));
632  }
633  } else {
634  if (useMultiGen) {
635  if (!_gs4) { _gs4 = pdf.prepareMultiGen(observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag)); }
636  data = pdf.generate(*_gs4) ;
637  } else {
638  data = pdf.generate (observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag));
639  }
640  }
641  }
642  }
643 
644  // in case of number counting print observables
645  // if (data->numEntries() == 1) {
646  // std::cout << "generate observables : ";
647  // RooStats::PrintListContent(*data->get(0), std::cout);
648  // }
649 
650  return data;
651 }
652 
653 
654 
655 
656 
657 // Extended interface to append to sampling distribution more samples
658 SamplingDistribution* ToyMCSampler::AppendSamplingDistribution(
659  RooArgSet& allParameters,
660  SamplingDistribution* last,
661  Int_t additionalMC)
662 {
663  Int_t tmp = fNToys;
664  fNToys = additionalMC;
665  SamplingDistribution* newSamples = GetSamplingDistribution(allParameters);
666  fNToys = tmp;
667 
668  if(last){
669  last->Add(newSamples);
670  delete newSamples;
671  return last;
672  }
673 
674  return newSamples;
675 }
676 
677 
678 
679 
680 void ToyMCSampler::ClearCache() {
681  // clear the cache obtained from the pdf used for speeding the toy and global observables generation
682  // needs to be called every time the model pdf (fPdf) changes
683 
684 
685  if (_gs1) delete _gs1;
686  _gs1 = 0;
687  if (_gs2) delete _gs2;
688  _gs2 = 0;
689  if (_gs3) delete _gs3;
690  _gs3 = 0;
691  if (_gs4) delete _gs4;
692  _gs4 = 0;
693 
694  // no need to delete the _pdfList since it is managed by the RooSimultaneous object
695  if (_pdfList.size() > 0) {
696  std::list<RooArgSet*>::iterator oiter = _obsList.begin();
697  for (std::list<RooAbsPdf::GenSpec*>::iterator giter = _gsList.begin(); giter != _gsList.end(); ++giter, ++oiter) {
698  delete *oiter;
699  delete *giter;
700  }
701  _pdfList.clear();
702  _obsList.clear();
703  _gsList.clear();
704  }
705 
706  //LM: is this set really needed ??
707  if (_allVars) delete _allVars;
708  _allVars = 0;
709 
710 }
711 
712 } // end namespace RooStats
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:52
RooCmdArg AllBinned()
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
virtual Bool_t setIndex(Int_t index, Bool_t printError=kTRUE)
Set value by specifying the index code of the desired state.
#define assert(cond)
Definition: unittest.h:542
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Definition: RooAbsArg.h:194
const RooAbsCategoryLValue & indexCat() const
#define oocoutI(o, a)
Definition: RooMsgService.h:45
RooCmdArg NumEvents(Int_t numEvents)
Basic string class.
Definition: TString.h:137
int Int_t
Definition: RtypesCore.h:41
StreamConfig & getStream(Int_t id)
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
GenSpec * prepareMultiGen(const RooArgSet &whatVars, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none())
Prepare GenSpec configuration object for efficient generation of multiple datasets from idetical spec...
Definition: RooAbsPdf.cxx:1855
void removeTopic(RooFit::MsgTopic oldTopic)
static RooMsgService & instance()
Return reference to singleton instance.
#define ooccoutP(o, a)
Definition: RooMsgService.h:53
STL namespace.
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events from this p.d.f for use in extended likelihood calculations.
Definition: RooAbsPdf.cxx:2930
RooCmdArg Extended(Bool_t flag=kTRUE)
#define ooccoutE(o, a)
Definition: RooMsgService.h:55
RooAbsArg * first() const
#define oocoutP(o, a)
Definition: RooMsgService.h:46
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
Definition: TString.cxx:2334
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:2081
#define oocoutE(o, a)
Definition: RooMsgService.h:48
Int_t numTypes(const char *=0) const
std::map< std::string, std::string >::const_iterator iter
Definition: TAlienJob.cxx:54
TString & Append(const char *cs)
Definition: TString.h:492
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:53
void setBins(Int_t nBins, const char *name=0)
Definition: RooRealVar.h:78
RooCmdArg ExpectedData(Bool_t flag=kTRUE)
ROOT::R::TRInterface & r
Definition: Object.C:4
TObject * Next()
Definition: TCollection.h:158
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:291
static Double_t infinity()
Return internal infinity representation.
Definition: RooNumber.cxx:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:99
RooCmdArg GenBinned(const char *tag)
Namespace for the RooStats classes.
Definition: Asimov.h:20
double Double_t
Definition: RtypesCore.h:55
RooCmdArg AutoBinned(Bool_t flag=kTRUE)
RooCmdArg ProtoData(const RooDataSet &protoData, Bool_t randomizeOrder=kFALSE, Bool_t resample=kFALSE)
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:216
#define name(a, b)
Definition: linkTestLib0.cpp:5
#define oocoutW(o, a)
Definition: RooMsgService.h:47
Mother of all ROOT objects.
Definition: TObject.h:58
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
RooDataSet * generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none())
Generate a new dataset containing the specified variables with events sampled from our distribution...
Definition: RooAbsPdf.cxx:1702
double ceil(double)
#define NULL
Definition: Rtypes.h:82
ClassImp(RooStats::ToyMCSampler) namespace RooStats
Int_t getSize() const
static void output(int code)
Definition: gifencode.c:226
const Bool_t kTRUE
Definition: Rtypes.h:91
virtual RooDataSet * generateSimGlobal(const RooArgSet &whatVars, Int_t nEvents)
Special generator interface for generation of 'global observables' – for RooStats tools...
Definition: RooAbsPdf.cxx:2383
virtual const RooArgSet * get(Int_t index) const
Return RooArgSet with coordinates of event 'index'.
virtual const char * getLabel() const
Return label string of current state.
Definition: RooCategory.h:40