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